35#include <boost/container/small_vector.hpp>
47template <
int dim,
int spacedim>
51 , polynomial_degree(fe.tensor_degree())
52 , n_shape_functions(fe.n_dofs_per_cell())
57template <
int dim,
int spacedim>
77template <
int dim,
int spacedim>
82 const unsigned int n_original_q_points)
86 this->update_each = update_flags;
88 const unsigned int n_q_points = q.
size();
91 covariant.resize(n_original_q_points);
94 contravariant.resize(n_original_q_points);
97 volume_elements.resize(n_original_q_points);
102 shape_values.resize(n_shape_functions * n_q_points);
104 if (this->update_each &
112 shape_derivatives.resize(n_shape_functions * n_q_points);
114 if (this->update_each &
116 shape_second_derivatives.resize(n_shape_functions * n_q_points);
120 shape_third_derivatives.resize(n_shape_functions * n_q_points);
124 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
127 compute_shape_function_values(q.
get_points());
135template <
int dim,
int spacedim>
140 const unsigned int n_original_q_points)
142 initialize(update_flags, q, n_original_q_points);
144 if (this->update_each &
152 const auto reference_cell = this->
fe.reference_cell();
153 const auto n_faces = reference_cell.n_faces();
155 for (
unsigned int i = 0; i < n_faces; ++i)
157 unit_tangentials[i].resize(n_original_q_points);
158 std::fill(unit_tangentials[i].begin(),
159 unit_tangentials[i].end(),
160 reference_cell.template unit_tangential_vectors<dim>(i, 0));
163 unit_tangentials[n_faces + i].resize(n_original_q_points);
165 unit_tangentials[n_faces + i].begin(),
166 unit_tangentials[n_faces + i].end(),
167 reference_cell.template unit_tangential_vectors<dim>(i, 1));
175template <
int dim,
int spacedim>
184 const auto &tensor_pols = fe_poly->get_poly_space();
186 const unsigned int n_shape_functions =
fe.n_dofs_per_cell();
187 const unsigned int n_points = unit_points.size();
189 std::vector<double> values;
190 std::vector<Tensor<1, dim>> grads;
191 if (shape_values.size() != 0)
193 Assert(shape_values.size() == n_shape_functions * n_points,
195 values.resize(n_shape_functions);
197 if (shape_derivatives.size() != 0)
199 Assert(shape_derivatives.size() == n_shape_functions * n_points,
201 grads.resize(n_shape_functions);
204 std::vector<Tensor<2, dim>> grad2;
205 if (shape_second_derivatives.size() != 0)
207 Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
209 grad2.resize(n_shape_functions);
212 std::vector<Tensor<3, dim>> grad3;
213 if (shape_third_derivatives.size() != 0)
215 Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
217 grad3.resize(n_shape_functions);
220 std::vector<Tensor<4, dim>> grad4;
221 if (shape_fourth_derivatives.size() != 0)
223 Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
225 grad4.resize(n_shape_functions);
229 if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
230 shape_second_derivatives.size() != 0 ||
231 shape_third_derivatives.size() != 0 ||
232 shape_fourth_derivatives.size() != 0)
233 for (
unsigned int point = 0; point < n_points; ++point)
235 tensor_pols.evaluate(
236 unit_points[point], values, grads, grad2, grad3, grad4);
238 if (shape_values.size() != 0)
239 for (
unsigned int i = 0; i < n_shape_functions; ++i)
240 shape(point, i) = values[i];
242 if (shape_derivatives.size() != 0)
243 for (
unsigned int i = 0; i < n_shape_functions; ++i)
244 derivative(point, i) = grads[i];
246 if (shape_second_derivatives.size() != 0)
247 for (
unsigned int i = 0; i < n_shape_functions; ++i)
248 second_derivative(point, i) = grad2[i];
250 if (shape_third_derivatives.size() != 0)
251 for (
unsigned int i = 0; i < n_shape_functions; ++i)
252 third_derivative(point, i) = grad3[i];
254 if (shape_fourth_derivatives.size() != 0)
255 for (
unsigned int i = 0; i < n_shape_functions; ++i)
256 fourth_derivative(point, i) = grad4[i];
263 namespace MappingFEImplementation
273 template <
int dim,
int spacedim>
275 maybe_compute_q_points(
277 const typename ::MappingFE<dim, spacedim>::InternalData &data,
279 const unsigned int n_q_points)
284 for (
unsigned int point = 0; point < n_q_points; ++point)
286 const double *shape = &data.shape(point + data_set, 0);
288 (shape[0] * data.mapping_support_points[0]);
289 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
290 for (
unsigned int i = 0; i < spacedim; ++i)
291 result[i] += shape[k] * data.mapping_support_points[k][i];
292 quadrature_points[point] = result;
306 template <
int dim,
int spacedim>
308 maybe_update_Jacobians(
310 const typename ::QProjector<dim>::DataSetDescriptor data_set,
311 const typename ::MappingFE<dim, spacedim>::InternalData &data,
312 const unsigned int n_q_points)
322 std::fill(data.contravariant.begin(),
323 data.contravariant.end(),
328 for (
unsigned int point = 0; point < n_q_points; ++point)
330 double result[spacedim][dim];
334 for (
unsigned int i = 0; i < spacedim; ++i)
335 for (
unsigned int j = 0; j < dim; ++j)
336 result[i][j] = data.derivative(point + data_set, 0)[j] *
337 data.mapping_support_points[0][i];
338 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
339 for (
unsigned int i = 0; i < spacedim; ++i)
340 for (
unsigned int j = 0; j < dim; ++j)
342 data.derivative(point + data_set, k)[j] *
343 data.mapping_support_points[k][i];
350 for (
unsigned int i = 0; i < spacedim; ++i)
351 for (
unsigned int j = 0; j < dim; ++j)
352 data.contravariant[point][i][j] = result[i][j];
359 for (
unsigned int point = 0; point < n_q_points; ++point)
361 data.covariant[point] =
362 (data.contravariant[point]).covariant_form();
369 for (
unsigned int point = 0; point < n_q_points; ++point)
370 data.volume_elements[point] =
371 data.contravariant[point].determinant();
381 template <
int dim,
int spacedim>
383 maybe_update_jacobian_grads(
386 const typename ::MappingFE<dim, spacedim>::InternalData &data,
388 const unsigned int n_q_points)
396 for (
unsigned int point = 0; point < n_q_points; ++point)
399 &data.second_derivative(point + data_set, 0);
400 double result[spacedim][dim][dim];
401 for (
unsigned int i = 0; i < spacedim; ++i)
402 for (
unsigned int j = 0; j < dim; ++j)
403 for (
unsigned int l = 0; l < dim; ++l)
405 (
second[0][j][l] * data.mapping_support_points[0][i]);
406 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
407 for (
unsigned int i = 0; i < spacedim; ++i)
408 for (
unsigned int j = 0; j < dim; ++j)
409 for (
unsigned int l = 0; l < dim; ++l)
412 data.mapping_support_points[k][i]);
414 for (
unsigned int i = 0; i < spacedim; ++i)
415 for (
unsigned int j = 0; j < dim; ++j)
416 for (
unsigned int l = 0; l < dim; ++l)
417 jacobian_grads[point][i][j][l] = result[i][j][l];
428 template <
int dim,
int spacedim>
430 maybe_update_jacobian_pushed_forward_grads(
433 const typename ::MappingFE<dim, spacedim>::InternalData &data,
435 const unsigned int n_q_points)
441 jacobian_pushed_forward_grads.size() + 1);
445 double tmp[spacedim][spacedim][spacedim];
446 for (
unsigned int point = 0; point < n_q_points; ++point)
449 &data.second_derivative(point + data_set, 0);
450 double result[spacedim][dim][dim];
451 for (
unsigned int i = 0; i < spacedim; ++i)
452 for (
unsigned int j = 0; j < dim; ++j)
453 for (
unsigned int l = 0; l < dim; ++l)
454 result[i][j][l] = (
second[0][j][l] *
455 data.mapping_support_points[0][i]);
456 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
457 for (
unsigned int i = 0; i < spacedim; ++i)
458 for (
unsigned int j = 0; j < dim; ++j)
459 for (
unsigned int l = 0; l < dim; ++l)
462 data.mapping_support_points[k][i]);
465 for (
unsigned int i = 0; i < spacedim; ++i)
466 for (
unsigned int j = 0; j < spacedim; ++j)
467 for (
unsigned int l = 0; l < dim; ++l)
470 result[i][0][l] * data.covariant[point][j][0];
471 for (
unsigned int jr = 1; jr < dim; ++jr)
473 tmp[i][j][l] += result[i][jr][l] *
474 data.covariant[point][j][jr];
479 for (
unsigned int i = 0; i < spacedim; ++i)
480 for (
unsigned int j = 0; j < spacedim; ++j)
481 for (
unsigned int l = 0; l < spacedim; ++l)
483 jacobian_pushed_forward_grads[point][i][j][l] =
484 tmp[i][j][0] * data.covariant[point][l][0];
485 for (
unsigned int lr = 1; lr < dim; ++lr)
487 jacobian_pushed_forward_grads[point][i][j][l] +=
488 tmp[i][j][lr] * data.covariant[point][l][lr];
502 template <
int dim,
int spacedim>
504 maybe_update_jacobian_2nd_derivatives(
507 const typename ::MappingFE<dim, spacedim>::InternalData &data,
509 const unsigned int n_q_points)
518 for (
unsigned int point = 0; point < n_q_points; ++point)
521 &data.third_derivative(point + data_set, 0);
522 double result[spacedim][dim][dim][dim];
523 for (
unsigned int i = 0; i < spacedim; ++i)
524 for (
unsigned int j = 0; j < dim; ++j)
525 for (
unsigned int l = 0; l < dim; ++l)
526 for (
unsigned int m = 0; m < dim; ++m)
529 data.mapping_support_points[0][i]);
530 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
531 for (
unsigned int i = 0; i < spacedim; ++i)
532 for (
unsigned int j = 0; j < dim; ++j)
533 for (
unsigned int l = 0; l < dim; ++l)
534 for (
unsigned int m = 0; m < dim; ++m)
535 result[i][j][l][m] +=
537 data.mapping_support_points[k][i]);
539 for (
unsigned int i = 0; i < spacedim; ++i)
540 for (
unsigned int j = 0; j < dim; ++j)
541 for (
unsigned int l = 0; l < dim; ++l)
542 for (
unsigned int m = 0; m < dim; ++m)
543 jacobian_2nd_derivatives[point][i][j][l][m] =
557 template <
int dim,
int spacedim>
559 maybe_update_jacobian_pushed_forward_2nd_derivatives(
562 const typename ::MappingFE<dim, spacedim>::InternalData &data,
564 &jacobian_pushed_forward_2nd_derivatives,
565 const unsigned int n_q_points)
571 jacobian_pushed_forward_2nd_derivatives.size() +
576 double tmp[spacedim][spacedim][spacedim][spacedim];
577 for (
unsigned int point = 0; point < n_q_points; ++point)
580 &data.third_derivative(point + data_set, 0);
581 double result[spacedim][dim][dim][dim];
582 for (
unsigned int i = 0; i < spacedim; ++i)
583 for (
unsigned int j = 0; j < dim; ++j)
584 for (
unsigned int l = 0; l < dim; ++l)
585 for (
unsigned int m = 0; m < dim; ++m)
588 data.mapping_support_points[0][i]);
589 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
590 for (
unsigned int i = 0; i < spacedim; ++i)
591 for (
unsigned int j = 0; j < dim; ++j)
592 for (
unsigned int l = 0; l < dim; ++l)
593 for (
unsigned int m = 0; m < dim; ++m)
594 result[i][j][l][m] +=
596 data.mapping_support_points[k][i]);
599 for (
unsigned int i = 0; i < spacedim; ++i)
600 for (
unsigned int j = 0; j < spacedim; ++j)
601 for (
unsigned int l = 0; l < dim; ++l)
602 for (
unsigned int m = 0; m < dim; ++m)
604 jacobian_pushed_forward_2nd_derivatives
605 [point][i][j][l][m] =
607 data.covariant[point][j][0];
608 for (
unsigned int jr = 1; jr < dim; ++jr)
609 jacobian_pushed_forward_2nd_derivatives[point]
612 result[i][jr][l][m] *
613 data.covariant[point][j][jr];
617 for (
unsigned int i = 0; i < spacedim; ++i)
618 for (
unsigned int j = 0; j < spacedim; ++j)
619 for (
unsigned int l = 0; l < spacedim; ++l)
620 for (
unsigned int m = 0; m < dim; ++m)
623 jacobian_pushed_forward_2nd_derivatives[point]
626 data.covariant[point][l][0];
627 for (
unsigned int lr = 1; lr < dim; ++lr)
629 jacobian_pushed_forward_2nd_derivatives
630 [point][i][j][lr][m] *
631 data.covariant[point][l][lr];
635 for (
unsigned int i = 0; i < spacedim; ++i)
636 for (
unsigned int j = 0; j < spacedim; ++j)
637 for (
unsigned int l = 0; l < spacedim; ++l)
638 for (
unsigned int m = 0; m < spacedim; ++m)
640 jacobian_pushed_forward_2nd_derivatives
641 [point][i][j][l][m] =
642 tmp[i][j][l][0] * data.covariant[point][m][0];
643 for (
unsigned int mr = 1; mr < dim; ++mr)
644 jacobian_pushed_forward_2nd_derivatives[point]
648 data.covariant[point][m][mr];
661 template <
int dim,
int spacedim>
663 maybe_update_jacobian_3rd_derivatives(
666 const typename ::MappingFE<dim, spacedim>::InternalData &data,
668 const unsigned int n_q_points)
677 for (
unsigned int point = 0; point < n_q_points; ++point)
680 &data.fourth_derivative(point + data_set, 0);
681 double result[spacedim][dim][dim][dim][dim];
682 for (
unsigned int i = 0; i < spacedim; ++i)
683 for (
unsigned int j = 0; j < dim; ++j)
684 for (
unsigned int l = 0; l < dim; ++l)
685 for (
unsigned int m = 0; m < dim; ++m)
686 for (
unsigned int n = 0; n < dim; ++n)
687 result[i][j][l][m][n] =
688 (fourth[0][j][l][m][n] *
689 data.mapping_support_points[0][i]);
690 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
691 for (
unsigned int i = 0; i < spacedim; ++i)
692 for (
unsigned int j = 0; j < dim; ++j)
693 for (
unsigned int l = 0; l < dim; ++l)
694 for (
unsigned int m = 0; m < dim; ++m)
695 for (
unsigned int n = 0; n < dim; ++n)
696 result[i][j][l][m][n] +=
697 (fourth[k][j][l][m][n] *
698 data.mapping_support_points[k][i]);
700 for (
unsigned int i = 0; i < spacedim; ++i)
701 for (
unsigned int j = 0; j < dim; ++j)
702 for (
unsigned int l = 0; l < dim; ++l)
703 for (
unsigned int m = 0; m < dim; ++m)
704 for (
unsigned int n = 0; n < dim; ++n)
705 jacobian_3rd_derivatives[point][i][j][l][m][n] =
706 result[i][j][l][m][n];
719 template <
int dim,
int spacedim>
721 maybe_update_jacobian_pushed_forward_3rd_derivatives(
724 const typename ::MappingFE<dim, spacedim>::InternalData &data,
726 &jacobian_pushed_forward_3rd_derivatives,
727 const unsigned int n_q_points)
733 jacobian_pushed_forward_3rd_derivatives.size() +
738 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
739 for (
unsigned int point = 0; point < n_q_points; ++point)
742 &data.fourth_derivative(point + data_set, 0);
743 double result[spacedim][dim][dim][dim][dim];
744 for (
unsigned int i = 0; i < spacedim; ++i)
745 for (
unsigned int j = 0; j < dim; ++j)
746 for (
unsigned int l = 0; l < dim; ++l)
747 for (
unsigned int m = 0; m < dim; ++m)
748 for (
unsigned int n = 0; n < dim; ++n)
749 result[i][j][l][m][n] =
750 (fourth[0][j][l][m][n] *
751 data.mapping_support_points[0][i]);
752 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
753 for (
unsigned int i = 0; i < spacedim; ++i)
754 for (
unsigned int j = 0; j < dim; ++j)
755 for (
unsigned int l = 0; l < dim; ++l)
756 for (
unsigned int m = 0; m < dim; ++m)
757 for (
unsigned int n = 0; n < dim; ++n)
758 result[i][j][l][m][n] +=
759 (fourth[k][j][l][m][n] *
760 data.mapping_support_points[k][i]);
763 for (
unsigned int i = 0; i < spacedim; ++i)
764 for (
unsigned int j = 0; j < spacedim; ++j)
765 for (
unsigned int l = 0; l < dim; ++l)
766 for (
unsigned int m = 0; m < dim; ++m)
767 for (
unsigned int n = 0; n < dim; ++n)
770 result[i][0][l][m][n] *
771 data.covariant[point][j][0];
772 for (
unsigned int jr = 1; jr < dim; ++jr)
773 tmp[i][j][l][m][n] +=
774 result[i][jr][l][m][n] *
775 data.covariant[point][j][jr];
779 for (
unsigned int i = 0; i < spacedim; ++i)
780 for (
unsigned int j = 0; j < spacedim; ++j)
781 for (
unsigned int l = 0; l < spacedim; ++l)
782 for (
unsigned int m = 0; m < dim; ++m)
783 for (
unsigned int n = 0; n < dim; ++n)
785 jacobian_pushed_forward_3rd_derivatives
786 [point][i][j][l][m][n] =
788 data.covariant[point][l][0];
789 for (
unsigned int lr = 1; lr < dim; ++lr)
790 jacobian_pushed_forward_3rd_derivatives
791 [point][i][j][l][m][n] +=
792 tmp[i][j][lr][m][n] *
793 data.covariant[point][l][lr];
797 for (
unsigned int i = 0; i < spacedim; ++i)
798 for (
unsigned int j = 0; j < spacedim; ++j)
799 for (
unsigned int l = 0; l < spacedim; ++l)
800 for (
unsigned int m = 0; m < spacedim; ++m)
801 for (
unsigned int n = 0; n < dim; ++n)
804 jacobian_pushed_forward_3rd_derivatives
805 [point][i][j][l][0][n] *
806 data.covariant[point][m][0];
807 for (
unsigned int mr = 1; mr < dim; ++mr)
808 tmp[i][j][l][m][n] +=
809 jacobian_pushed_forward_3rd_derivatives
810 [point][i][j][l][mr][n] *
811 data.covariant[point][m][mr];
815 for (
unsigned int i = 0; i < spacedim; ++i)
816 for (
unsigned int j = 0; j < spacedim; ++j)
817 for (
unsigned int l = 0; l < spacedim; ++l)
818 for (
unsigned int m = 0; m < spacedim; ++m)
819 for (
unsigned int n = 0; n < spacedim; ++n)
821 jacobian_pushed_forward_3rd_derivatives
822 [point][i][j][l][m][n] =
824 data.covariant[point][n][0];
825 for (
unsigned int nr = 1; nr < dim; ++nr)
826 jacobian_pushed_forward_3rd_derivatives
827 [point][i][j][l][m][n] +=
828 tmp[i][j][l][m][nr] *
829 data.covariant[point][n][nr];
841template <
int dim,
int spacedim>
847 ExcMessage(
"It only makes sense to create polynomial mappings "
848 "with a polynomial degree greater or equal to one."));
853 const auto &mapping_support_points =
fe.get_unit_support_points();
855 const auto reference_cell =
fe.reference_cell();
857 const unsigned int n_points = mapping_support_points.size();
858 const unsigned int n_shape_functions = reference_cell.n_vertices();
863 for (
unsigned int point = 0; point < n_points; ++point)
864 for (
unsigned int i = 0; i < n_shape_functions; ++i)
866 reference_cell.d_linear_shape_function(mapping_support_points[point],
872template <
int dim,
int spacedim>
874 : fe(mapping.fe->clone())
875 , polynomial_degree(mapping.polynomial_degree)
876 , mapping_support_point_weights(mapping.mapping_support_point_weights)
881template <
int dim,
int spacedim>
882std::unique_ptr<Mapping<dim, spacedim>>
885 return std::make_unique<MappingFE<dim, spacedim>>(*this);
890template <
int dim,
int spacedim>
894 return polynomial_degree;
899template <
int dim,
int spacedim>
905 const std::vector<Point<spacedim>> support_points =
906 this->compute_mapping_support_points(cell);
910 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
911 mapped_point += support_points[i] * this->fe->shape_value(i, p);
918template <
int dim,
int spacedim>
924 const std::vector<Point<spacedim>> support_points =
925 this->compute_mapping_support_points(cell);
927 const double eps = 1.e-12 * cell->diameter();
928 const unsigned int loop_limit = 10;
932 unsigned int loop = 0;
948 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
950 mapped_point += support_points[i] * this->fe->shape_value(i, p_unit);
951 const auto grad_F_i = this->fe->shape_grad(i, p_unit);
952 const auto hessian_F_i = this->fe->shape_grad_grad(i, p_unit);
953 for (
unsigned int j = 0; j < dim; ++j)
955 grad_FT[j] += grad_F_i[j] * support_points[i];
956 for (
unsigned int l = 0; l < dim; ++l)
957 hess_FT[j][l] += hessian_F_i[j][l] * support_points[i];
962 const auto residual = p - mapped_point;
969 if (grad_FT_residual.norm() <= eps)
974 for (
unsigned int j = 0; j < dim; ++j)
975 for (
unsigned int l = 0; l < dim; ++l)
976 corrected_metric_tensor[j][l] =
977 -grad_FT[j] * grad_FT[l] + hess_FT[j][l] * residual;
980 const auto g_inverse =
invert(corrected_metric_tensor);
981 p_unit -=
Point<dim>(g_inverse * grad_FT_residual);
985 while (loop < loop_limit);
999template <
int dim,
int spacedim>
1010 for (
unsigned int i = 0; i < 5; ++i)
1055template <
int dim,
int spacedim>
1056std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1060 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1061 std::make_unique<InternalData>(*this->fe);
1063 data.
initialize(this->requires_update_flags(update_flags), q, q.
size());
1070template <
int dim,
int spacedim>
1071std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1076 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1077 std::make_unique<InternalData>(*this->fe);
1081 this->fe->reference_cell(), quadrature),
1089template <
int dim,
int spacedim>
1090std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1095 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1096 std::make_unique<InternalData>(*this->fe);
1100 this->fe->reference_cell(), quadrature),
1108template <
int dim,
int spacedim>
1123 const unsigned int n_q_points = quadrature.
size();
1145 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1151 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1152 computed_cell_similarity,
1157 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1158 computed_cell_similarity,
1164 internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
1166 spacedim>(computed_cell_similarity,
1172 internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
1174 spacedim>(computed_cell_similarity,
1180 internal::MappingFEImplementation::
1181 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1182 computed_cell_similarity,
1188 internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
1190 spacedim>(computed_cell_similarity,
1196 internal::MappingFEImplementation::
1197 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1198 computed_cell_similarity,
1204 const UpdateFlags update_flags = data.update_each;
1205 const std::vector<double> &weights = quadrature.
get_weights();
1221 for (
unsigned int point = 0; point < n_q_points; ++point)
1223 if (dim == spacedim)
1233 1e-12 * Utilities::fixed_power<dim>(
1234 cell->diameter() /
std::sqrt(
double(dim))),
1236 cell->center(), det, point)));
1238 output_data.
JxW_values[point] = weights[point] * det;
1246 for (
unsigned int i = 0; i < spacedim; ++i)
1247 for (
unsigned int j = 0; j < dim; ++j)
1251 for (
unsigned int i = 0; i < dim; ++i)
1252 for (
unsigned int j = 0; j < dim; ++j)
1253 G[i][j] = DX_t[i] * DX_t[j];
1258 if (computed_cell_similarity ==
1269 Assert(spacedim == dim + 1,
1271 "There is no (unique) cell normal for " +
1273 "-dimensional cells in " +
1275 "-dimensional space. This only works if the "
1276 "space dimension is one greater than the "
1277 "dimensionality of the mesh cells."));
1281 cross_product_2d(-DX_t[0]);
1284 cross_product_3d(DX_t[0], DX_t[1]);
1289 if (cell->direction_flag() ==
false)
1304 for (
unsigned int point = 0; point < n_q_points; ++point)
1313 for (
unsigned int point = 0; point < n_q_points; ++point)
1318 return computed_cell_similarity;
1325 namespace MappingFEImplementation
1339 template <
int dim,
int spacedim>
1341 maybe_compute_face_data(
1342 const ::MappingFE<dim, spacedim> &mapping,
1343 const typename ::Triangulation<dim, spacedim>::cell_iterator
1345 const unsigned int face_no,
1346 const unsigned int subface_no,
1347 const unsigned int n_q_points,
1349 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1353 const UpdateFlags update_flags = data.update_each;
1376 for (
unsigned int d = 0; d != dim - 1; ++d)
1378 Assert(face_no + cell->n_faces() * d <
1379 data.unit_tangentials.size(),
1382 data.aux[d].size() <=
1383 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1388 data.unit_tangentials[face_no + cell->n_faces() * d]),
1399 if (dim == spacedim)
1401 for (
unsigned int i = 0; i < n_q_points; ++i)
1411 (face_no == 0 ? -1 : +1);
1415 cross_product_2d(data.aux[0][i]);
1419 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1435 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1441 data.contravariant[
point].transpose()[0];
1443 (face_no == 0 ? -1. : +1.) *
1453 cross_product_3d(DX_t[0], DX_t[1]);
1454 cell_normal /= cell_normal.
norm();
1459 cross_product_3d(data.aux[0][point], cell_normal);
1466 for (
unsigned int i = 0; i < n_q_points; ++i)
1470 data.quadrature_weights[i + data_set];
1475 const double area_ratio =
1477 cell->subface_case(face_no), subface_no);
1486 for (
unsigned int i = 0; i < n_q_points; ++i)
1492 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1493 output_data.
jacobians[point] = data.contravariant[point];
1496 for (
unsigned int point = 0;
point < n_q_points; ++
point)
1498 data.covariant[point].transpose();
1509 template <
int dim,
int spacedim>
1512 const ::MappingFE<dim, spacedim> &mapping,
1513 const typename ::Triangulation<dim, spacedim>::cell_iterator
1515 const unsigned int face_no,
1516 const unsigned int subface_no,
1519 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1523 const unsigned int n_q_points = quadrature.
size();
1525 maybe_compute_q_points<dim, spacedim>(data_set,
1538 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1544 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1550 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1556 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1562 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1584template <
int dim,
int spacedim>
1588 const unsigned int face_no,
1604 (&cell->get_triangulation() !=
1612 internal::MappingFEImplementation::do_fill_fe_face_values(
1619 cell->face_orientation(face_no),
1620 cell->face_flip(face_no),
1621 cell->face_rotation(face_no),
1623 quadrature[quadrature.
size() == 1 ? 0 : face_no],
1630template <
int dim,
int spacedim>
1634 const unsigned int face_no,
1635 const unsigned int subface_no,
1651 (&cell->get_triangulation() !=
1659 internal::MappingFEImplementation::do_fill_fe_face_values(
1667 cell->face_orientation(face_no),
1668 cell->face_flip(face_no),
1669 cell->face_rotation(face_no),
1671 cell->subface_case(face_no)),
1681 namespace MappingFEImplementation
1685 template <
int dim,
int spacedim,
int rank>
1702 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1703 &mapping_data) !=
nullptr),
1705 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1707 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1710 switch (mapping_kind)
1717 "update_contravariant_transformation"));
1719 for (
unsigned int i = 0; i < input.size(); ++i)
1731 "update_contravariant_transformation"));
1735 "update_volume_elements"));
1740 for (
unsigned int i = 0; i < input.size(); ++i)
1744 output[i] /= data.volume_elements[i];
1756 "update_covariant_transformation"));
1758 for (
unsigned int i = 0; i < input.size(); ++i)
1770 template <
int dim,
int spacedim,
int rank>
1781 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1782 &mapping_data) !=
nullptr),
1784 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1786 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1789 switch (mapping_kind)
1796 "update_covariant_transformation"));
1800 "update_contravariant_transformation"));
1803 for (
unsigned int i = 0; i < output.size(); ++i)
1820 "update_covariant_transformation"));
1823 for (
unsigned int i = 0; i < output.size(); ++i)
1840 "update_covariant_transformation"));
1844 "update_contravariant_transformation"));
1848 "update_volume_elements"));
1851 for (
unsigned int i = 0; i < output.size(); ++i)
1860 output[i] /= data.volume_elements[i];
1873 template <
int dim,
int spacedim>
1884 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1885 &mapping_data) !=
nullptr),
1887 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1889 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1892 switch (mapping_kind)
1899 "update_covariant_transformation"));
1903 "update_contravariant_transformation"));
1905 for (
unsigned int q = 0; q < output.size(); ++q)
1906 for (
unsigned int i = 0; i < spacedim; ++i)
1908 double tmp1[dim][dim];
1909 for (
unsigned int J = 0; J < dim; ++J)
1910 for (
unsigned int K = 0;
K < dim; ++
K)
1913 data.contravariant[q][i][0] * input[q][0][J][
K];
1914 for (
unsigned int I = 1; I < dim; ++I)
1916 data.contravariant[q][i][I] * input[q][I][J][K];
1918 for (
unsigned int j = 0; j < spacedim; ++j)
1921 for (
unsigned int K = 0;
K < dim; ++
K)
1923 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1924 for (
unsigned int J = 1; J < dim; ++J)
1925 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1927 for (
unsigned int k = 0; k < spacedim; ++k)
1929 output[q][i][j][k] =
1930 data.covariant[q][k][0] * tmp2[0];
1931 for (
unsigned int K = 1;
K < dim; ++
K)
1932 output[q][i][j][k] +=
1933 data.covariant[q][k][K] * tmp2[K];
1945 "update_covariant_transformation"));
1947 for (
unsigned int q = 0; q < output.size(); ++q)
1948 for (
unsigned int i = 0; i < spacedim; ++i)
1950 double tmp1[dim][dim];
1951 for (
unsigned int J = 0; J < dim; ++J)
1952 for (
unsigned int K = 0;
K < dim; ++
K)
1955 data.covariant[q][i][0] * input[q][0][J][
K];
1956 for (
unsigned int I = 1; I < dim; ++I)
1958 data.covariant[q][i][I] * input[q][I][J][K];
1960 for (
unsigned int j = 0; j < spacedim; ++j)
1963 for (
unsigned int K = 0;
K < dim; ++
K)
1965 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1966 for (
unsigned int J = 1; J < dim; ++J)
1967 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1969 for (
unsigned int k = 0; k < spacedim; ++k)
1971 output[q][i][j][k] =
1972 data.covariant[q][k][0] * tmp2[0];
1973 for (
unsigned int K = 1;
K < dim; ++
K)
1974 output[q][i][j][k] +=
1975 data.covariant[q][k][K] * tmp2[K];
1988 "update_covariant_transformation"));
1992 "update_contravariant_transformation"));
1996 "update_volume_elements"));
1998 for (
unsigned int q = 0; q < output.size(); ++q)
1999 for (
unsigned int i = 0; i < spacedim; ++i)
2002 for (
unsigned int I = 0; I < dim; ++I)
2004 data.contravariant[q][i][I] / data.volume_elements[q];
2005 double tmp1[dim][dim];
2006 for (
unsigned int J = 0; J < dim; ++J)
2007 for (
unsigned int K = 0;
K < dim; ++
K)
2009 tmp1[J][
K] = factor[0] * input[q][0][J][
K];
2010 for (
unsigned int I = 1; I < dim; ++I)
2011 tmp1[J][K] += factor[I] * input[q][I][J][K];
2013 for (
unsigned int j = 0; j < spacedim; ++j)
2016 for (
unsigned int K = 0;
K < dim; ++
K)
2018 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
2019 for (
unsigned int J = 1; J < dim; ++J)
2020 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2022 for (
unsigned int k = 0; k < spacedim; ++k)
2024 output[q][i][j][k] =
2025 data.covariant[q][k][0] * tmp2[0];
2026 for (
unsigned int K = 1;
K < dim; ++
K)
2027 output[q][i][j][k] +=
2028 data.covariant[q][k][K] * tmp2[K];
2043 template <
int dim,
int spacedim,
int rank>
2054 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
2055 &mapping_data) !=
nullptr),
2057 const typename ::MappingFE<dim, spacedim>::InternalData &data =
2059 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
2062 switch (mapping_kind)
2069 "update_covariant_transformation"));
2071 for (
unsigned int i = 0; i < output.size(); ++i)
2086template <
int dim,
int spacedim>
2094 internal::MappingFEImplementation::transform_fields(input,
2102template <
int dim,
int spacedim>
2110 internal::MappingFEImplementation::transform_differential_forms(input,
2118template <
int dim,
int spacedim>
2126 switch (mapping_kind)
2129 internal::MappingFEImplementation::transform_fields(input,
2138 internal::MappingFEImplementation::transform_gradients(input,
2150template <
int dim,
int spacedim>
2163 switch (mapping_kind)
2169 "update_covariant_transformation"));
2171 for (
unsigned int q = 0; q < output.size(); ++q)
2172 for (
unsigned int i = 0; i < spacedim; ++i)
2173 for (
unsigned int j = 0; j < spacedim; ++j)
2176 for (
unsigned int K = 0; K < dim; ++K)
2178 tmp[K] = data.
covariant[q][j][0] * input[q][i][0][K];
2179 for (
unsigned int J = 1; J < dim; ++J)
2180 tmp[K] += data.
covariant[q][j][J] * input[q][i][J][K];
2182 for (
unsigned int k = 0; k < spacedim; ++k)
2184 output[q][i][j][k] = data.
covariant[q][k][0] * tmp[0];
2185 for (
unsigned int K = 1; K < dim; ++K)
2186 output[q][i][j][k] += data.
covariant[q][k][K] * tmp[K];
2199template <
int dim,
int spacedim>
2207 switch (mapping_kind)
2212 internal::MappingFEImplementation::transform_hessians(input,
2226 template <
int spacedim>
2228 check_all_manifold_ids_identical(
2236 template <
int spacedim>
2238 check_all_manifold_ids_identical(
2241 const auto m_id = cell->manifold_id();
2243 for (
const auto f : cell->face_indices())
2252 template <
int spacedim>
2254 check_all_manifold_ids_identical(
2257 const auto m_id = cell->manifold_id();
2259 for (
const auto f : cell->face_indices())
2263 for (
const auto l : cell->line_indices())
2273template <
int dim,
int spacedim>
2274std::vector<Point<spacedim>>
2279 check_all_manifold_ids_identical(cell),
2281 "All entities of a cell need to have the same manifold id as the cell has."));
2283 std::vector<Point<spacedim>>
vertices(cell->n_vertices());
2285 for (
const unsigned int i : cell->vertex_indices())
2288 std::vector<Point<spacedim>> mapping_support_points(
2289 fe->get_unit_support_points().size());
2291 cell->get_manifold().get_new_points(
vertices,
2292 mapping_support_point_weights,
2293 mapping_support_points);
2295 return mapping_support_points;
2300template <
int dim,
int spacedim>
2310template <
int dim,
int spacedim>
2315 Assert(dim == reference_cell.get_dimension(),
2316 ExcMessage(
"The dimension of your mapping (" +
2318 ") and the reference cell cell_type (" +
2320 " ) do not agree."));
2322 return fe->reference_cell() == reference_cell;
2328#include "mapping_fe.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
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
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
void compute_shape_function_values(const std::vector< Point< dim > > &unit_points)
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
virtual std::size_t memory_consumption() const override
InternalData(const FiniteElement< dim, spacedim > &fe)
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
Abstract base class for mapping classes.
static DataSetDescriptor cell()
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
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
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)
@ 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
#define DEAL_II_NOT_IMPLEMENTED()
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)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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 > &)