669 namespace MappingFEFieldImplementation
679 template <
int dim,
int spacedim,
typename VectorType>
681 maybe_compute_q_points(
682 const typename ::QProjector<dim>::DataSetDescriptor data_set,
683 const typename ::MappingFEField<dim, spacedim, VectorType>::
687 const std::vector<unsigned int> &fe_to_real,
688 std::vector<Point<spacedim>> &quadrature_points)
694 for (
unsigned int point = 0; point < quadrature_points.size();
698 const double *shape = &data.shape(point + data_set, 0);
700 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
702 const unsigned int comp_k =
703 fe.system_to_component_index(k).first;
705 result[fe_to_real[comp_k]] +=
706 data.local_dof_values[k] * shape[k];
709 quadrature_points[point] = result;
721 template <
int dim,
int spacedim,
typename VectorType>
723 maybe_update_Jacobians(
724 const typename ::QProjector<dim>::DataSetDescriptor data_set,
725 const typename ::MappingFEField<dim, spacedim, VectorType>::
729 const std::vector<unsigned int> &fe_to_real)
736 const unsigned int n_q_points = data.contravariant.size();
740 for (
unsigned int point = 0; point < n_q_points; ++point)
743 &data.derivative(point + data_set, 0);
747 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
749 const unsigned int comp_k =
750 fe.system_to_component_index(k).first;
752 result[fe_to_real[comp_k]] +=
753 data.local_dof_values[k] * data_derv[k];
757 for (
unsigned int i = 0; i < spacedim; ++i)
759 data.contravariant[point][i] = result[i];
767 for (
unsigned int point = 0; point < data.contravariant.size();
769 data.covariant[point] =
770 (data.contravariant[point]).covariant_form();
776 data.volume_elements.size());
777 for (
unsigned int point = 0; point < data.contravariant.size();
779 data.volume_elements[point] =
780 data.contravariant[point].determinant();
790 template <
int dim,
int spacedim,
typename VectorType>
792 maybe_update_jacobian_grads(
793 const typename ::QProjector<dim>::DataSetDescriptor data_set,
794 const typename ::MappingFEField<dim, spacedim, VectorType>::
798 const std::vector<unsigned int> &fe_to_real,
799 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
804 const unsigned int n_q_points = jacobian_grads.size();
806 for (
unsigned int point = 0; point < n_q_points; ++point)
809 &data.second_derivative(point + data_set, 0);
813 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
815 const unsigned int comp_k =
816 fe.system_to_component_index(k).first;
818 for (
unsigned int j = 0; j < dim; ++j)
819 for (
unsigned int l = 0; l < dim; ++l)
820 result[fe_to_real[comp_k]][j][l] +=
821 (
second[k][j][l] * data.local_dof_values[k]);
826 for (
unsigned int i = 0; i < spacedim; ++i)
827 for (
unsigned int j = 0; j < dim; ++j)
828 for (
unsigned int l = 0; l < dim; ++l)
829 jacobian_grads[point][i][j][l] = result[i][j][l];
840 template <
int dim,
int spacedim,
typename VectorType>
842 maybe_update_jacobian_pushed_forward_grads(
843 const typename ::QProjector<dim>::DataSetDescriptor data_set,
844 const typename ::MappingFEField<dim, spacedim, VectorType>::
848 const std::vector<unsigned int> &fe_to_real,
849 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
854 const unsigned int n_q_points =
855 jacobian_pushed_forward_grads.size();
857 double tmp[spacedim][spacedim][spacedim];
858 for (
unsigned int point = 0; point < n_q_points; ++point)
861 &data.second_derivative(point + data_set, 0);
865 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
867 const unsigned int comp_k =
868 fe.system_to_component_index(k).first;
870 for (
unsigned int j = 0; j < dim; ++j)
871 for (
unsigned int l = 0; l < dim; ++l)
872 result[fe_to_real[comp_k]][j][l] +=
873 (
second[k][j][l] * data.local_dof_values[k]);
877 for (
unsigned int i = 0; i < spacedim; ++i)
878 for (
unsigned int j = 0; j < spacedim; ++j)
879 for (
unsigned int l = 0; l < dim; ++l)
882 result[i][0][l] * data.covariant[point][j][0];
883 for (
unsigned int jr = 1; jr < dim; ++jr)
886 result[i][jr][l] * data.covariant[point][j][jr];
891 for (
unsigned int i = 0; i < spacedim; ++i)
892 for (
unsigned int j = 0; j < spacedim; ++j)
893 for (
unsigned int l = 0; l < spacedim; ++l)
895 jacobian_pushed_forward_grads[point][i][j][l] =
896 tmp[i][j][0] * data.covariant[point][l][0];
897 for (
unsigned int lr = 1; lr < dim; ++lr)
899 jacobian_pushed_forward_grads[point][i][j][l] +=
900 tmp[i][j][lr] * data.covariant[point][l][lr];
913 template <
int dim,
int spacedim,
typename VectorType>
915 maybe_update_jacobian_2nd_derivatives(
916 const typename ::QProjector<dim>::DataSetDescriptor data_set,
917 const typename ::MappingFEField<dim, spacedim, VectorType>::
921 const std::vector<unsigned int> &fe_to_real,
922 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
927 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
929 for (
unsigned int point = 0; point < n_q_points; ++point)
932 &data.third_derivative(point + data_set, 0);
936 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
938 const unsigned int comp_k =
939 fe.system_to_component_index(k).first;
941 for (
unsigned int j = 0; j < dim; ++j)
942 for (
unsigned int l = 0; l < dim; ++l)
943 for (
unsigned int m = 0; m < dim; ++m)
944 result[fe_to_real[comp_k]][j][l][m] +=
945 (third[k][j][l][m] * data.local_dof_values[k]);
950 for (
unsigned int i = 0; i < spacedim; ++i)
951 for (
unsigned int j = 0; j < dim; ++j)
952 for (
unsigned int l = 0; l < dim; ++l)
953 for (
unsigned int m = 0; m < dim; ++m)
954 jacobian_2nd_derivatives[point][i][j][l][m] =
967 template <
int dim,
int spacedim,
typename VectorType>
969 maybe_update_jacobian_pushed_forward_2nd_derivatives(
970 const typename ::QProjector<dim>::DataSetDescriptor data_set,
971 const typename ::MappingFEField<dim, spacedim, VectorType>::
975 const std::vector<unsigned int> &fe_to_real,
976 std::vector<Tensor<4, spacedim>>
977 &jacobian_pushed_forward_2nd_derivatives)
982 const unsigned int n_q_points =
983 jacobian_pushed_forward_2nd_derivatives.size();
985 double tmp[spacedim][spacedim][spacedim][spacedim];
986 for (
unsigned int point = 0; point < n_q_points; ++point)
989 &data.third_derivative(point + data_set, 0);
993 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
995 const unsigned int comp_k =
996 fe.system_to_component_index(k).first;
998 for (
unsigned int j = 0; j < dim; ++j)
999 for (
unsigned int l = 0; l < dim; ++l)
1000 for (
unsigned int m = 0; m < dim; ++m)
1001 result[fe_to_real[comp_k]][j][l][m] +=
1002 (third[k][j][l][m] * data.local_dof_values[k]);
1006 for (
unsigned int i = 0; i < spacedim; ++i)
1007 for (
unsigned int j = 0; j < spacedim; ++j)
1008 for (
unsigned int l = 0; l < dim; ++l)
1009 for (
unsigned int m = 0; m < dim; ++m)
1011 jacobian_pushed_forward_2nd_derivatives
1012 [point][i][j][l][m] =
1013 result[i][0][l][m] * data.covariant[point][j][0];
1014 for (
unsigned int jr = 1; jr < dim; ++jr)
1015 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1017 result[i][jr][l][m] *
1018 data.covariant[point][j][jr];
1022 for (
unsigned int i = 0; i < spacedim; ++i)
1023 for (
unsigned int j = 0; j < spacedim; ++j)
1024 for (
unsigned int l = 0; l < spacedim; ++l)
1025 for (
unsigned int m = 0; m < dim; ++m)
1028 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1030 data.covariant[point][l][0];
1031 for (
unsigned int lr = 1; lr < dim; ++lr)
1033 jacobian_pushed_forward_2nd_derivatives[point][i]
1036 data.covariant[point][l][lr];
1040 for (
unsigned int i = 0; i < spacedim; ++i)
1041 for (
unsigned int j = 0; j < spacedim; ++j)
1042 for (
unsigned int l = 0; l < spacedim; ++l)
1043 for (
unsigned int m = 0; m < spacedim; ++m)
1045 jacobian_pushed_forward_2nd_derivatives
1046 [point][i][j][l][m] =
1047 tmp[i][j][l][0] * data.covariant[point][m][0];
1048 for (
unsigned int mr = 1; mr < dim; ++mr)
1049 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1051 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1063 template <
int dim,
int spacedim,
typename VectorType>
1065 maybe_update_jacobian_3rd_derivatives(
1066 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1067 const typename ::MappingFEField<dim, spacedim, VectorType>::
1071 const std::vector<unsigned int> &fe_to_real,
1072 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1074 const UpdateFlags update_flags = data.update_each;
1077 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1079 for (
unsigned int point = 0; point < n_q_points; ++point)
1082 &data.fourth_derivative(point + data_set, 0);
1086 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1088 const unsigned int comp_k =
1089 fe.system_to_component_index(k).first;
1090 if (fe_mask[comp_k])
1091 for (
unsigned int j = 0; j < dim; ++j)
1092 for (
unsigned int l = 0; l < dim; ++l)
1093 for (
unsigned int m = 0; m < dim; ++m)
1094 for (
unsigned int n = 0; n < dim; ++n)
1095 result[fe_to_real[comp_k]][j][l][m][n] +=
1096 (fourth[k][j][l][m][n] *
1097 data.local_dof_values[k]);
1103 for (
unsigned int i = 0; i < spacedim; ++i)
1104 for (
unsigned int j = 0; j < dim; ++j)
1105 for (
unsigned int l = 0; l < dim; ++l)
1106 for (
unsigned int m = 0; m < dim; ++m)
1107 for (
unsigned int n = 0; n < dim; ++n)
1108 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1109 result[i][j][l][m][n];
1121 template <
int dim,
int spacedim,
typename VectorType>
1123 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1124 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1125 const typename ::MappingFEField<dim, spacedim, VectorType>::
1129 const std::vector<unsigned int> &fe_to_real,
1130 std::vector<Tensor<5, spacedim>>
1131 &jacobian_pushed_forward_3rd_derivatives)
1133 const UpdateFlags update_flags = data.update_each;
1136 const unsigned int n_q_points =
1137 jacobian_pushed_forward_3rd_derivatives.size();
1139 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1140 for (
unsigned int point = 0; point < n_q_points; ++point)
1143 &data.fourth_derivative(point + data_set, 0);
1147 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1149 const unsigned int comp_k =
1150 fe.system_to_component_index(k).first;
1151 if (fe_mask[comp_k])
1152 for (
unsigned int j = 0; j < dim; ++j)
1153 for (
unsigned int l = 0; l < dim; ++l)
1154 for (
unsigned int m = 0; m < dim; ++m)
1155 for (
unsigned int n = 0; n < dim; ++n)
1156 result[fe_to_real[comp_k]][j][l][m][n] +=
1157 (fourth[k][j][l][m][n] *
1158 data.local_dof_values[k]);
1162 for (
unsigned int i = 0; i < spacedim; ++i)
1163 for (
unsigned int j = 0; j < spacedim; ++j)
1164 for (
unsigned int l = 0; l < dim; ++l)
1165 for (
unsigned int m = 0; m < dim; ++m)
1166 for (
unsigned int n = 0; n < dim; ++n)
1168 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1169 data.covariant[point][j][0];
1170 for (
unsigned int jr = 1; jr < dim; ++jr)
1171 tmp[i][j][l][m][n] +=
1172 result[i][jr][l][m][n] *
1173 data.covariant[point][j][jr];
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 < spacedim; ++l)
1180 for (
unsigned int m = 0; m < dim; ++m)
1181 for (
unsigned int n = 0; n < dim; ++n)
1183 jacobian_pushed_forward_3rd_derivatives
1184 [point][i][j][l][m][n] =
1185 tmp[i][j][0][m][n] *
1186 data.covariant[point][l][0];
1187 for (
unsigned int lr = 1; lr < dim; ++lr)
1188 jacobian_pushed_forward_3rd_derivatives[point][i]
1191 tmp[i][j][lr][m][n] *
1192 data.covariant[point][l][lr];
1196 for (
unsigned int i = 0; i < spacedim; ++i)
1197 for (
unsigned int j = 0; j < spacedim; ++j)
1198 for (
unsigned int l = 0; l < spacedim; ++l)
1199 for (
unsigned int m = 0; m < spacedim; ++m)
1200 for (
unsigned int n = 0; n < dim; ++n)
1202 tmp[i][j][l][m][n] =
1203 jacobian_pushed_forward_3rd_derivatives[point][i]
1206 data.covariant[point][m][0];
1207 for (
unsigned int mr = 1; mr < dim; ++mr)
1208 tmp[i][j][l][m][n] +=
1209 jacobian_pushed_forward_3rd_derivatives[point]
1212 data.covariant[point][m][mr];
1216 for (
unsigned int i = 0; i < spacedim; ++i)
1217 for (
unsigned int j = 0; j < spacedim; ++j)
1218 for (
unsigned int l = 0; l < spacedim; ++l)
1219 for (
unsigned int m = 0; m < spacedim; ++m)
1220 for (
unsigned int n = 0; n < spacedim; ++n)
1222 jacobian_pushed_forward_3rd_derivatives
1223 [point][i][j][l][m][n] =
1224 tmp[i][j][l][m][0] *
1225 data.covariant[point][n][0];
1226 for (
unsigned int nr = 1; nr < dim; ++nr)
1227 jacobian_pushed_forward_3rd_derivatives[point][i]
1230 tmp[i][j][l][m][nr] *
1231 data.covariant[point][n][nr];
1247 template <
int dim,
int spacedim,
typename VectorType>
1249 maybe_compute_face_data(
1250 const ::Mapping<dim, spacedim> &mapping,
1251 const typename ::Triangulation<dim, spacedim>::cell_iterator
1253 const unsigned int face_no,
1254 const unsigned int subface_no,
1256 const typename ::MappingFEField<dim, spacedim, VectorType>::
1261 const UpdateFlags update_flags = data.update_each;
1265 const unsigned int n_q_points = output_data.boundary_forms.size();
1274 for (
unsigned int d = 0; d != dim - 1; ++d)
1276 Assert(face_no + cell->n_faces() * d <
1277 data.unit_tangentials.size(),
1280 data.aux[d].size() <=
1281 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1286 data.unit_tangentials[face_no + cell->n_faces() * d]),
1294 if (dim == spacedim)
1296 for (
unsigned int i = 0; i < n_q_points; ++i)
1304 output_data.boundary_forms[i][0] =
1305 (face_no == 0 ? -1 : +1);
1308 output_data.boundary_forms[i] =
1309 cross_product_2d(data.aux[0][i]);
1312 output_data.boundary_forms[i] =
1313 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1329 for (
unsigned int point = 0; point < n_q_points; ++point)
1334 output_data.boundary_forms[point] =
1335 data.contravariant[point].transpose()[0];
1336 output_data.boundary_forms[point] /=
1337 (face_no == 0 ? -1. : +1.) *
1338 output_data.boundary_forms[point].norm();
1347 cross_product_3d(DX_t[0], DX_t[1]);
1348 cell_normal /= cell_normal.
norm();
1352 output_data.boundary_forms[point] =
1353 cross_product_3d(data.aux[0][point], cell_normal);
1359 for (
unsigned int i = 0; i < output_data.boundary_forms.size();
1364 output_data.JxW_values[i] =
1365 output_data.boundary_forms[i].norm() *
1366 data.quadrature_weights[i + data_set];
1371 const double area_ratio =
1373 cell->subface_case(face_no), subface_no);
1374 output_data.JxW_values[i] *= area_ratio;
1379 output_data.normal_vectors[i] =
1381 output_data.boundary_forms[i].norm());
1392 template <
int dim,
int spacedim,
typename VectorType>
1394 do_fill_fe_face_values(
1395 const ::Mapping<dim, spacedim> &mapping,
1396 const typename ::Triangulation<dim, spacedim>::cell_iterator
1398 const unsigned int face_no,
1399 const unsigned int subface_no,
1400 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1401 const typename ::MappingFEField<dim, spacedim, VectorType>::
1405 const std::vector<unsigned int> &fe_to_real,
1409 maybe_compute_q_points<dim, spacedim, VectorType>(
1415 output_data.quadrature_points);
1417 maybe_update_Jacobians<dim, spacedim, VectorType>(
1418 data_set, data, fe, fe_mask, fe_to_real);
1420 const UpdateFlags update_flags = data.update_each;
1421 const unsigned int n_q_points = data.contravariant.size();
1424 for (
unsigned int point = 0; point < n_q_points; ++point)
1425 output_data.jacobians[point] = data.contravariant[point];
1428 for (
unsigned int point = 0; point < n_q_points; ++point)
1429 output_data.inverse_jacobians[point] =
1430 data.covariant[point].transpose();
1432 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1433 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1435 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1441 output_data.jacobian_pushed_forward_grads);
1443 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1449 output_data.jacobian_2nd_derivatives);
1451 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1459 output_data.jacobian_pushed_forward_2nd_derivatives);
1461 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1467 output_data.jacobian_3rd_derivatives);
1469 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1477 output_data.jacobian_pushed_forward_3rd_derivatives);
1479 maybe_compute_face_data<dim, spacedim, VectorType>(
1480 mapping, cell, face_no, subface_no, data_set, data, output_data);