37 #include <boost/container/small_vector.hpp>
50 template <
int dim,
int spacedim>
54 , polynomial_degree(fe.tensor_degree())
55 , n_shape_functions(fe.n_dofs_per_cell())
60 template <
int dim,
int spacedim>
80 template <
int dim,
int spacedim>
85 const unsigned int n_original_q_points)
89 this->update_each = update_flags;
91 const unsigned int n_q_points = q.
size();
94 covariant.resize(n_original_q_points);
97 contravariant.resize(n_original_q_points);
100 volume_elements.resize(n_original_q_points);
105 shape_values.resize(n_shape_functions * n_q_points);
107 if (this->update_each &
115 shape_derivatives.resize(n_shape_functions * n_q_points);
117 if (this->update_each &
119 shape_second_derivatives.resize(n_shape_functions * n_q_points);
123 shape_third_derivatives.resize(n_shape_functions * n_q_points);
127 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
130 compute_shape_function_values(q.
get_points());
138 template <
int dim,
int spacedim>
143 const unsigned int n_original_q_points)
145 initialize(update_flags, q, n_original_q_points);
147 if (this->update_each &
158 for (
unsigned int i = 0; i < n_faces; ++i)
160 unit_tangentials[i].resize(n_original_q_points);
161 std::fill(unit_tangentials[i].
begin(),
162 unit_tangentials[i].
end(),
166 unit_tangentials[n_faces + i].resize(n_original_q_points);
168 unit_tangentials[n_faces + i].
begin(),
169 unit_tangentials[n_faces + i].
end(),
178 template <
int dim,
int spacedim>
187 const auto &tensor_pols = fe_poly->get_poly_space();
189 const unsigned int n_shape_functions =
fe.n_dofs_per_cell();
190 const unsigned int n_points = unit_points.size();
192 std::vector<double>
values;
193 std::vector<Tensor<1, dim>> grads;
194 if (shape_values.size() != 0)
196 Assert(shape_values.size() == n_shape_functions * n_points,
198 values.resize(n_shape_functions);
200 if (shape_derivatives.size() != 0)
202 Assert(shape_derivatives.size() == n_shape_functions * n_points,
204 grads.resize(n_shape_functions);
207 std::vector<Tensor<2, dim>> grad2;
208 if (shape_second_derivatives.size() != 0)
210 Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
212 grad2.resize(n_shape_functions);
215 std::vector<Tensor<3, dim>> grad3;
216 if (shape_third_derivatives.size() != 0)
218 Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
220 grad3.resize(n_shape_functions);
223 std::vector<Tensor<4, dim>> grad4;
224 if (shape_fourth_derivatives.size() != 0)
226 Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
228 grad4.resize(n_shape_functions);
232 if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
233 shape_second_derivatives.size() != 0 ||
234 shape_third_derivatives.size() != 0 ||
235 shape_fourth_derivatives.size() != 0)
238 tensor_pols.evaluate(
239 unit_points[
point],
values, grads, grad2, grad3, grad4);
241 if (shape_values.size() != 0)
242 for (
unsigned int i = 0; i < n_shape_functions; ++i)
245 if (shape_derivatives.size() != 0)
246 for (
unsigned int i = 0; i < n_shape_functions; ++i)
247 derivative(
point, i) = grads[i];
249 if (shape_second_derivatives.size() != 0)
250 for (
unsigned int i = 0; i < n_shape_functions; ++i)
251 second_derivative(
point, i) = grad2[i];
253 if (shape_third_derivatives.size() != 0)
254 for (
unsigned int i = 0; i < n_shape_functions; ++i)
255 third_derivative(
point, i) = grad3[i];
257 if (shape_fourth_derivatives.size() != 0)
258 for (
unsigned int i = 0; i < n_shape_functions; ++i)
259 fourth_derivative(
point, i) = grad4[i];
266 namespace MappingFEImplementation
276 template <
int dim,
int spacedim>
280 const typename ::MappingFE<dim, spacedim>::InternalData &data,
282 const unsigned int n_q_points)
289 const double * shape = &data.shape(
point + data_set, 0);
291 (shape[0] * data.mapping_support_points[0]);
292 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
293 for (
unsigned int i = 0; i < spacedim; ++i)
294 result[i] += shape[k] * data.mapping_support_points[k][i];
309 template <
int dim,
int spacedim>
313 const typename ::QProjector<dim>::DataSetDescriptor data_set,
314 const typename ::MappingFE<dim, spacedim>::InternalData &data,
315 const unsigned int n_q_points)
325 std::fill(data.contravariant.begin(),
326 data.contravariant.end(),
333 double result[spacedim][dim];
337 for (
unsigned int i = 0; i < spacedim; ++i)
338 for (
unsigned int j = 0; j < dim; ++j)
339 result[i][j] = data.derivative(
point + data_set, 0)[j] *
340 data.mapping_support_points[0][i];
341 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
342 for (
unsigned int i = 0; i < spacedim; ++i)
343 for (
unsigned int j = 0; j < dim; ++j)
345 data.derivative(
point + data_set, k)[j] *
346 data.mapping_support_points[k][i];
353 for (
unsigned int i = 0; i < spacedim; ++i)
354 for (
unsigned int j = 0; j < dim; ++j)
355 data.contravariant[
point][i][j] = result[i][j];
364 data.covariant[
point] =
365 (data.contravariant[
point]).covariant_form();
373 data.volume_elements[
point] =
374 data.contravariant[
point].determinant();
384 template <
int dim,
int spacedim>
389 const typename ::MappingFE<dim, spacedim>::InternalData &data,
391 const unsigned int n_q_points)
402 &data.second_derivative(
point + data_set, 0);
403 double result[spacedim][dim][dim];
404 for (
unsigned int i = 0; i < spacedim; ++i)
405 for (
unsigned int j = 0; j < dim; ++j)
406 for (
unsigned int l = 0;
l < dim; ++
l)
408 (
second[0][j][
l] * data.mapping_support_points[0][i]);
409 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
410 for (
unsigned int i = 0; i < spacedim; ++i)
411 for (
unsigned int j = 0; j < dim; ++j)
412 for (
unsigned int l = 0;
l < dim; ++
l)
415 data.mapping_support_points[k][i]);
417 for (
unsigned int i = 0; i < spacedim; ++i)
418 for (
unsigned int j = 0; j < dim; ++j)
419 for (
unsigned int l = 0;
l < dim; ++
l)
420 jacobian_grads[
point][i][j][
l] = result[i][j][
l];
431 template <
int dim,
int spacedim>
436 const typename ::MappingFE<dim, spacedim>::InternalData &data,
438 const unsigned int n_q_points)
444 jacobian_pushed_forward_grads.size() + 1);
448 double tmp[spacedim][spacedim][spacedim];
452 &data.second_derivative(
point + data_set, 0);
453 double result[spacedim][dim][dim];
454 for (
unsigned int i = 0; i < spacedim; ++i)
455 for (
unsigned int j = 0; j < dim; ++j)
456 for (
unsigned int l = 0;
l < dim; ++
l)
457 result[i][j][
l] = (
second[0][j][
l] *
458 data.mapping_support_points[0][i]);
459 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
460 for (
unsigned int i = 0; i < spacedim; ++i)
461 for (
unsigned int j = 0; j < dim; ++j)
462 for (
unsigned int l = 0;
l < dim; ++
l)
465 data.mapping_support_points[k][i]);
468 for (
unsigned int i = 0; i < spacedim; ++i)
469 for (
unsigned int j = 0; j < spacedim; ++j)
470 for (
unsigned int l = 0;
l < dim; ++
l)
473 result[i][0][
l] * data.covariant[
point][j][0];
474 for (
unsigned int jr = 1; jr < dim; ++jr)
476 tmp[i][j][
l] += result[i][jr][
l] *
477 data.covariant[
point][j][jr];
482 for (
unsigned int i = 0; i < spacedim; ++i)
483 for (
unsigned int j = 0; j < spacedim; ++j)
484 for (
unsigned int l = 0;
l < spacedim; ++
l)
486 jacobian_pushed_forward_grads[
point][i][j][
l] =
487 tmp[i][j][0] * data.covariant[
point][
l][0];
488 for (
unsigned int lr = 1; lr < dim; ++lr)
490 jacobian_pushed_forward_grads[
point][i][j][
l] +=
491 tmp[i][j][lr] * data.covariant[
point][
l][lr];
505 template <
int dim,
int spacedim>
510 const typename ::MappingFE<dim, spacedim>::InternalData &data,
512 const unsigned int n_q_points)
524 &data.third_derivative(
point + data_set, 0);
525 double result[spacedim][dim][dim][dim];
526 for (
unsigned int i = 0; i < spacedim; ++i)
527 for (
unsigned int j = 0; j < dim; ++j)
528 for (
unsigned int l = 0;
l < dim; ++
l)
529 for (
unsigned int m = 0; m < dim; ++m)
532 data.mapping_support_points[0][i]);
533 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
534 for (
unsigned int i = 0; i < spacedim; ++i)
535 for (
unsigned int j = 0; j < dim; ++j)
536 for (
unsigned int l = 0;
l < dim; ++
l)
537 for (
unsigned int m = 0; m < dim; ++m)
538 result[i][j][
l][m] +=
540 data.mapping_support_points[k][i]);
542 for (
unsigned int i = 0; i < spacedim; ++i)
543 for (
unsigned int j = 0; j < dim; ++j)
544 for (
unsigned int l = 0;
l < dim; ++
l)
545 for (
unsigned int m = 0; m < dim; ++m)
546 jacobian_2nd_derivatives[
point][i][j][
l][m] =
560 template <
int dim,
int spacedim>
565 const typename ::MappingFE<dim, spacedim>::InternalData &data,
567 & jacobian_pushed_forward_2nd_derivatives,
568 const unsigned int n_q_points)
574 jacobian_pushed_forward_2nd_derivatives.size() +
579 double tmp[spacedim][spacedim][spacedim][spacedim];
583 &data.third_derivative(
point + data_set, 0);
584 double result[spacedim][dim][dim][dim];
585 for (
unsigned int i = 0; i < spacedim; ++i)
586 for (
unsigned int j = 0; j < dim; ++j)
587 for (
unsigned int l = 0;
l < dim; ++
l)
588 for (
unsigned int m = 0; m < dim; ++m)
591 data.mapping_support_points[0][i]);
592 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
593 for (
unsigned int i = 0; i < spacedim; ++i)
594 for (
unsigned int j = 0; j < dim; ++j)
595 for (
unsigned int l = 0;
l < dim; ++
l)
596 for (
unsigned int m = 0; m < dim; ++m)
597 result[i][j][
l][m] +=
599 data.mapping_support_points[k][i]);
602 for (
unsigned int i = 0; i < spacedim; ++i)
603 for (
unsigned int j = 0; j < spacedim; ++j)
604 for (
unsigned int l = 0;
l < dim; ++
l)
605 for (
unsigned int m = 0; m < dim; ++m)
607 jacobian_pushed_forward_2nd_derivatives
610 data.covariant[
point][j][0];
611 for (
unsigned int jr = 1; jr < dim; ++jr)
612 jacobian_pushed_forward_2nd_derivatives[
point]
615 result[i][jr][
l][m] *
616 data.covariant[
point][j][jr];
620 for (
unsigned int i = 0; i < spacedim; ++i)
621 for (
unsigned int j = 0; j < spacedim; ++j)
622 for (
unsigned int l = 0;
l < spacedim; ++
l)
623 for (
unsigned int m = 0; m < dim; ++m)
626 jacobian_pushed_forward_2nd_derivatives[
point]
629 data.covariant[
point][
l][0];
630 for (
unsigned int lr = 1; lr < dim; ++lr)
632 jacobian_pushed_forward_2nd_derivatives
633 [
point][i][j][lr][m] *
634 data.covariant[
point][
l][lr];
638 for (
unsigned int i = 0; i < spacedim; ++i)
639 for (
unsigned int j = 0; j < spacedim; ++j)
640 for (
unsigned int l = 0;
l < spacedim; ++
l)
641 for (
unsigned int m = 0; m < spacedim; ++m)
643 jacobian_pushed_forward_2nd_derivatives
645 tmp[i][j][
l][0] * data.covariant[
point][m][0];
646 for (
unsigned int mr = 1; mr < dim; ++mr)
647 jacobian_pushed_forward_2nd_derivatives[
point]
651 data.covariant[
point][m][mr];
664 template <
int dim,
int spacedim>
669 const typename ::MappingFE<dim, spacedim>::InternalData &data,
671 const unsigned int n_q_points)
683 &data.fourth_derivative(
point + data_set, 0);
684 double result[spacedim][dim][dim][dim][dim];
685 for (
unsigned int i = 0; i < spacedim; ++i)
686 for (
unsigned int j = 0; j < dim; ++j)
687 for (
unsigned int l = 0;
l < dim; ++
l)
688 for (
unsigned int m = 0; m < dim; ++m)
689 for (
unsigned int n = 0; n < dim; ++n)
690 result[i][j][
l][m][n] =
691 (fourth[0][j][
l][m][n] *
692 data.mapping_support_points[0][i]);
693 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
694 for (
unsigned int i = 0; i < spacedim; ++i)
695 for (
unsigned int j = 0; j < dim; ++j)
696 for (
unsigned int l = 0;
l < dim; ++
l)
697 for (
unsigned int m = 0; m < dim; ++m)
698 for (
unsigned int n = 0; n < dim; ++n)
699 result[i][j][
l][m][n] +=
700 (fourth[k][j][
l][m][n] *
701 data.mapping_support_points[k][i]);
703 for (
unsigned int i = 0; i < spacedim; ++i)
704 for (
unsigned int j = 0; j < dim; ++j)
705 for (
unsigned int l = 0;
l < dim; ++
l)
706 for (
unsigned int m = 0; m < dim; ++m)
707 for (
unsigned int n = 0; n < dim; ++n)
708 jacobian_3rd_derivatives[
point][i][j][
l][m][n] =
709 result[i][j][
l][m][n];
722 template <
int dim,
int spacedim>
727 const typename ::MappingFE<dim, spacedim>::InternalData &data,
729 & jacobian_pushed_forward_3rd_derivatives,
730 const unsigned int n_q_points)
736 jacobian_pushed_forward_3rd_derivatives.size() +
741 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
745 &data.fourth_derivative(
point + data_set, 0);
746 double result[spacedim][dim][dim][dim][dim];
747 for (
unsigned int i = 0; i < spacedim; ++i)
748 for (
unsigned int j = 0; j < dim; ++j)
749 for (
unsigned int l = 0;
l < dim; ++
l)
750 for (
unsigned int m = 0; m < dim; ++m)
751 for (
unsigned int n = 0; n < dim; ++n)
752 result[i][j][
l][m][n] =
753 (fourth[0][j][
l][m][n] *
754 data.mapping_support_points[0][i]);
755 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
756 for (
unsigned int i = 0; i < spacedim; ++i)
757 for (
unsigned int j = 0; j < dim; ++j)
758 for (
unsigned int l = 0;
l < dim; ++
l)
759 for (
unsigned int m = 0; m < dim; ++m)
760 for (
unsigned int n = 0; n < dim; ++n)
761 result[i][j][
l][m][n] +=
762 (fourth[k][j][
l][m][n] *
763 data.mapping_support_points[k][i]);
766 for (
unsigned int i = 0; i < spacedim; ++i)
767 for (
unsigned int j = 0; j < spacedim; ++j)
768 for (
unsigned int l = 0;
l < dim; ++
l)
769 for (
unsigned int m = 0; m < dim; ++m)
770 for (
unsigned int n = 0; n < dim; ++n)
773 result[i][0][
l][m][n] *
774 data.covariant[
point][j][0];
775 for (
unsigned int jr = 1; jr < dim; ++jr)
776 tmp[i][j][
l][m][n] +=
777 result[i][jr][
l][m][n] *
778 data.covariant[
point][j][jr];
782 for (
unsigned int i = 0; i < spacedim; ++i)
783 for (
unsigned int j = 0; j < spacedim; ++j)
784 for (
unsigned int l = 0;
l < spacedim; ++
l)
785 for (
unsigned int m = 0; m < dim; ++m)
786 for (
unsigned int n = 0; n < dim; ++n)
788 jacobian_pushed_forward_3rd_derivatives
791 data.covariant[
point][
l][0];
792 for (
unsigned int lr = 1; lr < dim; ++lr)
793 jacobian_pushed_forward_3rd_derivatives
795 tmp[i][j][lr][m][n] *
796 data.covariant[
point][
l][lr];
800 for (
unsigned int i = 0; i < spacedim; ++i)
801 for (
unsigned int j = 0; j < spacedim; ++j)
802 for (
unsigned int l = 0;
l < spacedim; ++
l)
803 for (
unsigned int m = 0; m < spacedim; ++m)
804 for (
unsigned int n = 0; n < dim; ++n)
807 jacobian_pushed_forward_3rd_derivatives
809 data.covariant[
point][m][0];
810 for (
unsigned int mr = 1; mr < dim; ++mr)
811 tmp[i][j][
l][m][n] +=
812 jacobian_pushed_forward_3rd_derivatives
814 data.covariant[
point][m][mr];
818 for (
unsigned int i = 0; i < spacedim; ++i)
819 for (
unsigned int j = 0; j < spacedim; ++j)
820 for (
unsigned int l = 0;
l < spacedim; ++
l)
821 for (
unsigned int m = 0; m < spacedim; ++m)
822 for (
unsigned int n = 0; n < spacedim; ++n)
824 jacobian_pushed_forward_3rd_derivatives
827 data.covariant[
point][n][0];
828 for (
unsigned int nr = 1; nr < dim; ++nr)
829 jacobian_pushed_forward_3rd_derivatives
831 tmp[i][j][
l][m][nr] *
832 data.covariant[
point][n][nr];
844 template <
int dim,
int spacedim>
850 ExcMessage(
"It only makes sense to create polynomial mappings "
851 "with a polynomial degree greater or equal to one."));
856 const auto &mapping_support_points =
fe.get_unit_support_points();
860 const unsigned int n_points = mapping_support_points.size();
861 const unsigned int n_shape_functions =
reference_cell.n_vertices();
867 for (
unsigned int i = 0; i < n_shape_functions; ++i)
875 template <
int dim,
int spacedim>
877 : fe(mapping.fe->clone())
878 , polynomial_degree(mapping.polynomial_degree)
879 , mapping_support_point_weights(mapping.mapping_support_point_weights)
884 template <
int dim,
int spacedim>
885 std::unique_ptr<Mapping<dim, spacedim>>
888 return std::make_unique<MappingFE<dim, spacedim>>(*this);
893 template <
int dim,
int spacedim>
897 return polynomial_degree;
902 template <
int dim,
int spacedim>
908 const std::vector<Point<spacedim>> support_points =
909 this->compute_mapping_support_points(cell);
913 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
914 mapped_point += support_points[i] * this->fe->shape_value(i, p);
921 template <
int dim,
int spacedim>
927 const std::vector<Point<spacedim>> support_points =
928 this->compute_mapping_support_points(cell);
930 const double eps = 1.e-12 * cell->diameter();
931 const unsigned int loop_limit = 10;
935 unsigned int loop = 0;
951 for (
unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
953 mapped_point += support_points[i] * this->fe->shape_value(i, p_unit);
954 const auto grad_F_i = this->fe->shape_grad(i, p_unit);
955 const auto hessian_F_i = this->fe->shape_grad_grad(i, p_unit);
956 for (
unsigned int j = 0; j < dim; ++j)
958 grad_FT[j] += grad_F_i[j] * support_points[i];
959 for (
unsigned int l = 0;
l < dim; ++
l)
960 hess_FT[j][
l] += hessian_F_i[j][
l] * support_points[i];
965 const auto residual = p - mapped_point;
972 if (grad_FT_residual.norm() <=
eps)
977 for (
unsigned int j = 0; j < dim; ++j)
978 for (
unsigned int l = 0;
l < dim; ++
l)
979 corrected_metric_tensor[j][
l] =
980 -grad_FT[j] * grad_FT[
l] + hess_FT[j][
l] * residual;
983 const auto g_inverse =
invert(corrected_metric_tensor);
984 p_unit -=
Point<dim>(g_inverse * grad_FT_residual);
988 while (
loop < loop_limit);
1002 template <
int dim,
int spacedim>
1013 for (
unsigned int i = 0; i < 5; ++i)
1058 template <
int dim,
int spacedim>
1059 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1063 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1064 std::make_unique<InternalData>(*this->fe);
1066 data.
initialize(this->requires_update_flags(update_flags), q, q.
size());
1073 template <
int dim,
int spacedim>
1074 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1079 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1080 std::make_unique<InternalData>(*this->fe);
1084 this->fe->reference_cell(), quadrature),
1092 template <
int dim,
int spacedim>
1093 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1098 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1099 std::make_unique<InternalData>(*this->fe);
1103 this->fe->reference_cell(), quadrature),
1111 template <
int dim,
int spacedim>
1126 const unsigned int n_q_points = quadrature.
size();
1148 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1154 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1155 computed_cell_similarity,
1160 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1161 computed_cell_similarity,
1169 spacedim>(computed_cell_similarity,
1177 spacedim>(computed_cell_similarity,
1183 internal::MappingFEImplementation::
1184 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1185 computed_cell_similarity,
1193 spacedim>(computed_cell_similarity,
1199 internal::MappingFEImplementation::
1200 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1201 computed_cell_similarity,
1207 const UpdateFlags update_flags = data.update_each;
1208 const std::vector<double> &weights = quadrature.
get_weights();
1226 if (dim == spacedim)
1236 1
e-12 * Utilities::fixed_power<dim>(
1237 cell->diameter() / std::sqrt(
double(dim))),
1239 cell->center(), det,
point)));
1249 for (
unsigned int i = 0; i < spacedim; ++i)
1250 for (
unsigned int j = 0; j < dim; ++j)
1254 for (
unsigned int i = 0; i < dim; ++i)
1255 for (
unsigned int j = 0; j < dim; ++j)
1256 G[i][j] = DX_t[i] * DX_t[j];
1261 if (computed_cell_similarity ==
1272 Assert(spacedim == dim + 1,
1274 "There is no (unique) cell normal for " +
1276 "-dimensional cells in " +
1278 "-dimensional space. This only works if the "
1279 "space dimension is one greater than the "
1280 "dimensionality of the mesh cells."));
1284 cross_product_2d(-DX_t[0]);
1287 cross_product_3d(DX_t[0], DX_t[1]);
1292 if (cell->direction_flag() ==
false)
1321 return computed_cell_similarity;
1328 namespace MappingFEImplementation
1342 template <
int dim,
int spacedim>
1345 const ::MappingFE<dim, spacedim> &mapping,
1346 const typename ::Triangulation<dim, spacedim>::cell_iterator
1348 const unsigned int face_no,
1349 const unsigned int subface_no,
1350 const unsigned int n_q_points,
1352 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1356 const UpdateFlags update_flags = data.update_each;
1379 for (
unsigned int d = 0;
d != dim - 1; ++
d)
1381 Assert(face_no + cell->n_faces() *
d <
1382 data.unit_tangentials.size(),
1385 data.aux[
d].size() <=
1386 data.unit_tangentials[face_no + cell->n_faces() *
d].size(),
1391 data.unit_tangentials[face_no + cell->n_faces() *
d]),
1402 if (dim == spacedim)
1404 for (
unsigned int i = 0; i < n_q_points; ++i)
1414 (face_no == 0 ? -1 : +1);
1418 cross_product_2d(data.aux[0][i]);
1422 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1444 data.contravariant[
point].transpose()[0];
1446 (face_no == 0 ? -1. : +1.) *
1456 cross_product_3d(DX_t[0], DX_t[1]);
1457 cell_normal /= cell_normal.
norm();
1462 cross_product_3d(data.aux[0][
point], cell_normal);
1469 for (
unsigned int i = 0; i < n_q_points; ++i)
1473 data.quadrature_weights[i + data_set];
1478 const double area_ratio =
1480 cell->subface_case(face_no), subface_no);
1489 for (
unsigned int i = 0; i < n_q_points; ++i)
1501 data.covariant[
point].transpose();
1512 template <
int dim,
int spacedim>
1515 const ::MappingFE<dim, spacedim> &mapping,
1516 const typename ::Triangulation<dim, spacedim>::cell_iterator
1518 const unsigned int face_no,
1519 const unsigned int subface_no,
1522 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1526 const unsigned int n_q_points = quadrature.
size();
1528 maybe_compute_q_points<dim, spacedim>(data_set,
1541 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1547 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1553 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1559 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1565 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1587 template <
int dim,
int spacedim>
1591 const unsigned int face_no,
1622 cell->face_orientation(face_no),
1623 cell->face_flip(face_no),
1624 cell->face_rotation(face_no),
1626 quadrature[quadrature.
size() == 1 ? 0 : face_no],
1633 template <
int dim,
int spacedim>
1637 const unsigned int face_no,
1638 const unsigned int subface_no,
1670 cell->face_orientation(face_no),
1671 cell->face_flip(face_no),
1672 cell->face_rotation(face_no),
1674 cell->subface_case(face_no)),
1684 namespace MappingFEImplementation
1688 template <
int dim,
int spacedim,
int rank>
1705 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1706 &mapping_data) !=
nullptr),
1708 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1710 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1713 switch (mapping_kind)
1720 "update_contravariant_transformation"));
1722 for (
unsigned int i = 0; i < input.size(); ++i)
1734 "update_contravariant_transformation"));
1738 "update_volume_elements"));
1743 for (
unsigned int i = 0; i < input.size(); ++i)
1747 output[i] /= data.volume_elements[i];
1759 "update_covariant_transformation"));
1761 for (
unsigned int i = 0; i < input.size(); ++i)
1773 template <
int dim,
int spacedim,
int rank>
1784 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1785 &mapping_data) !=
nullptr),
1787 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1789 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1792 switch (mapping_kind)
1799 "update_covariant_transformation"));
1803 "update_contravariant_transformation"));
1806 for (
unsigned int i = 0; i < output.size(); ++i)
1823 "update_covariant_transformation"));
1826 for (
unsigned int i = 0; i < output.size(); ++i)
1843 "update_covariant_transformation"));
1847 "update_contravariant_transformation"));
1851 "update_volume_elements"));
1854 for (
unsigned int i = 0; i < output.size(); ++i)
1863 output[i] /= data.volume_elements[i];
1876 template <
int dim,
int spacedim>
1887 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
1888 &mapping_data) !=
nullptr),
1890 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1892 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
1895 switch (mapping_kind)
1902 "update_covariant_transformation"));
1906 "update_contravariant_transformation"));
1908 for (
unsigned int q = 0; q < output.size(); ++q)
1909 for (
unsigned int i = 0; i < spacedim; ++i)
1911 double tmp1[dim][dim];
1912 for (
unsigned int J = 0; J < dim; ++J)
1913 for (
unsigned int K = 0;
K < dim; ++
K)
1916 data.contravariant[q][i][0] * input[q][0][J][
K];
1917 for (
unsigned int I = 1; I < dim; ++I)
1919 data.contravariant[q][i][I] * input[q][I][J][
K];
1921 for (
unsigned int j = 0; j < spacedim; ++j)
1924 for (
unsigned int K = 0;
K < dim; ++
K)
1926 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1927 for (
unsigned int J = 1; J < dim; ++J)
1928 tmp2[
K] += data.covariant[q][j][J] * tmp1[J][
K];
1930 for (
unsigned int k = 0; k < spacedim; ++k)
1932 output[q][i][j][k] =
1933 data.covariant[q][k][0] * tmp2[0];
1934 for (
unsigned int K = 1;
K < dim; ++
K)
1935 output[q][i][j][k] +=
1936 data.covariant[q][k][
K] * tmp2[
K];
1948 "update_covariant_transformation"));
1950 for (
unsigned int q = 0; q < output.size(); ++q)
1951 for (
unsigned int i = 0; i < spacedim; ++i)
1953 double tmp1[dim][dim];
1954 for (
unsigned int J = 0; J < dim; ++J)
1955 for (
unsigned int K = 0;
K < dim; ++
K)
1958 data.covariant[q][i][0] * input[q][0][J][
K];
1959 for (
unsigned int I = 1; I < dim; ++I)
1961 data.covariant[q][i][I] * input[q][I][J][
K];
1963 for (
unsigned int j = 0; j < spacedim; ++j)
1966 for (
unsigned int K = 0;
K < dim; ++
K)
1968 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
1969 for (
unsigned int J = 1; J < dim; ++J)
1970 tmp2[
K] += data.covariant[q][j][J] * tmp1[J][
K];
1972 for (
unsigned int k = 0; k < spacedim; ++k)
1974 output[q][i][j][k] =
1975 data.covariant[q][k][0] * tmp2[0];
1976 for (
unsigned int K = 1;
K < dim; ++
K)
1977 output[q][i][j][k] +=
1978 data.covariant[q][k][
K] * tmp2[
K];
1991 "update_covariant_transformation"));
1995 "update_contravariant_transformation"));
1999 "update_volume_elements"));
2001 for (
unsigned int q = 0; q < output.size(); ++q)
2002 for (
unsigned int i = 0; i < spacedim; ++i)
2005 for (
unsigned int I = 0; I < dim; ++I)
2007 data.contravariant[q][i][I] / data.volume_elements[q];
2008 double tmp1[dim][dim];
2009 for (
unsigned int J = 0; J < dim; ++J)
2010 for (
unsigned int K = 0;
K < dim; ++
K)
2012 tmp1[J][
K] = factor[0] * input[q][0][J][
K];
2013 for (
unsigned int I = 1; I < dim; ++I)
2014 tmp1[J][
K] += factor[I] * input[q][I][J][
K];
2016 for (
unsigned int j = 0; j < spacedim; ++j)
2019 for (
unsigned int K = 0;
K < dim; ++
K)
2021 tmp2[
K] = data.covariant[q][j][0] * tmp1[0][
K];
2022 for (
unsigned int J = 1; J < dim; ++J)
2023 tmp2[
K] += data.covariant[q][j][J] * tmp1[J][
K];
2025 for (
unsigned int k = 0; k < spacedim; ++k)
2027 output[q][i][j][k] =
2028 data.covariant[q][k][0] * tmp2[0];
2029 for (
unsigned int K = 1;
K < dim; ++
K)
2030 output[q][i][j][k] +=
2031 data.covariant[q][k][
K] * tmp2[
K];
2046 template <
int dim,
int spacedim,
int rank>
2057 const typename ::MappingFE<dim, spacedim>::InternalData *
>(
2058 &mapping_data) !=
nullptr),
2060 const typename ::MappingFE<dim, spacedim>::InternalData &data =
2062 const typename ::MappingFE<dim, spacedim>::InternalData &
>(
2065 switch (mapping_kind)
2072 "update_covariant_transformation"));
2074 for (
unsigned int i = 0; i < output.size(); ++i)
2089 template <
int dim,
int spacedim>
2105 template <
int dim,
int spacedim>
2121 template <
int dim,
int spacedim>
2129 switch (mapping_kind)
2153 template <
int dim,
int spacedim>
2166 switch (mapping_kind)
2172 "update_covariant_transformation"));
2174 for (
unsigned int q = 0; q < output.size(); ++q)
2175 for (
unsigned int i = 0; i < spacedim; ++i)
2176 for (
unsigned int j = 0; j < spacedim; ++j)
2179 for (
unsigned int K = 0;
K < dim; ++
K)
2181 tmp[
K] = data.
covariant[q][j][0] * input[q][i][0][
K];
2182 for (
unsigned int J = 1; J < dim; ++J)
2183 tmp[
K] += data.
covariant[q][j][J] * input[q][i][J][
K];
2185 for (
unsigned int k = 0; k < spacedim; ++k)
2187 output[q][i][j][k] = data.
covariant[q][k][0] * tmp[0];
2188 for (
unsigned int K = 1;
K < dim; ++
K)
2189 output[q][i][j][k] += data.
covariant[q][k][
K] * tmp[
K];
2202 template <
int dim,
int spacedim>
2210 switch (mapping_kind)
2229 template <
int spacedim>
2231 check_all_manifold_ids_identical(
2239 template <
int spacedim>
2241 check_all_manifold_ids_identical(
2244 const auto m_id = cell->manifold_id();
2246 for (
const auto f : cell->face_indices())
2247 if (m_id != cell->face(f)->manifold_id())
2255 template <
int spacedim>
2257 check_all_manifold_ids_identical(
2260 const auto m_id = cell->manifold_id();
2262 for (
const auto f : cell->face_indices())
2263 if (m_id != cell->face(f)->manifold_id())
2266 for (
const auto l : cell->line_indices())
2267 if (m_id != cell->line(
l)->manifold_id())
2276 template <
int dim,
int spacedim>
2277 std::vector<Point<spacedim>>
2282 check_all_manifold_ids_identical(cell),
2284 "All entities of a cell need to have the same manifold id as the cell has."));
2288 for (
const unsigned int i : cell->vertex_indices())
2291 std::vector<Point<spacedim>> mapping_support_points(
2292 fe->get_unit_support_points().size());
2295 mapping_support_point_weights,
2296 mapping_support_points);
2298 return mapping_support_points;
2303 template <
int dim,
int spacedim>
2313 template <
int dim,
int spacedim>
2319 ExcMessage(
"The dimension of your mapping (" +
2321 ") and the reference cell cell_type (" +
2323 " ) do not agree."));
2331 #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)
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
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)
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
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 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 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 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< Point< dim > > & get_points() const
const std::vector< double > & get_weights() const
unsigned int size() const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
unsigned int size() const
unsigned int max_n_quadrature_points() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
@ 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
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
@ mapping_covariant_gradient
@ mapping_contravariant_hessian
@ mapping_covariant_hessian
@ mapping_contravariant_gradient
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
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_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_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 maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 5, spacedim >> &jacobian_pushed_forward_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 4, spacedim >> &jacobian_pushed_forward_2nd_derivatives)
void maybe_compute_q_points(const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Point< spacedim >> &quadrature_points)
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data)
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, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
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)
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 transform_fields(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 maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
static const unsigned int invalid_unsigned_int
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)