deal.II version GIT relicensing-3321-g4b434670a0 2025-05-16 17:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mapping_fe.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
22#include <deal.II/base/table.h>
24
25#include <deal.II/fe/fe_poly.h>
28
31#include <deal.II/grid/tria.h>
33
35
36#include <boost/container/small_vector.hpp>
37
38#include <algorithm>
39#include <array>
40#include <cmath>
41#include <memory>
42#include <numeric>
43
44
46
47
48template <int dim, int spacedim>
51 : fe(fe)
52 , polynomial_degree(fe.tensor_degree())
53 , n_shape_functions(fe.n_dofs_per_cell())
54{}
55
56
57
58template <int dim, int spacedim>
59std::size_t
76
77
78
79template <int dim, int spacedim>
80void
82 const Quadrature<dim> &q)
83{
84 // store the flags in the internal data object so we can access them
85 // in fill_fe_*_values()
86 this->update_each = update_flags;
87
88 const unsigned int n_q_points = q.size();
89
90 if (this->update_each & update_covariant_transformation)
91 covariant.resize(n_q_points);
92
93 if (this->update_each & update_contravariant_transformation)
94 contravariant.resize(n_q_points);
95
96 if (this->update_each & update_volume_elements)
97 volume_elements.resize(n_q_points);
98
99 // see if we need the (transformation) shape function values
100 // and/or gradients and resize the necessary arrays
101 if (this->update_each & update_quadrature_points)
102 shape_values.resize(n_shape_functions * n_q_points);
103
104 if (this->update_each &
112 shape_derivatives.resize(n_shape_functions * n_q_points);
113
114 if (this->update_each &
116 shape_second_derivatives.resize(n_shape_functions * n_q_points);
117
118 if (this->update_each & (update_jacobian_2nd_derivatives |
120 shape_third_derivatives.resize(n_shape_functions * n_q_points);
121
122 if (this->update_each & (update_jacobian_3rd_derivatives |
124 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
125
126 // now also fill the various fields with their correct values
127 compute_shape_function_values(q.get_points());
128
129 // copy (projected) quadrature weights
130 quadrature_weights = q.get_weights();
131}
132
133
134
135template <int dim, int spacedim>
136void
138 const UpdateFlags update_flags,
139 const Quadrature<dim> &q,
140 const unsigned int n_original_q_points)
141{
142 reinit(update_flags, q);
143
144 if (this->update_each &
147 {
148 aux.resize(dim - 1,
149 std::vector<Tensor<1, spacedim>>(n_original_q_points));
150
151 // Compute tangentials to the unit cell.
152 const auto reference_cell = this->fe.reference_cell();
153 const auto n_faces = reference_cell.n_faces();
154
155 for (unsigned int i = 0; i < n_faces; ++i)
156 {
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 face_tangent_vector<dim>(i, 0));
161 if (dim > 2)
162 {
163 unit_tangentials[n_faces + i].resize(n_original_q_points);
164 std::fill(unit_tangentials[n_faces + i].begin(),
165 unit_tangentials[n_faces + i].end(),
166 reference_cell.template face_tangent_vector<dim>(i, 1));
167 }
168 }
169 }
170}
171
172
173
174template <int dim, int spacedim>
175void
177 const std::vector<Point<dim>> &unit_points)
178{
179 const auto fe_poly = dynamic_cast<const FE_Poly<dim, spacedim> *>(&this->fe);
180
181 Assert(fe_poly != nullptr, ExcNotImplemented());
182
183 const auto &tensor_pols = fe_poly->get_poly_space();
184
185 const unsigned int n_shape_functions = fe.n_dofs_per_cell();
186 const unsigned int n_points = unit_points.size();
187
188 std::vector<double> values;
189 std::vector<Tensor<1, dim>> grads;
190 if (shape_values.size() != 0)
191 {
192 Assert(shape_values.size() == n_shape_functions * n_points,
194 values.resize(n_shape_functions);
195 }
196 if (shape_derivatives.size() != 0)
197 {
198 Assert(shape_derivatives.size() == n_shape_functions * n_points,
200 grads.resize(n_shape_functions);
201 }
202
203 std::vector<Tensor<2, dim>> grad2;
204 if (shape_second_derivatives.size() != 0)
205 {
206 Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
208 grad2.resize(n_shape_functions);
209 }
210
211 std::vector<Tensor<3, dim>> grad3;
212 if (shape_third_derivatives.size() != 0)
213 {
214 Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
216 grad3.resize(n_shape_functions);
217 }
218
219 std::vector<Tensor<4, dim>> grad4;
220 if (shape_fourth_derivatives.size() != 0)
221 {
222 Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
224 grad4.resize(n_shape_functions);
225 }
226
227
228 if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
229 shape_second_derivatives.size() != 0 ||
230 shape_third_derivatives.size() != 0 ||
231 shape_fourth_derivatives.size() != 0)
232 for (unsigned int point = 0; point < n_points; ++point)
233 {
234 tensor_pols.evaluate(
235 unit_points[point], values, grads, grad2, grad3, grad4);
236
237 if (shape_values.size() != 0)
238 for (unsigned int i = 0; i < n_shape_functions; ++i)
239 shape(point, i) = values[i];
240
241 if (shape_derivatives.size() != 0)
242 for (unsigned int i = 0; i < n_shape_functions; ++i)
243 derivative(point, i) = grads[i];
244
245 if (shape_second_derivatives.size() != 0)
246 for (unsigned int i = 0; i < n_shape_functions; ++i)
247 second_derivative(point, i) = grad2[i];
248
249 if (shape_third_derivatives.size() != 0)
250 for (unsigned int i = 0; i < n_shape_functions; ++i)
251 third_derivative(point, i) = grad3[i];
252
253 if (shape_fourth_derivatives.size() != 0)
254 for (unsigned int i = 0; i < n_shape_functions; ++i)
255 fourth_derivative(point, i) = grad4[i];
256 }
257}
258
259
260namespace internal
261{
262 namespace MappingFEImplementation
263 {
264 namespace
265 {
272 template <int dim, int spacedim>
273 void
274 maybe_compute_q_points(
275 const typename QProjector<dim>::DataSetDescriptor data_set,
276 const typename ::MappingFE<dim, spacedim>::InternalData &data,
277 std::vector<Point<spacedim>> &quadrature_points,
278 const unsigned int n_q_points)
279 {
280 const UpdateFlags update_flags = data.update_each;
281
282 if (update_flags & update_quadrature_points)
283 for (unsigned int point = 0; point < n_q_points; ++point)
284 {
285 const double *shape = &data.shape(point + data_set, 0);
286 Point<spacedim> result =
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;
292 }
293 }
294
295
296
305 template <int dim, int spacedim>
306 void
307 maybe_update_Jacobians(
308 const CellSimilarity::Similarity cell_similarity,
309 const typename ::QProjector<dim>::DataSetDescriptor data_set,
310 const typename ::MappingFE<dim, spacedim>::InternalData &data,
311 const unsigned int n_q_points)
312 {
313 const UpdateFlags update_flags = data.update_each;
314
315 if (update_flags & update_contravariant_transformation)
316 // if the current cell is just a
317 // translation of the previous one, no
318 // need to recompute jacobians...
319 if (cell_similarity != CellSimilarity::translation)
320 {
321 std::fill(data.contravariant.begin(),
322 data.contravariant.end(),
324
325 Assert(data.n_shape_functions > 0, ExcInternalError());
326
327 for (unsigned int point = 0; point < n_q_points; ++point)
328 {
329 double result[spacedim][dim];
330
331 // peel away part of sum to avoid zeroing the
332 // entries and adding for the first time
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)
340 result[i][j] +=
341 data.derivative(point + data_set, k)[j] *
342 data.mapping_support_points[k][i];
343
344 // write result into contravariant data. for
345 // j=dim in the case dim<spacedim, there will
346 // never be any nonzero data that arrives in
347 // here, so it is ok anyway because it was
348 // initialized to zero at the initialization
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];
352 }
353 }
354
355 if (update_flags & update_covariant_transformation)
356 if (cell_similarity != CellSimilarity::translation)
357 {
358 for (unsigned int point = 0; point < n_q_points; ++point)
359 {
360 data.covariant[point] =
361 (data.contravariant[point]).covariant_form();
362 }
363 }
364
365 if (update_flags & update_volume_elements)
366 if (cell_similarity != CellSimilarity::translation)
367 {
368 for (unsigned int point = 0; point < n_q_points; ++point)
369 data.volume_elements[point] =
370 data.contravariant[point].determinant();
371 }
372 }
373
380 template <int dim, int spacedim>
381 void
382 maybe_update_jacobian_grads(
383 const CellSimilarity::Similarity cell_similarity,
384 const typename QProjector<dim>::DataSetDescriptor data_set,
385 const typename ::MappingFE<dim, spacedim>::InternalData &data,
386 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads,
387 const unsigned int n_q_points)
388 {
389 const UpdateFlags update_flags = data.update_each;
390 if (update_flags & update_jacobian_grads)
391 {
392 AssertIndexRange(n_q_points, jacobian_grads.size() + 1);
393
394 if (cell_similarity != CellSimilarity::translation)
395 for (unsigned int point = 0; point < n_q_points; ++point)
396 {
397 const Tensor<2, dim> *second =
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)
403 result[i][j][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)
409 result[i][j][l] +=
410 (second[k][j][l] *
411 data.mapping_support_points[k][i]);
412
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];
417 }
418 }
419 }
420
427 template <int dim, int spacedim>
428 void
429 maybe_update_jacobian_pushed_forward_grads(
430 const CellSimilarity::Similarity cell_similarity,
431 const typename QProjector<dim>::DataSetDescriptor data_set,
432 const typename ::MappingFE<dim, spacedim>::InternalData &data,
433 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads,
434 const unsigned int n_q_points)
435 {
436 const UpdateFlags update_flags = data.update_each;
437 if (update_flags & update_jacobian_pushed_forward_grads)
438 {
439 AssertIndexRange(n_q_points,
440 jacobian_pushed_forward_grads.size() + 1);
441
442 if (cell_similarity != CellSimilarity::translation)
443 {
444 double tmp[spacedim][spacedim][spacedim];
445 for (unsigned int point = 0; point < n_q_points; ++point)
446 {
447 const Tensor<2, dim> *second =
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)
459 result[i][j][l] +=
460 (second[k][j][l] *
461 data.mapping_support_points[k][i]);
462
463 // first push forward the j-components
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)
467 {
468 tmp[i][j][l] =
469 result[i][0][l] * data.covariant[point][j][0];
470 for (unsigned int jr = 1; jr < dim; ++jr)
471 {
472 tmp[i][j][l] += result[i][jr][l] *
473 data.covariant[point][j][jr];
474 }
475 }
476
477 // now, pushing forward the l-components
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)
481 {
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)
485 {
486 jacobian_pushed_forward_grads[point][i][j][l] +=
487 tmp[i][j][lr] * data.covariant[point][l][lr];
488 }
489 }
490 }
491 }
492 }
493 }
494
501 template <int dim, int spacedim>
502 void
503 maybe_update_jacobian_2nd_derivatives(
504 const CellSimilarity::Similarity cell_similarity,
505 const typename QProjector<dim>::DataSetDescriptor data_set,
506 const typename ::MappingFE<dim, spacedim>::InternalData &data,
507 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives,
508 const unsigned int n_q_points)
509 {
510 const UpdateFlags update_flags = data.update_each;
511 if (update_flags & update_jacobian_2nd_derivatives)
512 {
513 AssertIndexRange(n_q_points, jacobian_2nd_derivatives.size() + 1);
514
515 if (cell_similarity != CellSimilarity::translation)
516 {
517 for (unsigned int point = 0; point < n_q_points; ++point)
518 {
519 const Tensor<3, dim> *third =
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)
526 result[i][j][l][m] =
527 (third[0][j][l][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] +=
535 (third[k][j][l][m] *
536 data.mapping_support_points[k][i]);
537
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] =
543 result[i][j][l][m];
544 }
545 }
546 }
547 }
548
556 template <int dim, int spacedim>
557 void
558 maybe_update_jacobian_pushed_forward_2nd_derivatives(
559 const CellSimilarity::Similarity cell_similarity,
560 const typename QProjector<dim>::DataSetDescriptor data_set,
561 const typename ::MappingFE<dim, spacedim>::InternalData &data,
562 std::vector<Tensor<4, spacedim>>
563 &jacobian_pushed_forward_2nd_derivatives,
564 const unsigned int n_q_points)
565 {
566 const UpdateFlags update_flags = data.update_each;
568 {
569 AssertIndexRange(n_q_points,
570 jacobian_pushed_forward_2nd_derivatives.size() +
571 1);
572
573 if (cell_similarity != CellSimilarity::translation)
574 {
575 double tmp[spacedim][spacedim][spacedim][spacedim];
576 for (unsigned int point = 0; point < n_q_points; ++point)
577 {
578 const Tensor<3, dim> *third =
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)
585 result[i][j][l][m] =
586 (third[0][j][l][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] +=
594 (third[k][j][l][m] *
595 data.mapping_support_points[k][i]);
596
597 // push forward the j-coordinate
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)
602 {
603 jacobian_pushed_forward_2nd_derivatives
604 [point][i][j][l][m] =
605 result[i][0][l][m] *
606 data.covariant[point][j][0];
607 for (unsigned int jr = 1; jr < dim; ++jr)
608 jacobian_pushed_forward_2nd_derivatives[point]
609 [i][j][l]
610 [m] +=
611 result[i][jr][l][m] *
612 data.covariant[point][j][jr];
613 }
614
615 // push forward the l-coordinate
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)
620 {
621 tmp[i][j][l][m] =
622 jacobian_pushed_forward_2nd_derivatives[point]
623 [i][j][0]
624 [m] *
625 data.covariant[point][l][0];
626 for (unsigned int lr = 1; lr < dim; ++lr)
627 tmp[i][j][l][m] +=
628 jacobian_pushed_forward_2nd_derivatives
629 [point][i][j][lr][m] *
630 data.covariant[point][l][lr];
631 }
632
633 // push forward the m-coordinate
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)
638 {
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]
644 [i][j][l]
645 [m] +=
646 tmp[i][j][l][mr] *
647 data.covariant[point][m][mr];
648 }
649 }
650 }
651 }
652 }
653
660 template <int dim, int spacedim>
661 void
662 maybe_update_jacobian_3rd_derivatives(
663 const CellSimilarity::Similarity cell_similarity,
664 const typename QProjector<dim>::DataSetDescriptor data_set,
665 const typename ::MappingFE<dim, spacedim>::InternalData &data,
666 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives,
667 const unsigned int n_q_points)
668 {
669 const UpdateFlags update_flags = data.update_each;
670 if (update_flags & update_jacobian_3rd_derivatives)
671 {
672 AssertIndexRange(n_q_points, jacobian_3rd_derivatives.size() + 1);
673
674 if (cell_similarity != CellSimilarity::translation)
675 {
676 for (unsigned int point = 0; point < n_q_points; ++point)
677 {
678 const Tensor<4, dim> *fourth =
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]);
698
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];
706 }
707 }
708 }
709 }
710
718 template <int dim, int spacedim>
719 void
720 maybe_update_jacobian_pushed_forward_3rd_derivatives(
721 const CellSimilarity::Similarity cell_similarity,
722 const typename QProjector<dim>::DataSetDescriptor data_set,
723 const typename ::MappingFE<dim, spacedim>::InternalData &data,
724 std::vector<Tensor<5, spacedim>>
725 &jacobian_pushed_forward_3rd_derivatives,
726 const unsigned int n_q_points)
727 {
728 const UpdateFlags update_flags = data.update_each;
730 {
731 AssertIndexRange(n_q_points,
732 jacobian_pushed_forward_3rd_derivatives.size() +
733 1);
734
735 if (cell_similarity != CellSimilarity::translation)
736 {
737 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
738 for (unsigned int point = 0; point < n_q_points; ++point)
739 {
740 const Tensor<4, dim> *fourth =
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]);
760
761 // push-forward the j-coordinate
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)
767 {
768 tmp[i][j][l][m][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];
775 }
776
777 // push-forward the l-coordinate
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)
783 {
784 jacobian_pushed_forward_3rd_derivatives
785 [point][i][j][l][m][n] =
786 tmp[i][j][0][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];
793 }
794
795 // push-forward the m-coordinate
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)
801 {
802 tmp[i][j][l][m][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];
811 }
812
813 // push-forward the n-coordinate
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)
819 {
820 jacobian_pushed_forward_3rd_derivatives
821 [point][i][j][l][m][n] =
822 tmp[i][j][l][m][0] *
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];
829 }
830 }
831 }
832 }
833 }
834 } // namespace
835 } // namespace MappingFEImplementation
836} // namespace internal
837
838
839
840template <int dim, int spacedim>
842 : fe(fe.clone())
843 , polynomial_degree(fe.tensor_degree())
844{
846 ExcMessage("It only makes sense to create polynomial mappings "
847 "with a polynomial degree greater or equal to one."));
848 Assert(fe.n_components() == 1, ExcNotImplemented());
849
850 Assert(fe.has_support_points(), ExcNotImplemented());
851
852 const auto &mapping_support_points = fe.get_unit_support_points();
853
854 const auto reference_cell = fe.reference_cell();
855
856 const unsigned int n_points = mapping_support_points.size();
857 const unsigned int n_shape_functions = reference_cell.n_vertices();
858
860 Table<2, double>(n_points, n_shape_functions);
861
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],
866 i);
867}
868
869
870
871template <int dim, int spacedim>
873 : fe(mapping.fe->clone())
874 , polynomial_degree(mapping.polynomial_degree)
875 , mapping_support_point_weights(mapping.mapping_support_point_weights)
876{}
877
878
879
880template <int dim, int spacedim>
881std::unique_ptr<Mapping<dim, spacedim>>
883{
884 return std::make_unique<MappingFE<dim, spacedim>>(*this);
885}
886
887
888
889template <int dim, int spacedim>
890unsigned int
892{
893 return polynomial_degree;
894}
895
896
897
898template <int dim, int spacedim>
902 const Point<dim> &p) const
903{
904 const std::vector<Point<spacedim>> support_points =
905 this->compute_mapping_support_points(cell);
906
907 Point<spacedim> mapped_point;
908
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);
911
912 return mapped_point;
913}
914
915
916
917template <int dim, int spacedim>
921 const Point<spacedim> &p) const
922{
923 const std::vector<Point<spacedim>> support_points =
924 this->compute_mapping_support_points(cell);
925
926 const double eps = 1.e-12 * cell->diameter();
927 const unsigned int loop_limit = 10;
928
929 Point<dim> p_unit;
930
931 unsigned int loop = 0;
932
933 // This loop solves the following problem:
934 // grad_F^T residual + (grad_F^T grad_F + grad_F^T hess_F^T dp) dp = 0
935 // where the term
936 // (grad_F^T hess_F dp) is approximated by (-hess_F * residual)
937 // This is basically a second order approximation of Newton method, where the
938 // Jacobian is corrected with a higher order term coming from the hessian.
939 do
940 {
941 Point<spacedim> mapped_point;
942
943 // Transpose of the gradient map
946
947 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
948 {
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)
953 {
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];
957 }
958 }
959
960 // Residual
961 const auto residual = p - mapped_point;
962 // Project the residual on the reference coordinate system
963 // to compute the error, and to filter components orthogonal to the
964 // manifold, and compute a 2nd order correction of the metric tensor
965 const auto grad_FT_residual = apply_transformation(grad_FT, residual);
966
967 // Do not invert nor compute the metric if not necessary.
968 if (grad_FT_residual.norm() <= eps)
969 break;
970
971 // Now compute the (corrected) metric tensor
972 Tensor<2, dim> corrected_metric_tensor;
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;
977
978 // And compute the update
979 const auto g_inverse = invert(corrected_metric_tensor);
980 p_unit -= Point<dim>(g_inverse * grad_FT_residual);
981
982 ++loop;
983 }
984 while (loop < loop_limit);
985
986 // Here we check that in the last execution of while the first
987 // condition was already wrong, meaning the residual was below
988 // eps. Only if the first condition failed, loop will have been
989 // increased and tested, and thus have reached the limit.
990 AssertThrow(loop < loop_limit,
992
993 return p_unit;
994}
995
996
997
998template <int dim, int spacedim>
1001{
1002 // add flags if the respective quantities are necessary to compute
1003 // what we need. note that some flags appear in both the conditions
1004 // and in subsequent set operations. this leads to some circular
1005 // logic. the only way to treat this is to iterate. since there are
1006 // 5 if-clauses in the loop, it will take at most 5 iterations to
1007 // converge. do them:
1008 UpdateFlags out = in;
1009 for (unsigned int i = 0; i < 5; ++i)
1010 {
1011 // The following is a little incorrect:
1012 // If not applied on a face,
1013 // update_boundary_forms does not
1014 // make sense. On the other hand,
1015 // it is necessary on a
1016 // face. Currently,
1017 // update_boundary_forms is simply
1018 // ignored for the interior of a
1019 // cell.
1021 out |= update_boundary_forms;
1022
1027
1028 if (out &
1033
1034 // The contravariant transformation is used in the Piola
1035 // transformation, which requires the determinant of the Jacobi
1036 // matrix of the transformation. Because we have no way of
1037 // knowing here whether the finite element wants to use the
1038 // contravariant or the Piola transforms, we add the JxW values
1039 // to the list of flags to be updated for each cell.
1042
1043 // the same is true when computing normal vectors: they require
1044 // the determinant of the Jacobian
1045 if (out & update_normal_vectors)
1047 }
1048
1049 return out;
1050}
1051
1052
1053
1054template <int dim, int spacedim>
1055std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1057 const Quadrature<dim> &q) const
1058{
1059 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1060 std::make_unique<InternalData>(*this->fe);
1061 data_ptr->reinit(this->requires_update_flags(update_flags), q);
1062 return data_ptr;
1063}
1064
1065
1066
1067template <int dim, int spacedim>
1068std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1070 const UpdateFlags update_flags,
1071 const hp::QCollection<dim - 1> &quadrature) const
1072{
1073 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1074 std::make_unique<InternalData>(*this->fe);
1075 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1076 data.initialize_face(this->requires_update_flags(update_flags),
1078 this->fe->reference_cell(), quadrature),
1079 quadrature.max_n_quadrature_points());
1080
1081 return data_ptr;
1082}
1083
1084
1085
1086template <int dim, int spacedim>
1087std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1089 const UpdateFlags update_flags,
1090 const Quadrature<dim - 1> &quadrature) const
1091{
1092 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1093 std::make_unique<InternalData>(*this->fe);
1094 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1095 data.initialize_face(this->requires_update_flags(update_flags),
1097 this->fe->reference_cell(), quadrature),
1098 quadrature.size());
1099
1100 return data_ptr;
1101}
1102
1103
1104
1105template <int dim, int spacedim>
1109 const CellSimilarity::Similarity cell_similarity,
1110 const Quadrature<dim> &quadrature,
1111 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1113 &output_data) const
1114{
1115 // ensure that the following static_cast is really correct:
1116 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1118 const InternalData &data = static_cast<const InternalData &>(internal_data);
1119
1120 const unsigned int n_q_points = quadrature.size();
1121
1122 // recompute the support points of the transformation of this
1123 // cell. we tried to be clever here in an earlier version of the
1124 // library by checking whether the cell is the same as the one we
1125 // had visited last, but it turns out to be difficult to determine
1126 // that because a cell for the purposes of a mapping is
1127 // characterized not just by its (triangulation, level, index)
1128 // triple, but also by the locations of its vertices, the manifold
1129 // object attached to the cell and all of its bounding faces/edges,
1130 // etc. to reliably test that the "cell" we are on is, therefore,
1131 // not easily done
1132 data.mapping_support_points = this->compute_mapping_support_points(cell);
1133 data.cell_of_current_support_points = cell;
1134
1135 // if the order of the mapping is greater than 1, then do not reuse any cell
1136 // similarity information. This is necessary because the cell similarity
1137 // value is computed with just cell vertices and does not take into account
1138 // cell curvature.
1139 const CellSimilarity::Similarity computed_cell_similarity =
1140 (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
1141
1142 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1144 data,
1145 output_data.quadrature_points,
1146 n_q_points);
1147
1148 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1149 computed_cell_similarity,
1151 data,
1152 n_q_points);
1153
1154 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1155 computed_cell_similarity,
1157 data,
1158 output_data.jacobian_grads,
1159 n_q_points);
1160
1161 internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
1162 dim,
1163 spacedim>(computed_cell_similarity,
1165 data,
1167 n_q_points);
1168
1169 internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
1170 dim,
1171 spacedim>(computed_cell_similarity,
1173 data,
1174 output_data.jacobian_2nd_derivatives,
1175 n_q_points);
1176
1177 internal::MappingFEImplementation::
1178 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1179 computed_cell_similarity,
1181 data,
1183 n_q_points);
1184
1185 internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
1186 dim,
1187 spacedim>(computed_cell_similarity,
1189 data,
1190 output_data.jacobian_3rd_derivatives,
1191 n_q_points);
1192
1193 internal::MappingFEImplementation::
1194 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1195 computed_cell_similarity,
1197 data,
1199 n_q_points);
1200
1201 const UpdateFlags update_flags = data.update_each;
1202 const std::vector<double> &weights = quadrature.get_weights();
1203
1204 // Multiply quadrature weights by absolute value of Jacobian determinants or
1205 // the area element g=sqrt(DX^t DX) in case of codim > 0
1206
1207 if (update_flags & (update_normal_vectors | update_JxW_values))
1208 {
1209 AssertDimension(output_data.JxW_values.size(), n_q_points);
1210
1211 Assert(!(update_flags & update_normal_vectors) ||
1212 (output_data.normal_vectors.size() == n_q_points),
1213 ExcDimensionMismatch(output_data.normal_vectors.size(),
1214 n_q_points));
1215
1216
1217 if (computed_cell_similarity != CellSimilarity::translation)
1218 for (unsigned int point = 0; point < n_q_points; ++point)
1219 {
1220 if (dim == spacedim)
1221 {
1222 const double det = data.contravariant[point].determinant();
1223
1224 // check for distorted cells.
1225
1226 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1227 // 1e12 in 2d. might want to find a finer
1228 // (dimension-independent) criterion
1229 Assert(det >
1230 1e-12 * Utilities::fixed_power<dim>(
1231 cell->diameter() / std::sqrt(double(dim))),
1233 cell->center(), det, point)));
1234
1235 output_data.JxW_values[point] = weights[point] * det;
1236 }
1237 // if dim==spacedim, then there is no cell normal to
1238 // compute. since this is for FEValues (and not FEFaceValues),
1239 // there are also no face normals to compute
1240 else // codim>0 case
1241 {
1242 Tensor<1, spacedim> DX_t[dim];
1243 for (unsigned int i = 0; i < spacedim; ++i)
1244 for (unsigned int j = 0; j < dim; ++j)
1245 DX_t[j][i] = data.contravariant[point][i][j];
1246
1247 Tensor<2, dim> G; // First fundamental form
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];
1251
1252 output_data.JxW_values[point] =
1253 std::sqrt(determinant(G)) * weights[point];
1254
1255 if (computed_cell_similarity ==
1257 {
1258 // we only need to flip the normal
1259 if (update_flags & update_normal_vectors)
1260 output_data.normal_vectors[point] *= -1.;
1261 }
1262 else
1263 {
1264 if (update_flags & update_normal_vectors)
1265 {
1266 Assert(spacedim == dim + 1,
1267 ExcMessage(
1268 "There is no (unique) cell normal for " +
1270 "-dimensional cells in " +
1271 Utilities::int_to_string(spacedim) +
1272 "-dimensional space. This only works if the "
1273 "space dimension is one greater than the "
1274 "dimensionality of the mesh cells."));
1275
1276 if (dim == 1)
1277 output_data.normal_vectors[point] =
1278 cross_product_2d(-DX_t[0]);
1279 else // dim == 2
1280 output_data.normal_vectors[point] =
1281 cross_product_3d(DX_t[0], DX_t[1]);
1282
1283 output_data.normal_vectors[point] /=
1284 output_data.normal_vectors[point].norm();
1285
1286 if (cell->direction_flag() == false)
1287 output_data.normal_vectors[point] *= -1.;
1288 }
1289 }
1290 } // codim>0 case
1291 }
1292 }
1293
1294
1295
1296 // copy values from InternalData to vector given by reference
1297 if (update_flags & update_jacobians)
1298 {
1299 AssertDimension(output_data.jacobians.size(), n_q_points);
1300 if (computed_cell_similarity != CellSimilarity::translation)
1301 for (unsigned int point = 0; point < n_q_points; ++point)
1302 output_data.jacobians[point] = data.contravariant[point];
1303 }
1304
1305 // copy values from InternalData to vector given by reference
1306 if (update_flags & update_inverse_jacobians)
1307 {
1308 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1309 if (computed_cell_similarity != CellSimilarity::translation)
1310 for (unsigned int point = 0; point < n_q_points; ++point)
1311 output_data.inverse_jacobians[point] =
1312 data.covariant[point].transpose();
1313 }
1314
1315 return computed_cell_similarity;
1316}
1317
1318
1319
1320namespace internal
1321{
1322 namespace MappingFEImplementation
1323 {
1324 namespace
1325 {
1336 template <int dim, int spacedim>
1337 void
1338 maybe_compute_face_data(
1339 const ::MappingFE<dim, spacedim> &mapping,
1340 const typename ::Triangulation<dim, spacedim>::cell_iterator
1341 &cell,
1342 const unsigned int face_no,
1343 const unsigned int subface_no,
1344 const unsigned int n_q_points,
1345 const typename QProjector<dim>::DataSetDescriptor data_set,
1346 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1348 &output_data)
1349 {
1350 const UpdateFlags update_flags = data.update_each;
1351
1352 if (update_flags &
1355 {
1356 if (update_flags & update_boundary_forms)
1357 AssertIndexRange(n_q_points,
1358 output_data.boundary_forms.size() + 1);
1359 if (update_flags & update_normal_vectors)
1360 AssertIndexRange(n_q_points,
1361 output_data.normal_vectors.size() + 1);
1362 if (update_flags & update_JxW_values)
1363 AssertIndexRange(n_q_points, output_data.JxW_values.size() + 1);
1364
1365 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1366
1367 // first compute some common data that is used for evaluating
1368 // all of the flags below
1369
1370 // map the unit tangentials to the real cell. checking for
1371 // d!=dim-1 eliminates compiler warnings regarding unsigned int
1372 // expressions < 0.
1373 for (unsigned int d = 0; d != dim - 1; ++d)
1374 {
1375 Assert(face_no + cell->n_faces() * d <
1376 data.unit_tangentials.size(),
1378 Assert(
1379 data.aux[d].size() <=
1380 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1382
1383 mapping.transform(
1385 data.unit_tangentials[face_no + cell->n_faces() * d]),
1387 data,
1388 make_array_view(data.aux[d]));
1389 }
1390
1391 if (update_flags & update_boundary_forms)
1392 {
1393 // if dim==spacedim, we can use the unit tangentials to
1394 // compute the boundary form by simply taking the cross
1395 // product
1396 if (dim == spacedim)
1397 {
1398 for (unsigned int i = 0; i < n_q_points; ++i)
1399 switch (dim)
1400 {
1401 case 1:
1402 // in 1d, we don't have access to any of the
1403 // data.aux fields (because it has only dim-1
1404 // components), but we can still compute the
1405 // boundary form by simply looking at the number
1406 // of the face
1407 output_data.boundary_forms[i][0] =
1408 (face_no == 0 ? -1 : +1);
1409 break;
1410 case 2:
1411 output_data.boundary_forms[i] =
1412 cross_product_2d(data.aux[0][i]);
1413 break;
1414 case 3:
1415 output_data.boundary_forms[i] =
1416 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1417 break;
1418 default:
1420 }
1421 }
1422 else //(dim < spacedim)
1423 {
1424 // in the codim-one case, the boundary form results from
1425 // the cross product of all the face tangential vectors
1426 // and the cell normal vector
1427 //
1428 // to compute the cell normal, use the same method used in
1429 // fill_fe_values for cells above
1430 AssertIndexRange(n_q_points, data.contravariant.size() + 1);
1431
1432 for (unsigned int point = 0; point < n_q_points; ++point)
1433 {
1434 if (dim == 1)
1435 {
1436 // J is a tangent vector
1437 output_data.boundary_forms[point] =
1438 data.contravariant[point].transpose()[0];
1439 output_data.boundary_forms[point] /=
1440 (face_no == 0 ? -1. : +1.) *
1441 output_data.boundary_forms[point].norm();
1442 }
1443
1444 if (dim == 2)
1445 {
1447 data.contravariant[point].transpose();
1448
1449 Tensor<1, spacedim> cell_normal =
1450 cross_product_3d(DX_t[0], DX_t[1]);
1451 cell_normal /= cell_normal.norm();
1452
1453 // then compute the face normal from the face
1454 // tangent and the cell normal:
1455 output_data.boundary_forms[point] =
1456 cross_product_3d(data.aux[0][point], cell_normal);
1457 }
1458 }
1459 }
1460 }
1461
1462 if (update_flags & update_JxW_values)
1463 for (unsigned int i = 0; i < n_q_points; ++i)
1464 {
1465 output_data.JxW_values[i] =
1466 output_data.boundary_forms[i].norm() *
1467 data.quadrature_weights[i + data_set];
1468
1469 if (subface_no != numbers::invalid_unsigned_int)
1470 {
1471 if (dim == 2)
1472 {
1473 const double area_ratio =
1474 1. / cell->reference_cell()
1475 .face_reference_cell(face_no)
1476 .n_isotropic_children();
1477 output_data.JxW_values[i] *= area_ratio;
1478 }
1479 else
1481 }
1482 }
1483
1484 if (update_flags & update_normal_vectors)
1485 for (unsigned int i = 0; i < n_q_points; ++i)
1486 output_data.normal_vectors[i] =
1487 Point<spacedim>(output_data.boundary_forms[i] /
1488 output_data.boundary_forms[i].norm());
1489
1490 if (update_flags & update_jacobians)
1491 for (unsigned int point = 0; point < n_q_points; ++point)
1492 output_data.jacobians[point] = data.contravariant[point];
1493
1494 if (update_flags & update_inverse_jacobians)
1495 for (unsigned int point = 0; point < n_q_points; ++point)
1496 output_data.inverse_jacobians[point] =
1497 data.covariant[point].transpose();
1498 }
1499 }
1500
1501
1508 template <int dim, int spacedim>
1509 void
1511 const ::MappingFE<dim, spacedim> &mapping,
1512 const typename ::Triangulation<dim, spacedim>::cell_iterator
1513 &cell,
1514 const unsigned int face_no,
1515 const unsigned int subface_no,
1516 const typename QProjector<dim>::DataSetDescriptor data_set,
1517 const Quadrature<dim - 1> &quadrature,
1518 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1520 &output_data)
1521 {
1522 const unsigned int n_q_points = quadrature.size();
1523
1524 maybe_compute_q_points<dim, spacedim>(data_set,
1525 data,
1526 output_data.quadrature_points,
1527 n_q_points);
1528 maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
1529 data_set,
1530 data,
1531 n_q_points);
1532 maybe_update_jacobian_grads<dim, spacedim>(CellSimilarity::none,
1533 data_set,
1534 data,
1535 output_data.jacobian_grads,
1536 n_q_points);
1537 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1539 data_set,
1540 data,
1542 n_q_points);
1543 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1545 data_set,
1546 data,
1547 output_data.jacobian_2nd_derivatives,
1548 n_q_points);
1549 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1551 data_set,
1552 data,
1554 n_q_points);
1555 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1557 data_set,
1558 data,
1559 output_data.jacobian_3rd_derivatives,
1560 n_q_points);
1561 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1563 data_set,
1564 data,
1566 n_q_points);
1567
1569 cell,
1570 face_no,
1571 subface_no,
1572 n_q_points,
1573 data_set,
1574 data,
1575 output_data);
1576 }
1577 } // namespace
1578 } // namespace MappingFEImplementation
1579} // namespace internal
1580
1581
1582
1583template <int dim, int spacedim>
1584void
1587 const unsigned int face_no,
1588 const hp::QCollection<dim - 1> &quadrature,
1589 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1591 &output_data) const
1592{
1593 // ensure that the following cast is really correct:
1594 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1596 const InternalData &data = static_cast<const InternalData &>(internal_data);
1597
1598 // if necessary, recompute the support points of the transformation of this
1599 // cell (note that we need to first check the triangulation pointer, since
1600 // otherwise the second test might trigger an exception if the
1601 // triangulations are not the same)
1602 if ((data.mapping_support_points.empty()) ||
1603 (&cell->get_triangulation() !=
1604 &data.cell_of_current_support_points->get_triangulation()) ||
1605 (cell != data.cell_of_current_support_points))
1606 {
1607 data.mapping_support_points = this->compute_mapping_support_points(cell);
1608 data.cell_of_current_support_points = cell;
1609 }
1610
1611 internal::MappingFEImplementation::do_fill_fe_face_values(
1612 *this,
1613 cell,
1614 face_no,
1616 QProjector<dim>::DataSetDescriptor::face(this->fe->reference_cell(),
1617 face_no,
1618 cell->combined_face_orientation(
1619 face_no),
1620 quadrature),
1621 quadrature[quadrature.size() == 1 ? 0 : face_no],
1622 data,
1623 output_data);
1624}
1625
1626
1627
1628template <int dim, int spacedim>
1629void
1632 const unsigned int face_no,
1633 const unsigned int subface_no,
1634 const Quadrature<dim - 1> &quadrature,
1635 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1637 &output_data) const
1638{
1639 // ensure that the following cast is really correct:
1640 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1642 const InternalData &data = static_cast<const InternalData &>(internal_data);
1643
1644 // if necessary, recompute the support points of the transformation of this
1645 // cell (note that we need to first check the triangulation pointer, since
1646 // otherwise the second test might trigger an exception if the
1647 // triangulations are not the same)
1648 if ((data.mapping_support_points.empty()) ||
1649 (&cell->get_triangulation() !=
1650 &data.cell_of_current_support_points->get_triangulation()) ||
1651 (cell != data.cell_of_current_support_points))
1652 {
1653 data.mapping_support_points = this->compute_mapping_support_points(cell);
1654 data.cell_of_current_support_points = cell;
1655 }
1656
1657 internal::MappingFEImplementation::do_fill_fe_face_values(
1658 *this,
1659 cell,
1660 face_no,
1661 subface_no,
1662 QProjector<dim>::DataSetDescriptor::subface(this->fe->reference_cell(),
1663 face_no,
1664 subface_no,
1665 cell->combined_face_orientation(
1666 face_no),
1667 quadrature.size(),
1668 cell->subface_case(face_no)),
1669 quadrature,
1670 data,
1671 output_data);
1672}
1673
1674
1675
1676namespace internal
1677{
1678 namespace MappingFEImplementation
1679 {
1680 namespace
1681 {
1682 template <int dim, int spacedim, int rank>
1683 void
1684 transform_fields(
1685 const ArrayView<const Tensor<rank, dim>> &input,
1686 const MappingKind mapping_kind,
1687 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1688 const ArrayView<Tensor<rank, spacedim>> &output)
1689 {
1690 // In the case of wedges and pyramids, faces might have different
1691 // numbers of quadrature points on each face with the result
1692 // that input and output have different sizes, since input has
1693 // the correct size but the size of output is the maximum of
1694 // all possible sizes.
1695 AssertIndexRange(input.size(), output.size() + 1);
1696
1697 Assert(
1698 (dynamic_cast<
1699 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1700 &mapping_data) != nullptr),
1702 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1703 static_cast<
1704 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1705 mapping_data);
1706
1707 switch (mapping_kind)
1708 {
1710 {
1711 Assert(
1714 "update_contravariant_transformation"));
1715
1716 for (unsigned int i = 0; i < input.size(); ++i)
1717 output[i] =
1718 apply_transformation(data.contravariant[i], input[i]);
1719
1720 return;
1721 }
1722
1723 case mapping_piola:
1724 {
1725 Assert(
1728 "update_contravariant_transformation"));
1729 Assert(
1730 data.update_each & update_volume_elements,
1732 "update_volume_elements"));
1733 Assert(rank == 1, ExcMessage("Only for rank 1"));
1734 if (rank != 1)
1735 return;
1736
1737 for (unsigned int i = 0; i < input.size(); ++i)
1738 {
1739 output[i] =
1740 apply_transformation(data.contravariant[i], input[i]);
1741 output[i] /= data.volume_elements[i];
1742 }
1743 return;
1744 }
1745 // We still allow this operation as in the
1746 // reference cell Derivatives are Tensor
1747 // rather than DerivativeForm
1748 case mapping_covariant:
1749 {
1750 Assert(
1753 "update_covariant_transformation"));
1754
1755 for (unsigned int i = 0; i < input.size(); ++i)
1756 output[i] = apply_transformation(data.covariant[i], input[i]);
1757
1758 return;
1759 }
1760
1761 default:
1763 }
1764 }
1765
1766
1767 template <int dim, int spacedim, int rank>
1768 void
1770 const ArrayView<const Tensor<rank, dim>> &input,
1771 const MappingKind mapping_kind,
1772 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1773 const ArrayView<Tensor<rank, spacedim>> &output)
1774 {
1775 AssertDimension(input.size(), output.size());
1776 Assert(
1777 (dynamic_cast<
1778 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1779 &mapping_data) != nullptr),
1781 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1782 static_cast<
1783 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1784 mapping_data);
1785
1786 switch (mapping_kind)
1787 {
1789 {
1790 Assert(
1793 "update_covariant_transformation"));
1794 Assert(
1797 "update_contravariant_transformation"));
1798 Assert(rank == 2, ExcMessage("Only for rank 2"));
1799
1800 for (unsigned int i = 0; i < output.size(); ++i)
1801 {
1803 apply_transformation(data.contravariant[i],
1804 transpose(input[i]));
1805 output[i] =
1806 apply_transformation(data.covariant[i], A.transpose());
1807 }
1808
1809 return;
1810 }
1811
1813 {
1814 Assert(
1817 "update_covariant_transformation"));
1818 Assert(rank == 2, ExcMessage("Only for rank 2"));
1819
1820 for (unsigned int i = 0; i < output.size(); ++i)
1821 {
1823 apply_transformation(data.covariant[i],
1824 transpose(input[i]));
1825 output[i] =
1826 apply_transformation(data.covariant[i], A.transpose());
1827 }
1828
1829 return;
1830 }
1831
1833 {
1834 Assert(
1837 "update_covariant_transformation"));
1838 Assert(
1841 "update_contravariant_transformation"));
1842 Assert(
1843 data.update_each & update_volume_elements,
1845 "update_volume_elements"));
1846 Assert(rank == 2, ExcMessage("Only for rank 2"));
1847
1848 for (unsigned int i = 0; i < output.size(); ++i)
1849 {
1851 apply_transformation(data.covariant[i], input[i]);
1852 const Tensor<2, spacedim> T =
1853 apply_transformation(data.contravariant[i],
1854 A.transpose());
1855
1856 output[i] = transpose(T);
1857 output[i] /= data.volume_elements[i];
1858 }
1859
1860 return;
1861 }
1862
1863 default:
1865 }
1866 }
1867
1868
1869
1870 template <int dim, int spacedim>
1871 void
1873 const ArrayView<const Tensor<3, dim>> &input,
1874 const MappingKind mapping_kind,
1875 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1876 const ArrayView<Tensor<3, spacedim>> &output)
1877 {
1878 AssertDimension(input.size(), output.size());
1879 Assert(
1880 (dynamic_cast<
1881 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1882 &mapping_data) != nullptr),
1884 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1885 static_cast<
1886 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1887 mapping_data);
1888
1889 switch (mapping_kind)
1890 {
1892 {
1893 Assert(
1896 "update_covariant_transformation"));
1897 Assert(
1900 "update_contravariant_transformation"));
1901
1902 for (unsigned int q = 0; q < output.size(); ++q)
1903 for (unsigned int i = 0; i < spacedim; ++i)
1904 {
1905 double tmp1[dim][dim];
1906 for (unsigned int J = 0; J < dim; ++J)
1907 for (unsigned int K = 0; K < dim; ++K)
1908 {
1909 tmp1[J][K] =
1910 data.contravariant[q][i][0] * input[q][0][J][K];
1911 for (unsigned int I = 1; I < dim; ++I)
1912 tmp1[J][K] +=
1913 data.contravariant[q][i][I] * input[q][I][J][K];
1914 }
1915 for (unsigned int j = 0; j < spacedim; ++j)
1916 {
1917 double tmp2[dim];
1918 for (unsigned int K = 0; K < dim; ++K)
1919 {
1920 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1921 for (unsigned int J = 1; J < dim; ++J)
1922 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1923 }
1924 for (unsigned int k = 0; k < spacedim; ++k)
1925 {
1926 output[q][i][j][k] =
1927 data.covariant[q][k][0] * tmp2[0];
1928 for (unsigned int K = 1; K < dim; ++K)
1929 output[q][i][j][k] +=
1930 data.covariant[q][k][K] * tmp2[K];
1931 }
1932 }
1933 }
1934 return;
1935 }
1936
1938 {
1939 Assert(
1942 "update_covariant_transformation"));
1943
1944 for (unsigned int q = 0; q < output.size(); ++q)
1945 for (unsigned int i = 0; i < spacedim; ++i)
1946 {
1947 double tmp1[dim][dim];
1948 for (unsigned int J = 0; J < dim; ++J)
1949 for (unsigned int K = 0; K < dim; ++K)
1950 {
1951 tmp1[J][K] =
1952 data.covariant[q][i][0] * input[q][0][J][K];
1953 for (unsigned int I = 1; I < dim; ++I)
1954 tmp1[J][K] +=
1955 data.covariant[q][i][I] * input[q][I][J][K];
1956 }
1957 for (unsigned int j = 0; j < spacedim; ++j)
1958 {
1959 double tmp2[dim];
1960 for (unsigned int K = 0; K < dim; ++K)
1961 {
1962 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1963 for (unsigned int J = 1; J < dim; ++J)
1964 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1965 }
1966 for (unsigned int k = 0; k < spacedim; ++k)
1967 {
1968 output[q][i][j][k] =
1969 data.covariant[q][k][0] * tmp2[0];
1970 for (unsigned int K = 1; K < dim; ++K)
1971 output[q][i][j][k] +=
1972 data.covariant[q][k][K] * tmp2[K];
1973 }
1974 }
1975 }
1976
1977 return;
1978 }
1979
1981 {
1982 Assert(
1985 "update_covariant_transformation"));
1986 Assert(
1989 "update_contravariant_transformation"));
1990 Assert(
1991 data.update_each & update_volume_elements,
1993 "update_volume_elements"));
1994
1995 for (unsigned int q = 0; q < output.size(); ++q)
1996 for (unsigned int i = 0; i < spacedim; ++i)
1997 {
1998 double factor[dim];
1999 for (unsigned int I = 0; I < dim; ++I)
2000 factor[I] =
2001 data.contravariant[q][i][I] / data.volume_elements[q];
2002 double tmp1[dim][dim];
2003 for (unsigned int J = 0; J < dim; ++J)
2004 for (unsigned int K = 0; K < dim; ++K)
2005 {
2006 tmp1[J][K] = factor[0] * input[q][0][J][K];
2007 for (unsigned int I = 1; I < dim; ++I)
2008 tmp1[J][K] += factor[I] * input[q][I][J][K];
2009 }
2010 for (unsigned int j = 0; j < spacedim; ++j)
2011 {
2012 double tmp2[dim];
2013 for (unsigned int K = 0; K < dim; ++K)
2014 {
2015 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2016 for (unsigned int J = 1; J < dim; ++J)
2017 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2018 }
2019 for (unsigned int k = 0; k < spacedim; ++k)
2020 {
2021 output[q][i][j][k] =
2022 data.covariant[q][k][0] * tmp2[0];
2023 for (unsigned int K = 1; K < dim; ++K)
2024 output[q][i][j][k] +=
2025 data.covariant[q][k][K] * tmp2[K];
2026 }
2027 }
2028 }
2029
2030 return;
2031 }
2032
2033 default:
2035 }
2036 }
2037
2038
2039
2040 template <int dim, int spacedim, int rank>
2041 void
2044 const MappingKind mapping_kind,
2045 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2047 {
2048 AssertDimension(input.size(), output.size());
2049 Assert(
2050 (dynamic_cast<
2051 const typename ::MappingFE<dim, spacedim>::InternalData *>(
2052 &mapping_data) != nullptr),
2054 const typename ::MappingFE<dim, spacedim>::InternalData &data =
2055 static_cast<
2056 const typename ::MappingFE<dim, spacedim>::InternalData &>(
2057 mapping_data);
2058
2059 switch (mapping_kind)
2060 {
2061 case mapping_covariant:
2062 {
2063 Assert(
2066 "update_covariant_transformation"));
2067
2068 for (unsigned int i = 0; i < output.size(); ++i)
2069 output[i] = apply_transformation(data.covariant[i], input[i]);
2070
2071 return;
2072 }
2073 default:
2075 }
2076 }
2077 } // namespace
2078 } // namespace MappingFEImplementation
2079} // namespace internal
2080
2081
2082
2083template <int dim, int spacedim>
2084void
2086 const ArrayView<const Tensor<1, dim>> &input,
2087 const MappingKind mapping_kind,
2088 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2089 const ArrayView<Tensor<1, spacedim>> &output) const
2090{
2091 internal::MappingFEImplementation::transform_fields(input,
2092 mapping_kind,
2093 mapping_data,
2094 output);
2095}
2096
2097
2098
2099template <int dim, int spacedim>
2100void
2102 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2103 const MappingKind mapping_kind,
2104 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2105 const ArrayView<Tensor<2, spacedim>> &output) const
2106{
2107 internal::MappingFEImplementation::transform_differential_forms(input,
2108 mapping_kind,
2109 mapping_data,
2110 output);
2111}
2112
2113
2114
2115template <int dim, int spacedim>
2116void
2118 const ArrayView<const Tensor<2, dim>> &input,
2119 const MappingKind mapping_kind,
2120 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2121 const ArrayView<Tensor<2, spacedim>> &output) const
2122{
2123 switch (mapping_kind)
2124 {
2126 internal::MappingFEImplementation::transform_fields(input,
2127 mapping_kind,
2128 mapping_data,
2129 output);
2130 return;
2131
2135 internal::MappingFEImplementation::transform_gradients(input,
2136 mapping_kind,
2137 mapping_data,
2138 output);
2139 return;
2140 default:
2142 }
2143}
2144
2145
2146
2147template <int dim, int spacedim>
2148void
2150 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2151 const MappingKind mapping_kind,
2152 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2153 const ArrayView<Tensor<3, spacedim>> &output) const
2154{
2155 AssertDimension(input.size(), output.size());
2156 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2158 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2159
2160 switch (mapping_kind)
2161 {
2163 {
2166 "update_covariant_transformation"));
2167
2168 for (unsigned int q = 0; q < output.size(); ++q)
2169 for (unsigned int i = 0; i < spacedim; ++i)
2170 for (unsigned int j = 0; j < spacedim; ++j)
2171 {
2172 double tmp[dim];
2173 for (unsigned int K = 0; K < dim; ++K)
2174 {
2175 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
2176 for (unsigned int J = 1; J < dim; ++J)
2177 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
2178 }
2179 for (unsigned int k = 0; k < spacedim; ++k)
2180 {
2181 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
2182 for (unsigned int K = 1; K < dim; ++K)
2183 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
2184 }
2185 }
2186 return;
2187 }
2188
2189 default:
2191 }
2192}
2193
2194
2195
2196template <int dim, int spacedim>
2197void
2199 const ArrayView<const Tensor<3, dim>> &input,
2200 const MappingKind mapping_kind,
2201 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2202 const ArrayView<Tensor<3, spacedim>> &output) const
2203{
2204 switch (mapping_kind)
2205 {
2209 internal::MappingFEImplementation::transform_hessians(input,
2210 mapping_kind,
2211 mapping_data,
2212 output);
2213 return;
2214 default:
2216 }
2217}
2218
2219
2220
2221namespace
2222{
2223 template <int spacedim>
2224 bool
2225 check_all_manifold_ids_identical(
2227 {
2228 return true;
2229 }
2230
2231
2232
2233 template <int spacedim>
2234 bool
2235 check_all_manifold_ids_identical(
2237 {
2238 const auto m_id = cell->manifold_id();
2239
2240 for (const auto f : cell->face_indices())
2241 if (m_id != cell->face(f)->manifold_id())
2242 return false;
2243
2244 return true;
2245 }
2246
2247
2248
2249 template <int spacedim>
2250 bool
2251 check_all_manifold_ids_identical(
2253 {
2254 const auto m_id = cell->manifold_id();
2255
2256 for (const auto f : cell->face_indices())
2257 if (m_id != cell->face(f)->manifold_id())
2258 return false;
2259
2260 for (const auto l : cell->line_indices())
2261 if (m_id != cell->line(l)->manifold_id())
2262 return false;
2263
2264 return true;
2265 }
2266} // namespace
2267
2268
2269
2270template <int dim, int spacedim>
2271std::vector<Point<spacedim>>
2273 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2274{
2275 std::vector<Point<spacedim>> mapping_support_points(
2276 fe->get_unit_support_points().size());
2277
2278 std::vector<Point<spacedim>> vertices(cell->n_vertices());
2279 for (const unsigned int i : cell->vertex_indices())
2280 vertices[i] = cell->vertex(i);
2281
2282 if (check_all_manifold_ids_identical(cell))
2283 {
2284 // version 1) all geometric entities have the same manifold:
2285
2286 cell->get_manifold().get_new_points(vertices,
2287 mapping_support_point_weights,
2288 mapping_support_points);
2289 }
2290 else
2291 {
2292 // version 1) geometric entities have different manifold
2293
2294 // helper function to compute mapped points on subentities
2295 // note[PM]: this function currently only uses the vertices
2296 // of cells to create new points on subentities; however,
2297 // one should use all bounding points to create new points
2298 // as in the case of MappingQ.
2299 const auto process = [&](const auto &manifold,
2300 const auto &indices,
2301 const unsigned n_points) {
2302 if ((indices.size() == 0) || (n_points == 0))
2303 return;
2304
2305 const unsigned int n_shape_functions =
2306 this->fe->reference_cell().n_vertices();
2307
2308 Table<2, double> mapping_support_point_weights_local(n_points,
2309 n_shape_functions);
2310 std::vector<Point<spacedim>> mapping_support_points_local(n_points);
2311
2312 for (unsigned int p = 0; p < n_points; ++p)
2313 for (unsigned int i = 0; i < n_shape_functions; ++i)
2314 mapping_support_point_weights_local(p, i) =
2315 mapping_support_point_weights(
2316 indices[p + (indices.size() - n_points)], i);
2317
2318 manifold.get_new_points(vertices,
2319 mapping_support_point_weights_local,
2320 mapping_support_points_local);
2321
2322 for (unsigned int p = 0; p < n_points; ++p)
2323 mapping_support_points[indices[p + (indices.size() - n_points)]] =
2324 mapping_support_points_local[p];
2325 };
2326
2327 // create dummy DoFHandler to extract indices on subobjects
2328 const auto &fe = *this->fe;
2330 GridGenerator::reference_cell(tria, fe.reference_cell());
2331 DoFHandler<dim, spacedim> dof_handler(tria);
2332 dof_handler.distribute_dofs(fe);
2333 const auto &cell_ref = dof_handler.begin_active();
2334
2335 std::vector<types::global_dof_index> indices;
2336
2337 // add vertices
2338 for (const unsigned int i : cell->vertex_indices())
2339 mapping_support_points[i] = cell->vertex(i);
2340
2341 // process and add line support points
2342 for (unsigned int l = 0; l < cell_ref->n_lines(); ++l)
2343 {
2344 const auto accessor = cell_ref->line(l);
2345 indices.resize(fe.n_dofs_per_line() + 2 * fe.n_dofs_per_vertex());
2346 accessor->get_dof_indices(indices);
2347 process(cell->line(l)->get_manifold(), indices, fe.n_dofs_per_line());
2348 }
2349
2350 // process and add face support points
2351 if constexpr (dim >= 3)
2352 {
2353 for (unsigned int f = 0; f < cell_ref->n_faces(); ++f)
2354 {
2355 const auto accessor = cell_ref->face(f);
2356 indices.resize(fe.n_dofs_per_face());
2357 accessor->get_dof_indices(indices);
2358 process(cell->face(f)->get_manifold(),
2359 indices,
2360 fe.n_dofs_per_quad());
2361 }
2362 }
2363
2364 // process and add volume support points
2365 if constexpr (dim >= 2)
2366 {
2367 indices.resize(fe.n_dofs_per_cell());
2368 cell_ref->get_dof_indices(indices);
2369 process(cell->get_manifold(),
2370 indices,
2371 (dim == 2) ? fe.n_dofs_per_quad() : fe.n_dofs_per_hex());
2372 }
2373 }
2374
2375 return mapping_support_points;
2376}
2377
2378
2379
2380template <int dim, int spacedim>
2383 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2384{
2385 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
2386}
2387
2388
2389
2390template <int dim, int spacedim>
2391bool
2393 const ReferenceCell &reference_cell) const
2394{
2395 Assert(dim == reference_cell.get_dimension(),
2396 ExcMessage("The dimension of your mapping (" +
2398 ") and the reference cell cell_type (" +
2399 Utilities::to_string(reference_cell.get_dimension()) +
2400 " ) do not agree."));
2401
2402 return fe->reference_cell() == reference_cell;
2403}
2404
2405
2406
2407//--------------------------- Explicit instantiations -----------------------
2408#include "fe/mapping_fe.inst"
2409
2410
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
DerivativeForm< 1, spacedim, dim, Number > transpose() const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
active_cell_iterator begin_active(const unsigned int level=0) 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
Definition mapping_fe.h:373
void compute_shape_function_values(const std::vector< Point< dim > > &unit_points)
virtual std::size_t memory_consumption() const override
Definition mapping_fe.cc:60
InternalData(const FiniteElement< dim, spacedim > &fe)
Definition mapping_fe.cc:49
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
Definition mapping_fe.cc:81
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
Definition mapping_fe.h:458
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
Definition mapping_fe.h:468
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
Definition mapping_fe.h:452
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
Definition qprojector.h:334
static DataSetDescriptor cell()
Definition qprojector.h:516
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
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
unsigned int size() const
Definition collection.h:316
unsigned int max_n_quadrature_points() const
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
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)
UpdateFlags
@ 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
MappingKind
Definition mapping.h:81
@ mapping_piola
Definition mapping.h:116
@ mapping_covariant_gradient
Definition mapping.h:102
@ mapping_covariant
Definition mapping.h:91
@ mapping_contravariant
Definition mapping.h:96
@ mapping_contravariant_hessian
Definition mapping.h:158
@ mapping_covariant_hessian
Definition mapping.h:152
@ mapping_contravariant_gradient
Definition mapping.h:108
@ mapping_piola_gradient
Definition mapping.h:122
@ mapping_piola_hessian
Definition mapping.h:164
std::vector< index_type > data
Definition mpi.cc:746
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
constexpr char T
constexpr char A
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)
Definition utilities.cc:193
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)
Definition utilities.cc:475
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
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)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int manifold_id
Definition types.h:173
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 > &)