35#include <boost/container/small_vector.hpp>
47template <
int dim,
int spacedim>
57template <
int dim,
int spacedim>
78template <
int dim,
int spacedim>
87 const unsigned int n_q_points = q.
size();
134template <
int dim,
int spacedim>
137 const UpdateFlags update_flags,
139 const unsigned int n_original_q_points)
151 const auto reference_cell = this->
fe.reference_cell();
152 const auto n_faces = reference_cell.n_faces();
154 for (
unsigned int i = 0; i < n_faces; ++i)
174template <
int dim,
int spacedim>
183 const auto &tensor_pols = fe_poly->get_poly_space();
186 const unsigned int n_points = unit_points.size();
188 std::vector<double> values;
189 std::vector<Tensor<1, dim>> grads;
203 std::vector<Tensor<2, dim>> grad2;
211 std::vector<Tensor<3, dim>> grad3;
219 std::vector<Tensor<4, dim>> grad4;
232 for (
unsigned int point = 0; point < n_points; ++point)
234 tensor_pols.evaluate(
235 unit_points[point], values, grads, grad2, grad3, grad4);
239 shape(point, i) = values[i];
272 template <
int dim,
int spacedim>
274 maybe_compute_q_points(
276 const typename ::MappingFE<dim, spacedim>::InternalData &data,
278 const unsigned int n_q_points)
280 const UpdateFlags update_flags = data.update_each;
283 for (
unsigned int point = 0; point < n_q_points; ++point)
285 const double *shape = &data.shape(point + data_set, 0);
287 (shape[0] * data.mapping_support_points[0]);
288 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
289 for (
unsigned int i = 0; i < spacedim; ++i)
290 result[i] += shape[k] * data.mapping_support_points[k][i];
291 quadrature_points[point] = result;
305 template <
int dim,
int spacedim>
307 maybe_update_Jacobians(
309 const typename ::QProjector<dim>::DataSetDescriptor data_set,
310 const typename ::MappingFE<dim, spacedim>::InternalData &data,
311 const unsigned int n_q_points)
313 const UpdateFlags update_flags = data.update_each;
321 std::fill(data.contravariant.begin(),
322 data.contravariant.end(),
327 for (
unsigned int point = 0; point < n_q_points; ++point)
329 double result[spacedim][dim];
333 for (
unsigned int i = 0; i < spacedim; ++i)
334 for (
unsigned int j = 0; j < dim; ++j)
335 result[i][j] = data.derivative(point + data_set, 0)[j] *
336 data.mapping_support_points[0][i];
337 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
338 for (
unsigned int i = 0; i < spacedim; ++i)
339 for (
unsigned int j = 0; j < dim; ++j)
341 data.derivative(point + data_set, k)[j] *
342 data.mapping_support_points[k][i];
349 for (
unsigned int i = 0; i < spacedim; ++i)
350 for (
unsigned int j = 0; j < dim; ++j)
351 data.contravariant[point][i][j] = result[i][j];
358 for (
unsigned int point = 0; point < n_q_points; ++point)
360 data.covariant[point] =
361 (data.contravariant[point]).covariant_form();
368 for (
unsigned int point = 0; point < n_q_points; ++point)
369 data.volume_elements[point] =
370 data.contravariant[point].determinant();
380 template <
int dim,
int spacedim>
382 maybe_update_jacobian_grads(
385 const typename ::MappingFE<dim, spacedim>::InternalData &data,
387 const unsigned int n_q_points)
389 const UpdateFlags update_flags = data.update_each;
395 for (
unsigned int point = 0; point < n_q_points; ++point)
398 &data.second_derivative(point + data_set, 0);
399 double result[spacedim][dim][dim];
400 for (
unsigned int i = 0; i < spacedim; ++i)
401 for (
unsigned int j = 0; j < dim; ++j)
402 for (
unsigned int l = 0; l < dim; ++l)
404 (second[0][j][l] * data.mapping_support_points[0][i]);
405 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
406 for (
unsigned int i = 0; i < spacedim; ++i)
407 for (
unsigned int j = 0; j < dim; ++j)
408 for (
unsigned int l = 0; l < dim; ++l)
411 data.mapping_support_points[k][i]);
413 for (
unsigned int i = 0; i < spacedim; ++i)
414 for (
unsigned int j = 0; j < dim; ++j)
415 for (
unsigned int l = 0; l < dim; ++l)
416 jacobian_grads[point][i][j][l] = result[i][j][l];
427 template <
int dim,
int spacedim>
429 maybe_update_jacobian_pushed_forward_grads(
432 const typename ::MappingFE<dim, spacedim>::InternalData &data,
434 const unsigned int n_q_points)
436 const UpdateFlags update_flags = data.update_each;
440 jacobian_pushed_forward_grads.size() + 1);
444 double tmp[spacedim][spacedim][spacedim];
445 for (
unsigned int point = 0; point < n_q_points; ++point)
448 &data.second_derivative(point + data_set, 0);
449 double result[spacedim][dim][dim];
450 for (
unsigned int i = 0; i < spacedim; ++i)
451 for (
unsigned int j = 0; j < dim; ++j)
452 for (
unsigned int l = 0; l < dim; ++l)
453 result[i][j][l] = (second[0][j][l] *
454 data.mapping_support_points[0][i]);
455 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
456 for (
unsigned int i = 0; i < spacedim; ++i)
457 for (
unsigned int j = 0; j < dim; ++j)
458 for (
unsigned int l = 0; l < dim; ++l)
461 data.mapping_support_points[k][i]);
464 for (
unsigned int i = 0; i < spacedim; ++i)
465 for (
unsigned int j = 0; j < spacedim; ++j)
466 for (
unsigned int l = 0; l < dim; ++l)
469 result[i][0][l] * data.covariant[point][j][0];
470 for (
unsigned int jr = 1; jr < dim; ++jr)
472 tmp[i][j][l] += result[i][jr][l] *
473 data.covariant[point][j][jr];
478 for (
unsigned int i = 0; i < spacedim; ++i)
479 for (
unsigned int j = 0; j < spacedim; ++j)
480 for (
unsigned int l = 0; l < spacedim; ++l)
482 jacobian_pushed_forward_grads[point][i][j][l] =
483 tmp[i][j][0] * data.covariant[point][l][0];
484 for (
unsigned int lr = 1; lr < dim; ++lr)
486 jacobian_pushed_forward_grads[point][i][j][l] +=
487 tmp[i][j][lr] * data.covariant[point][l][lr];
501 template <
int dim,
int spacedim>
503 maybe_update_jacobian_2nd_derivatives(
506 const typename ::MappingFE<dim, spacedim>::InternalData &data,
508 const unsigned int n_q_points)
510 const UpdateFlags update_flags = data.update_each;
517 for (
unsigned int point = 0; point < n_q_points; ++point)
520 &data.third_derivative(point + data_set, 0);
521 double result[spacedim][dim][dim][dim];
522 for (
unsigned int i = 0; i < spacedim; ++i)
523 for (
unsigned int j = 0; j < dim; ++j)
524 for (
unsigned int l = 0; l < dim; ++l)
525 for (
unsigned int m = 0; m < dim; ++m)
528 data.mapping_support_points[0][i]);
529 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
530 for (
unsigned int i = 0; i < spacedim; ++i)
531 for (
unsigned int j = 0; j < dim; ++j)
532 for (
unsigned int l = 0; l < dim; ++l)
533 for (
unsigned int m = 0; m < dim; ++m)
534 result[i][j][l][m] +=
536 data.mapping_support_points[k][i]);
538 for (
unsigned int i = 0; i < spacedim; ++i)
539 for (
unsigned int j = 0; j < dim; ++j)
540 for (
unsigned int l = 0; l < dim; ++l)
541 for (
unsigned int m = 0; m < dim; ++m)
542 jacobian_2nd_derivatives[point][i][j][l][m] =
556 template <
int dim,
int spacedim>
558 maybe_update_jacobian_pushed_forward_2nd_derivatives(
561 const typename ::MappingFE<dim, spacedim>::InternalData &data,
563 &jacobian_pushed_forward_2nd_derivatives,
564 const unsigned int n_q_points)
566 const UpdateFlags update_flags = data.update_each;
570 jacobian_pushed_forward_2nd_derivatives.size() +
575 double tmp[spacedim][spacedim][spacedim][spacedim];
576 for (
unsigned int point = 0; point < n_q_points; ++point)
579 &data.third_derivative(point + data_set, 0);
580 double result[spacedim][dim][dim][dim];
581 for (
unsigned int i = 0; i < spacedim; ++i)
582 for (
unsigned int j = 0; j < dim; ++j)
583 for (
unsigned int l = 0; l < dim; ++l)
584 for (
unsigned int m = 0; m < dim; ++m)
587 data.mapping_support_points[0][i]);
588 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
589 for (
unsigned int i = 0; i < spacedim; ++i)
590 for (
unsigned int j = 0; j < dim; ++j)
591 for (
unsigned int l = 0; l < dim; ++l)
592 for (
unsigned int m = 0; m < dim; ++m)
593 result[i][j][l][m] +=
595 data.mapping_support_points[k][i]);
598 for (
unsigned int i = 0; i < spacedim; ++i)
599 for (
unsigned int j = 0; j < spacedim; ++j)
600 for (
unsigned int l = 0; l < dim; ++l)
601 for (
unsigned int m = 0; m < dim; ++m)
603 jacobian_pushed_forward_2nd_derivatives
604 [point][i][j][l][m] =
606 data.covariant[point][j][0];
607 for (
unsigned int jr = 1; jr < dim; ++jr)
608 jacobian_pushed_forward_2nd_derivatives[point]
611 result[i][jr][l][m] *
612 data.covariant[point][j][jr];
616 for (
unsigned int i = 0; i < spacedim; ++i)
617 for (
unsigned int j = 0; j < spacedim; ++j)
618 for (
unsigned int l = 0; l < spacedim; ++l)
619 for (
unsigned int m = 0; m < dim; ++m)
622 jacobian_pushed_forward_2nd_derivatives[point]
625 data.covariant[point][l][0];
626 for (
unsigned int lr = 1; lr < dim; ++lr)
628 jacobian_pushed_forward_2nd_derivatives
629 [point][i][j][lr][m] *
630 data.covariant[point][l][lr];
634 for (
unsigned int i = 0; i < spacedim; ++i)
635 for (
unsigned int j = 0; j < spacedim; ++j)
636 for (
unsigned int l = 0; l < spacedim; ++l)
637 for (
unsigned int m = 0; m < spacedim; ++m)
639 jacobian_pushed_forward_2nd_derivatives
640 [point][i][j][l][m] =
641 tmp[i][j][l][0] * data.covariant[point][m][0];
642 for (
unsigned int mr = 1; mr < dim; ++mr)
643 jacobian_pushed_forward_2nd_derivatives[point]
647 data.covariant[point][m][mr];
660 template <
int dim,
int spacedim>
662 maybe_update_jacobian_3rd_derivatives(
665 const typename ::MappingFE<dim, spacedim>::InternalData &data,
667 const unsigned int n_q_points)
669 const UpdateFlags update_flags = data.update_each;
676 for (
unsigned int point = 0; point < n_q_points; ++point)
679 &data.fourth_derivative(point + data_set, 0);
680 double result[spacedim][dim][dim][dim][dim];
681 for (
unsigned int i = 0; i < spacedim; ++i)
682 for (
unsigned int j = 0; j < dim; ++j)
683 for (
unsigned int l = 0; l < dim; ++l)
684 for (
unsigned int m = 0; m < dim; ++m)
685 for (
unsigned int n = 0; n < dim; ++n)
686 result[i][j][l][m][n] =
687 (fourth[0][j][l][m][n] *
688 data.mapping_support_points[0][i]);
689 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
690 for (
unsigned int i = 0; i < spacedim; ++i)
691 for (
unsigned int j = 0; j < dim; ++j)
692 for (
unsigned int l = 0; l < dim; ++l)
693 for (
unsigned int m = 0; m < dim; ++m)
694 for (
unsigned int n = 0; n < dim; ++n)
695 result[i][j][l][m][n] +=
696 (fourth[k][j][l][m][n] *
697 data.mapping_support_points[k][i]);
699 for (
unsigned int i = 0; i < spacedim; ++i)
700 for (
unsigned int j = 0; j < dim; ++j)
701 for (
unsigned int l = 0; l < dim; ++l)
702 for (
unsigned int m = 0; m < dim; ++m)
703 for (
unsigned int n = 0; n < dim; ++n)
704 jacobian_3rd_derivatives[point][i][j][l][m][n] =
705 result[i][j][l][m][n];
718 template <
int dim,
int spacedim>
720 maybe_update_jacobian_pushed_forward_3rd_derivatives(
723 const typename ::MappingFE<dim, spacedim>::InternalData &data,
725 &jacobian_pushed_forward_3rd_derivatives,
726 const unsigned int n_q_points)
728 const UpdateFlags update_flags = data.update_each;
732 jacobian_pushed_forward_3rd_derivatives.size() +
737 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
738 for (
unsigned int point = 0; point < n_q_points; ++point)
741 &data.fourth_derivative(point + data_set, 0);
742 double result[spacedim][dim][dim][dim][dim];
743 for (
unsigned int i = 0; i < spacedim; ++i)
744 for (
unsigned int j = 0; j < dim; ++j)
745 for (
unsigned int l = 0; l < dim; ++l)
746 for (
unsigned int m = 0; m < dim; ++m)
747 for (
unsigned int n = 0; n < dim; ++n)
748 result[i][j][l][m][n] =
749 (fourth[0][j][l][m][n] *
750 data.mapping_support_points[0][i]);
751 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
752 for (
unsigned int i = 0; i < spacedim; ++i)
753 for (
unsigned int j = 0; j < dim; ++j)
754 for (
unsigned int l = 0; l < dim; ++l)
755 for (
unsigned int m = 0; m < dim; ++m)
756 for (
unsigned int n = 0; n < dim; ++n)
757 result[i][j][l][m][n] +=
758 (fourth[k][j][l][m][n] *
759 data.mapping_support_points[k][i]);
762 for (
unsigned int i = 0; i < spacedim; ++i)
763 for (
unsigned int j = 0; j < spacedim; ++j)
764 for (
unsigned int l = 0; l < dim; ++l)
765 for (
unsigned int m = 0; m < dim; ++m)
766 for (
unsigned int n = 0; n < dim; ++n)
769 result[i][0][l][m][n] *
770 data.covariant[point][j][0];
771 for (
unsigned int jr = 1; jr < dim; ++jr)
772 tmp[i][j][l][m][n] +=
773 result[i][jr][l][m][n] *
774 data.covariant[point][j][jr];
778 for (
unsigned int i = 0; i < spacedim; ++i)
779 for (
unsigned int j = 0; j < spacedim; ++j)
780 for (
unsigned int l = 0; l < spacedim; ++l)
781 for (
unsigned int m = 0; m < dim; ++m)
782 for (
unsigned int n = 0; n < dim; ++n)
784 jacobian_pushed_forward_3rd_derivatives
785 [point][i][j][l][m][n] =
787 data.covariant[point][l][0];
788 for (
unsigned int lr = 1; lr < dim; ++lr)
789 jacobian_pushed_forward_3rd_derivatives
790 [point][i][j][l][m][n] +=
791 tmp[i][j][lr][m][n] *
792 data.covariant[point][l][lr];
796 for (
unsigned int i = 0; i < spacedim; ++i)
797 for (
unsigned int j = 0; j < spacedim; ++j)
798 for (
unsigned int l = 0; l < spacedim; ++l)
799 for (
unsigned int m = 0; m < spacedim; ++m)
800 for (
unsigned int n = 0; n < dim; ++n)
803 jacobian_pushed_forward_3rd_derivatives
804 [point][i][j][l][0][n] *
805 data.covariant[point][m][0];
806 for (
unsigned int mr = 1; mr < dim; ++mr)
807 tmp[i][j][l][m][n] +=
808 jacobian_pushed_forward_3rd_derivatives
809 [point][i][j][l][mr][n] *
810 data.covariant[point][m][mr];
814 for (
unsigned int i = 0; i < spacedim; ++i)
815 for (
unsigned int j = 0; j < spacedim; ++j)
816 for (
unsigned int l = 0; l < spacedim; ++l)
817 for (
unsigned int m = 0; m < spacedim; ++m)
818 for (
unsigned int n = 0; n < spacedim; ++n)
820 jacobian_pushed_forward_3rd_derivatives
821 [point][i][j][l][m][n] =
823 data.covariant[point][n][0];
824 for (
unsigned int nr = 1; nr < dim; ++nr)
825 jacobian_pushed_forward_3rd_derivatives
826 [point][i][j][l][m][n] +=
827 tmp[i][j][l][m][nr] *
828 data.covariant[point][n][nr];
840template <
int dim,
int spacedim>
846 ExcMessage(
"It only makes sense to create polynomial mappings "
847 "with a polynomial degree greater or equal to one."));
852 const auto &mapping_support_points =
fe.get_unit_support_points();
854 const auto reference_cell =
fe.reference_cell();
856 const unsigned int n_points = mapping_support_points.size();
857 const unsigned int n_shape_functions = reference_cell.n_vertices();
862 for (
unsigned int point = 0; point < n_points; ++point)
863 for (
unsigned int i = 0; i < n_shape_functions; ++i)
865 reference_cell.d_linear_shape_function(mapping_support_points[point],
871template <
int dim,
int spacedim>
880template <
int dim,
int spacedim>
881std::unique_ptr<Mapping<dim, spacedim>>
884 return std::make_unique<MappingFE<dim, spacedim>>(*this);
889template <
int dim,
int spacedim>
898template <
int dim,
int spacedim>
904 const std::vector<Point<spacedim>> support_points =
909 for (
unsigned int i = 0; i < this->
fe->n_dofs_per_cell(); ++i)
910 mapped_point += support_points[i] * this->
fe->shape_value(i, p);
917template <
int dim,
int spacedim>
923 const std::vector<Point<spacedim>> support_points =
926 const double eps = 1.e-12 * cell->
diameter();
927 const unsigned int loop_limit = 10;
931 unsigned int loop = 0;
947 for (
unsigned int i = 0; i < this->
fe->n_dofs_per_cell(); ++i)
949 mapped_point += support_points[i] * this->
fe->shape_value(i, p_unit);
950 const auto grad_F_i = this->
fe->shape_grad(i, p_unit);
951 const auto hessian_F_i = this->
fe->shape_grad_grad(i, p_unit);
952 for (
unsigned int j = 0; j < dim; ++j)
954 grad_FT[j] += grad_F_i[j] * support_points[i];
955 for (
unsigned int l = 0; l < dim; ++l)
956 hess_FT[j][l] += hessian_F_i[j][l] * support_points[i];
961 const auto residual = p - mapped_point;
968 if (grad_FT_residual.norm() <= eps)
973 for (
unsigned int j = 0; j < dim; ++j)
974 for (
unsigned int l = 0; l < dim; ++l)
975 corrected_metric_tensor[j][l] =
976 -grad_FT[j] * grad_FT[l] + hess_FT[j][l] * residual;
979 const auto g_inverse =
invert(corrected_metric_tensor);
980 p_unit -=
Point<dim>(g_inverse * grad_FT_residual);
984 while (loop < loop_limit);
998template <
int dim,
int spacedim>
1008 UpdateFlags out = in;
1009 for (
unsigned int i = 0; i < 5; ++i)
1054template <
int dim,
int spacedim>
1055std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1059 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1060 std::make_unique<InternalData>(*this->
fe);
1067template <
int dim,
int spacedim>
1068std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1070 const UpdateFlags update_flags,
1073 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1074 std::make_unique<InternalData>(*this->
fe);
1078 this->
fe->reference_cell(), quadrature),
1086template <
int dim,
int spacedim>
1087std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1089 const UpdateFlags update_flags,
1092 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1093 std::make_unique<InternalData>(*this->
fe);
1097 this->
fe->reference_cell(), quadrature),
1105template <
int dim,
int spacedim>
1120 const unsigned int n_q_points = quadrature.
size();
1142 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1148 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1149 computed_cell_similarity,
1154 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1155 computed_cell_similarity,
1161 internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
1163 spacedim>(computed_cell_similarity,
1169 internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
1171 spacedim>(computed_cell_similarity,
1177 internal::MappingFEImplementation::
1178 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1179 computed_cell_similarity,
1185 internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
1187 spacedim>(computed_cell_similarity,
1193 internal::MappingFEImplementation::
1194 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1195 computed_cell_similarity,
1201 const UpdateFlags update_flags = data.
update_each;
1202 const std::vector<double> &weights = quadrature.
get_weights();
1218 for (
unsigned int point = 0; point < n_q_points; ++point)
1220 if (dim == spacedim)
1233 cell->
center(), det, point)));
1235 output_data.
JxW_values[point] = weights[point] * det;
1243 for (
unsigned int i = 0; i < spacedim; ++i)
1244 for (
unsigned int j = 0; j < dim; ++j)
1248 for (
unsigned int i = 0; i < dim; ++i)
1249 for (
unsigned int j = 0; j < dim; ++j)
1250 G[i][j] = DX_t[i] * DX_t[j];
1255 if (computed_cell_similarity ==
1266 Assert(spacedim == dim + 1,
1268 "There is no (unique) cell normal for " +
1270 "-dimensional cells in " +
1272 "-dimensional space. This only works if the "
1273 "space dimension is one greater than the "
1274 "dimensionality of the mesh cells."));
1301 for (
unsigned int point = 0; point < n_q_points; ++point)
1310 for (
unsigned int point = 0; point < n_q_points; ++point)
1315 return computed_cell_similarity;
1322 namespace MappingFEImplementation
1336 template <
int dim,
int spacedim>
1338 maybe_compute_face_data(
1339 const ::MappingFE<dim, spacedim> &
mapping,
1340 const typename ::Triangulation<dim, spacedim>::cell_iterator
1342 const unsigned int face_no,
1343 const unsigned int subface_no,
1344 const unsigned int n_q_points,
1346 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1350 const UpdateFlags update_flags = data.update_each;
1373 for (
unsigned int d = 0; d != dim - 1; ++d)
1375 Assert(face_no + cell->n_faces() * d <
1376 data.unit_tangentials.size(),
1379 data.aux[d].size() <=
1380 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1385 data.unit_tangentials[face_no + cell->n_faces() * d]),
1396 if (dim == spacedim)
1398 for (
unsigned int i = 0; i < n_q_points; ++i)
1408 (face_no == 0 ? -1 : +1);
1432 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1438 data.contravariant[
point].transpose()[0];
1440 (face_no == 0 ? -1. : +1.) *
1446 const DerivativeForm<1, spacedim, dim> DX_t =
1447 data.contravariant[
point].transpose();
1449 Tensor<1, spacedim> cell_normal =
1451 cell_normal /= cell_normal.
norm();
1463 for (
unsigned int i = 0; i < n_q_points; ++i)
1467 data.quadrature_weights[i + data_set];
1472 const double area_ratio =
1474 cell->subface_case(face_no), subface_no);
1483 for (
unsigned int i = 0; i < n_q_points; ++i)
1489 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1490 output_data.
jacobians[point] = data.contravariant[point];
1493 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1495 data.covariant[point].transpose();
1506 template <
int dim,
int spacedim>
1509 const ::MappingFE<dim, spacedim> &
mapping,
1510 const typename ::Triangulation<dim, spacedim>::cell_iterator
1512 const unsigned int face_no,
1513 const unsigned int subface_no,
1514 const typename QProjector<dim>::DataSetDescriptor data_set,
1515 const Quadrature<dim - 1> &quadrature,
1516 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1517 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1520 const unsigned int n_q_points = quadrature.
size();
1522 maybe_compute_q_points<dim, spacedim>(data_set,
1535 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1541 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1547 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1553 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1559 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1581template <
int dim,
int spacedim>
1585 const unsigned int face_no,
1609 internal::MappingFEImplementation::do_fill_fe_face_values(
1619 quadrature[quadrature.
size() == 1 ? 0 : face_no],
1626template <
int dim,
int spacedim>
1630 const unsigned int face_no,
1631 const unsigned int subface_no,
1655 internal::MappingFEImplementation::do_fill_fe_face_values(
1676 namespace MappingFEImplementation
1680 template <
int dim,
int spacedim,
int rank>
1697 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1698 &mapping_data) !=
nullptr),
1700 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1702 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1705 switch (mapping_kind)
1712 "update_contravariant_transformation"));
1714 for (
unsigned int i = 0; i < input.size(); ++i)
1726 "update_contravariant_transformation"));
1730 "update_volume_elements"));
1735 for (
unsigned int i = 0; i < input.size(); ++i)
1739 output[i] /= data.volume_elements[i];
1751 "update_covariant_transformation"));
1753 for (
unsigned int i = 0; i < input.size(); ++i)
1765 template <
int dim,
int spacedim,
int rank>
1768 const ArrayView<
const Tensor<rank, dim>> &input,
1770 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1771 const ArrayView<Tensor<rank, spacedim>> &output)
1776 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1777 &mapping_data) !=
nullptr),
1779 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1781 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1784 switch (mapping_kind)
1791 "update_covariant_transformation"));
1795 "update_contravariant_transformation"));
1798 for (
unsigned int i = 0; i < output.size(); ++i)
1800 const DerivativeForm<1, spacedim, dim> A =
1815 "update_covariant_transformation"));
1818 for (
unsigned int i = 0; i < output.size(); ++i)
1820 const DerivativeForm<1, spacedim, dim> A =
1835 "update_covariant_transformation"));
1839 "update_contravariant_transformation"));
1843 "update_volume_elements"));
1846 for (
unsigned int i = 0; i < output.size(); ++i)
1848 const DerivativeForm<1, spacedim, dim> A =
1850 const Tensor<2, spacedim> T =
1855 output[i] /= data.volume_elements[i];
1868 template <
int dim,
int spacedim>
1871 const ArrayView<
const Tensor<3, dim>> &input,
1873 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1874 const ArrayView<Tensor<3, spacedim>> &output)
1879 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1880 &mapping_data) !=
nullptr),
1882 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1884 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1887 switch (mapping_kind)
1894 "update_covariant_transformation"));
1898 "update_contravariant_transformation"));
1900 for (
unsigned int q = 0; q < output.size(); ++q)
1901 for (
unsigned int i = 0; i < spacedim; ++i)
1903 double tmp1[dim][dim];
1904 for (
unsigned int J = 0; J < dim; ++J)
1905 for (
unsigned int K = 0;
K < dim; ++
K)
1908 data.contravariant[q][i][0] * input[q][0][J][
K];
1909 for (
unsigned int I = 1; I < dim; ++I)
1911 data.contravariant[q][i][I] * input[q][I][J][K];
1913 for (
unsigned int j = 0; j < spacedim; ++j)
1916 for (
unsigned int K = 0;
K < dim; ++
K)
1918 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1919 for (
unsigned int J = 1; J < dim; ++J)
1920 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1922 for (
unsigned int k = 0; k < spacedim; ++k)
1924 output[q][i][j][k] =
1925 data.covariant[q][k][0] * tmp2[0];
1926 for (
unsigned int K = 1;
K < dim; ++
K)
1927 output[q][i][j][k] +=
1928 data.covariant[q][k][K] * tmp2[K];
1940 "update_covariant_transformation"));
1942 for (
unsigned int q = 0; q < output.size(); ++q)
1943 for (
unsigned int i = 0; i < spacedim; ++i)
1945 double tmp1[dim][dim];
1946 for (
unsigned int J = 0; J < dim; ++J)
1947 for (
unsigned int K = 0;
K < dim; ++
K)
1950 data.covariant[q][i][0] * input[q][0][J][
K];
1951 for (
unsigned int I = 1; I < dim; ++I)
1953 data.covariant[q][i][I] * input[q][I][J][K];
1955 for (
unsigned int j = 0; j < spacedim; ++j)
1958 for (
unsigned int K = 0;
K < dim; ++
K)
1960 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1961 for (
unsigned int J = 1; J < dim; ++J)
1962 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1964 for (
unsigned int k = 0; k < spacedim; ++k)
1966 output[q][i][j][k] =
1967 data.covariant[q][k][0] * tmp2[0];
1968 for (
unsigned int K = 1;
K < dim; ++
K)
1969 output[q][i][j][k] +=
1970 data.covariant[q][k][K] * tmp2[K];
1983 "update_covariant_transformation"));
1987 "update_contravariant_transformation"));
1991 "update_volume_elements"));
1993 for (
unsigned int q = 0; q < output.size(); ++q)
1994 for (
unsigned int i = 0; i < spacedim; ++i)
1997 for (
unsigned int I = 0; I < dim; ++I)
1999 data.contravariant[q][i][I] / data.volume_elements[q];
2000 double tmp1[dim][dim];
2001 for (
unsigned int J = 0; J < dim; ++J)
2002 for (
unsigned int K = 0;
K < dim; ++
K)
2004 tmp1[J][
K] = factor[0] * input[q][0][J][
K];
2005 for (
unsigned int I = 1; I < dim; ++I)
2006 tmp1[J][K] += factor[I] * input[q][I][J][K];
2008 for (
unsigned int j = 0; j < spacedim; ++j)
2011 for (
unsigned int K = 0;
K < dim; ++
K)
2013 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
2014 for (
unsigned int J = 1; J < dim; ++J)
2015 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2017 for (
unsigned int k = 0; k < spacedim; ++k)
2019 output[q][i][j][k] =
2020 data.covariant[q][k][0] * tmp2[0];
2021 for (
unsigned int K = 1;
K < dim; ++
K)
2022 output[q][i][j][k] +=
2023 data.covariant[q][k][K] * tmp2[K];
2038 template <
int dim,
int spacedim,
int rank>
2041 const ArrayView<
const DerivativeForm<rank, dim, spacedim>> &input,
2043 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2044 const ArrayView<Tensor<rank + 1, spacedim>> &output)
2049 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
2050 &mapping_data) !=
nullptr),
2052 const typename ::MappingFE<dim, spacedim>::InternalData &data =
2054 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
2057 switch (mapping_kind)
2064 "update_covariant_transformation"));
2066 for (
unsigned int i = 0; i < output.size(); ++i)
2081template <
int dim,
int spacedim>
2085 const MappingKind mapping_kind,
2089 internal::MappingFEImplementation::transform_fields(input,
2097template <
int dim,
int spacedim>
2101 const MappingKind mapping_kind,
2105 internal::MappingFEImplementation::transform_differential_forms(input,
2113template <
int dim,
int spacedim>
2117 const MappingKind mapping_kind,
2121 switch (mapping_kind)
2124 internal::MappingFEImplementation::transform_fields(input,
2133 internal::MappingFEImplementation::transform_gradients(input,
2145template <
int dim,
int spacedim>
2149 const MappingKind mapping_kind,
2158 switch (mapping_kind)
2164 "update_covariant_transformation"));
2166 for (
unsigned int q = 0; q < output.size(); ++q)
2167 for (
unsigned int i = 0; i < spacedim; ++i)
2168 for (
unsigned int j = 0; j < spacedim; ++j)
2171 for (
unsigned int K = 0; K < dim; ++K)
2173 tmp[K] = data.
covariant[q][j][0] * input[q][i][0][K];
2174 for (
unsigned int J = 1; J < dim; ++J)
2175 tmp[K] += data.
covariant[q][j][J] * input[q][i][J][K];
2177 for (
unsigned int k = 0; k < spacedim; ++k)
2179 output[q][i][j][k] = data.
covariant[q][k][0] * tmp[0];
2180 for (
unsigned int K = 1; K < dim; ++K)
2181 output[q][i][j][k] += data.
covariant[q][k][K] * tmp[K];
2194template <
int dim,
int spacedim>
2198 const MappingKind mapping_kind,
2202 switch (mapping_kind)
2207 internal::MappingFEImplementation::transform_hessians(input,
2221 template <
int spacedim>
2223 check_all_manifold_ids_identical(
2231 template <
int spacedim>
2233 check_all_manifold_ids_identical(
2236 const auto m_id = cell->manifold_id();
2238 for (
const auto f : cell->face_indices())
2239 if (m_id != cell->face(f)->manifold_id())
2247 template <
int spacedim>
2249 check_all_manifold_ids_identical(
2252 const auto m_id = cell->manifold_id();
2254 for (
const auto f : cell->face_indices())
2255 if (m_id != cell->face(f)->manifold_id())
2258 for (
const auto l : cell->line_indices())
2259 if (m_id != cell->line(l)->manifold_id())
2268template <
int dim,
int spacedim>
2269std::vector<Point<spacedim>>
2274 check_all_manifold_ids_identical(cell),
2276 "All entities of a cell need to have the same manifold id as the cell has."));
2278 std::vector<Point<spacedim>> vertices(cell->
n_vertices());
2281 vertices[i] = cell->
vertex(i);
2283 std::vector<Point<spacedim>> mapping_support_points(
2284 fe->get_unit_support_points().size());
2288 mapping_support_points);
2290 return mapping_support_points;
2295template <
int dim,
int spacedim>
2305template <
int dim,
int spacedim>
2310 Assert(dim == reference_cell.get_dimension(),
2311 ExcMessage(
"The dimension of your mapping (" +
2313 ") and the reference cell cell_type (" +
2315 " ) do not agree."));
2317 return fe->reference_cell() == reference_cell;
2323#include "mapping_fe.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
double diameter(const Mapping< dim, spacedim > &mapping) const
bool direction_flag() const
::internal::SubfaceCase< dim > subface_case(const unsigned int face_no) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< Point< spacedim > > mapping_support_points
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< std::vector< Tensor< 1, spacedim > > > aux
std::vector< Tensor< 3, dim > > shape_third_derivatives
std::vector< Tensor< 2, dim > > shape_second_derivatives
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
void compute_shape_function_values(const std::vector< Point< dim > > &unit_points)
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
std::vector< double > quadrature_weights
const FiniteElement< dim, spacedim > & fe
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< double > shape_values
std::vector< Tensor< 1, dim > > shape_derivatives
virtual std::size_t memory_consumption() const override
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
InternalData(const FiniteElement< dim, spacedim > &fe)
const unsigned int n_shape_functions
std::vector< double > volume_elements
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
const unsigned int polynomial_degree
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
const unsigned int polynomial_degree
MappingFE(const FiniteElement< dim, spacedim > &fe)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Table< 2, double > mapping_support_point_weights
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual void transform(const ArrayView< const Tensor< 1, dim > > &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim > > &output) const override
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
unsigned int get_degree() const
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
const std::unique_ptr< FiniteElement< dim, spacedim > > fe
virtual std::size_t memory_consumption() const
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static DataSetDescriptor cell()
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
static Quadrature< dim > project_to_all_subfaces(const ReferenceCell &reference_cell, const SubQuadrature &quadrature)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
const Triangulation< dim, spacedim > & get_triangulation() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
const Manifold< dim, spacedim > & get_manifold() const
unsigned int n_vertices() const
Point< spacedim > center(const bool respect_manifold=false, const bool interpolate_from_surrounding=false) const
Point< spacedim > & vertex(const unsigned int i) const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int size() const
unsigned int max_n_quadrature_points() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
constexpr ndarray< Tensor< 1, dim >, GeometryInfo< dim >::faces_per_cell, dim - 1 > GeometryInfo< dim >::unit_tangential_vectors
static ::ExceptionBase & ExcAccessToUninitializedField(std::string arg1)
static ::ExceptionBase & ExcDistortedMappedCell(Point< spacedim > arg1, double arg2, int arg3)
static ::ExceptionBase & ExcTransformationFailed()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_volume_elements
Determinant of the Jacobian.
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_jacobian_3rd_derivatives
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
@ mapping_covariant_gradient
@ mapping_contravariant_hessian
@ mapping_covariant_hessian
@ mapping_contravariant_gradient
unsigned char combined_face_orientation(const unsigned int face) const
#define DEAL_II_NOT_IMPLEMENTED()
MappingQ< dim, spacedim > StaticMappingQ1< dim, spacedim >::mapping
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
constexpr T fixed_power(const T t)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void transform_gradients(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void transform_hessians(const ArrayView< const Tensor< 3, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim > > &output)
void maybe_compute_face_data(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int n_q_points, const std::vector< double > &weights, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)