deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50: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>
29
32#include <deal.II/grid/tria.h>
34
36
37#include <boost/container/small_vector.hpp>
38
39#include <algorithm>
40#include <array>
41#include <cmath>
42#include <memory>
43#include <numeric>
44
45
47
48
49template <int dim, int spacedim>
52 : fe(fe)
53 , polynomial_degree(fe.tensor_degree())
54 , n_shape_functions(fe.n_dofs_per_cell())
55{}
56
57
58
59template <int dim, int spacedim>
60std::size_t
77
78
79
80template <int dim, int spacedim>
81void
83 const Quadrature<dim> &q)
84{
85 // store the flags in the internal data object so we can access them
86 // in fill_fe_*_values()
87 this->update_each = update_flags;
88
89 const unsigned int n_q_points = q.size();
90
91 if (this->update_each & update_covariant_transformation)
92 covariant.resize(n_q_points);
93
94 if (this->update_each & update_contravariant_transformation)
95 contravariant.resize(n_q_points);
96
97 if (this->update_each & update_volume_elements)
98 volume_elements.resize(n_q_points);
99
100 // see if we need the (transformation) shape function values
101 // and/or gradients and resize the necessary arrays
102 if (this->update_each & update_quadrature_points)
103 shape_values.resize(n_shape_functions * n_q_points);
104
105 if (this->update_each &
113 shape_derivatives.resize(n_shape_functions * n_q_points);
114
115 if (this->update_each &
117 shape_second_derivatives.resize(n_shape_functions * n_q_points);
118
119 if (this->update_each & (update_jacobian_2nd_derivatives |
121 shape_third_derivatives.resize(n_shape_functions * n_q_points);
122
123 if (this->update_each & (update_jacobian_3rd_derivatives |
125 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
126
127 // now also fill the various fields with their correct values
128 compute_shape_function_values(q.get_points());
129
130 // copy (projected) quadrature weights
131 quadrature_weights = q.get_weights();
132}
133
134
135
136template <int dim, int spacedim>
137void
139 const UpdateFlags update_flags,
140 const Quadrature<dim> &q,
141 const unsigned int n_original_q_points)
142{
143 reinit(update_flags, q);
144
145 if (this->update_each &
148 {
149 aux.resize(dim - 1,
150 std::vector<Tensor<1, spacedim>>(n_original_q_points));
151
152 // Compute tangentials to the unit cell.
153 const auto reference_cell = this->fe.reference_cell();
154 const auto n_faces = reference_cell.n_faces();
155
156 for (unsigned int i = 0; i < n_faces; ++i)
157 {
158 unit_tangentials[i].resize(n_original_q_points);
159 std::fill(unit_tangentials[i].begin(),
160 unit_tangentials[i].end(),
161 reference_cell.template face_tangent_vector<dim>(i, 0));
162 if (dim > 2)
163 {
164 unit_tangentials[n_faces + i].resize(n_original_q_points);
165 std::fill(unit_tangentials[n_faces + i].begin(),
166 unit_tangentials[n_faces + i].end(),
167 reference_cell.template face_tangent_vector<dim>(i, 1));
168 }
169 }
170 }
171}
172
173
174
175template <int dim, int spacedim>
176void
178 const std::vector<Point<dim>> &unit_points)
179{
180 const auto fe_poly = dynamic_cast<const FE_Poly<dim, spacedim> *>(&this->fe);
181
182 Assert(fe_poly != nullptr, ExcNotImplemented());
183
184 const auto &tensor_pols = fe_poly->get_poly_space();
185
186 const unsigned int n_shape_functions = fe.n_dofs_per_cell();
187 const unsigned int n_points = unit_points.size();
188
189 std::vector<double> values;
190 std::vector<Tensor<1, dim>> grads;
191 if (shape_values.size() != 0)
192 {
193 Assert(shape_values.size() == n_shape_functions * n_points,
195 values.resize(n_shape_functions);
196 }
197 if (shape_derivatives.size() != 0)
198 {
199 Assert(shape_derivatives.size() == n_shape_functions * n_points,
201 grads.resize(n_shape_functions);
202 }
203
204 std::vector<Tensor<2, dim>> grad2;
205 if (shape_second_derivatives.size() != 0)
206 {
207 Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
209 grad2.resize(n_shape_functions);
210 }
211
212 std::vector<Tensor<3, dim>> grad3;
213 if (shape_third_derivatives.size() != 0)
214 {
215 Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
217 grad3.resize(n_shape_functions);
218 }
219
220 std::vector<Tensor<4, dim>> grad4;
221 if (shape_fourth_derivatives.size() != 0)
222 {
223 Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
225 grad4.resize(n_shape_functions);
226 }
227
228
229 if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
230 shape_second_derivatives.size() != 0 ||
231 shape_third_derivatives.size() != 0 ||
232 shape_fourth_derivatives.size() != 0)
233 for (unsigned int point = 0; point < n_points; ++point)
234 {
235 tensor_pols.evaluate(
236 unit_points[point], values, grads, grad2, grad3, grad4);
237
238 if (shape_values.size() != 0)
239 for (unsigned int i = 0; i < n_shape_functions; ++i)
240 shape(point, i) = values[i];
241
242 if (shape_derivatives.size() != 0)
243 for (unsigned int i = 0; i < n_shape_functions; ++i)
244 derivative(point, i) = grads[i];
245
246 if (shape_second_derivatives.size() != 0)
247 for (unsigned int i = 0; i < n_shape_functions; ++i)
248 second_derivative(point, i) = grad2[i];
249
250 if (shape_third_derivatives.size() != 0)
251 for (unsigned int i = 0; i < n_shape_functions; ++i)
252 third_derivative(point, i) = grad3[i];
253
254 if (shape_fourth_derivatives.size() != 0)
255 for (unsigned int i = 0; i < n_shape_functions; ++i)
256 fourth_derivative(point, i) = grad4[i];
257 }
258}
259
260
261namespace internal
262{
263 namespace MappingFEImplementation
264 {
265 namespace
266 {
273 template <int dim, int spacedim>
274 void
275 maybe_compute_q_points(
276 const typename QProjector<dim>::DataSetDescriptor data_set,
277 const typename ::MappingFE<dim, spacedim>::InternalData &data,
278 std::vector<Point<spacedim>> &quadrature_points,
279 const unsigned int n_q_points)
280 {
281 const UpdateFlags update_flags = data.update_each;
282
283 if (update_flags & update_quadrature_points)
284 for (unsigned int point = 0; point < n_q_points; ++point)
285 {
286 const double *shape = &data.shape(point + data_set, 0);
287 Point<spacedim> result =
288 (shape[0] * data.mapping_support_points[0]);
289 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
290 for (unsigned int i = 0; i < spacedim; ++i)
291 result[i] += shape[k] * data.mapping_support_points[k][i];
292 quadrature_points[point] = result;
293 }
294 }
295
296
297
306 template <int dim, int spacedim>
307 void
308 maybe_update_Jacobians(
309 const CellSimilarity::Similarity cell_similarity,
310 const typename ::QProjector<dim>::DataSetDescriptor data_set,
311 const typename ::MappingFE<dim, spacedim>::InternalData &data,
312 const unsigned int n_q_points)
313 {
314 const UpdateFlags update_flags = data.update_each;
315
316 if (update_flags & update_contravariant_transformation)
317 // if the current cell is just a
318 // translation of the previous one, no
319 // need to recompute jacobians...
320 if (cell_similarity != CellSimilarity::translation)
321 {
322 std::fill(data.contravariant.begin(),
323 data.contravariant.end(),
325
326 Assert(data.n_shape_functions > 0, ExcInternalError());
327
328 for (unsigned int point = 0; point < n_q_points; ++point)
329 {
330 double result[spacedim][dim];
331
332 // peel away part of sum to avoid zeroing the
333 // entries and adding for the first time
334 for (unsigned int i = 0; i < spacedim; ++i)
335 for (unsigned int j = 0; j < dim; ++j)
336 result[i][j] = data.derivative(point + data_set, 0)[j] *
337 data.mapping_support_points[0][i];
338 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
339 for (unsigned int i = 0; i < spacedim; ++i)
340 for (unsigned int j = 0; j < dim; ++j)
341 result[i][j] +=
342 data.derivative(point + data_set, k)[j] *
343 data.mapping_support_points[k][i];
344
345 // write result into contravariant data. for
346 // j=dim in the case dim<spacedim, there will
347 // never be any nonzero data that arrives in
348 // here, so it is ok anyway because it was
349 // initialized to zero at the initialization
350 for (unsigned int i = 0; i < spacedim; ++i)
351 for (unsigned int j = 0; j < dim; ++j)
352 data.contravariant[point][i][j] = result[i][j];
353 }
354 }
355
356 if (update_flags & update_covariant_transformation)
357 if (cell_similarity != CellSimilarity::translation)
358 {
359 for (unsigned int point = 0; point < n_q_points; ++point)
360 {
361 data.covariant[point] =
362 (data.contravariant[point]).covariant_form();
363 }
364 }
365
366 if (update_flags & update_volume_elements)
367 if (cell_similarity != CellSimilarity::translation)
368 {
369 for (unsigned int point = 0; point < n_q_points; ++point)
370 data.volume_elements[point] =
371 data.contravariant[point].determinant();
372 }
373 }
374
381 template <int dim, int spacedim>
382 void
383 maybe_update_jacobian_grads(
384 const CellSimilarity::Similarity cell_similarity,
385 const typename QProjector<dim>::DataSetDescriptor data_set,
386 const typename ::MappingFE<dim, spacedim>::InternalData &data,
387 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads,
388 const unsigned int n_q_points)
389 {
390 const UpdateFlags update_flags = data.update_each;
391 if (update_flags & update_jacobian_grads)
392 {
393 AssertIndexRange(n_q_points, jacobian_grads.size() + 1);
394
395 if (cell_similarity != CellSimilarity::translation)
396 for (unsigned int point = 0; point < n_q_points; ++point)
397 {
398 const Tensor<2, dim> *second =
399 &data.second_derivative(point + data_set, 0);
400 double result[spacedim][dim][dim];
401 for (unsigned int i = 0; i < spacedim; ++i)
402 for (unsigned int j = 0; j < dim; ++j)
403 for (unsigned int l = 0; l < dim; ++l)
404 result[i][j][l] =
405 (second[0][j][l] * data.mapping_support_points[0][i]);
406 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
407 for (unsigned int i = 0; i < spacedim; ++i)
408 for (unsigned int j = 0; j < dim; ++j)
409 for (unsigned int l = 0; l < dim; ++l)
410 result[i][j][l] +=
411 (second[k][j][l] *
412 data.mapping_support_points[k][i]);
413
414 for (unsigned int i = 0; i < spacedim; ++i)
415 for (unsigned int j = 0; j < dim; ++j)
416 for (unsigned int l = 0; l < dim; ++l)
417 jacobian_grads[point][i][j][l] = result[i][j][l];
418 }
419 }
420 }
421
428 template <int dim, int spacedim>
429 void
430 maybe_update_jacobian_pushed_forward_grads(
431 const CellSimilarity::Similarity cell_similarity,
432 const typename QProjector<dim>::DataSetDescriptor data_set,
433 const typename ::MappingFE<dim, spacedim>::InternalData &data,
434 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads,
435 const unsigned int n_q_points)
436 {
437 const UpdateFlags update_flags = data.update_each;
438 if (update_flags & update_jacobian_pushed_forward_grads)
439 {
440 AssertIndexRange(n_q_points,
441 jacobian_pushed_forward_grads.size() + 1);
442
443 if (cell_similarity != CellSimilarity::translation)
444 {
445 double tmp[spacedim][spacedim][spacedim];
446 for (unsigned int point = 0; point < n_q_points; ++point)
447 {
448 const Tensor<2, dim> *second =
449 &data.second_derivative(point + data_set, 0);
450 double result[spacedim][dim][dim];
451 for (unsigned int i = 0; i < spacedim; ++i)
452 for (unsigned int j = 0; j < dim; ++j)
453 for (unsigned int l = 0; l < dim; ++l)
454 result[i][j][l] = (second[0][j][l] *
455 data.mapping_support_points[0][i]);
456 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
457 for (unsigned int i = 0; i < spacedim; ++i)
458 for (unsigned int j = 0; j < dim; ++j)
459 for (unsigned int l = 0; l < dim; ++l)
460 result[i][j][l] +=
461 (second[k][j][l] *
462 data.mapping_support_points[k][i]);
463
464 // first push forward the j-components
465 for (unsigned int i = 0; i < spacedim; ++i)
466 for (unsigned int j = 0; j < spacedim; ++j)
467 for (unsigned int l = 0; l < dim; ++l)
468 {
469 tmp[i][j][l] =
470 result[i][0][l] * data.covariant[point][j][0];
471 for (unsigned int jr = 1; jr < dim; ++jr)
472 {
473 tmp[i][j][l] += result[i][jr][l] *
474 data.covariant[point][j][jr];
475 }
476 }
477
478 // now, pushing forward the l-components
479 for (unsigned int i = 0; i < spacedim; ++i)
480 for (unsigned int j = 0; j < spacedim; ++j)
481 for (unsigned int l = 0; l < spacedim; ++l)
482 {
483 jacobian_pushed_forward_grads[point][i][j][l] =
484 tmp[i][j][0] * data.covariant[point][l][0];
485 for (unsigned int lr = 1; lr < dim; ++lr)
486 {
487 jacobian_pushed_forward_grads[point][i][j][l] +=
488 tmp[i][j][lr] * data.covariant[point][l][lr];
489 }
490 }
491 }
492 }
493 }
494 }
495
502 template <int dim, int spacedim>
503 void
504 maybe_update_jacobian_2nd_derivatives(
505 const CellSimilarity::Similarity cell_similarity,
506 const typename QProjector<dim>::DataSetDescriptor data_set,
507 const typename ::MappingFE<dim, spacedim>::InternalData &data,
508 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives,
509 const unsigned int n_q_points)
510 {
511 const UpdateFlags update_flags = data.update_each;
512 if (update_flags & update_jacobian_2nd_derivatives)
513 {
514 AssertIndexRange(n_q_points, jacobian_2nd_derivatives.size() + 1);
515
516 if (cell_similarity != CellSimilarity::translation)
517 {
518 for (unsigned int point = 0; point < n_q_points; ++point)
519 {
520 const Tensor<3, dim> *third =
521 &data.third_derivative(point + data_set, 0);
522 double result[spacedim][dim][dim][dim];
523 for (unsigned int i = 0; i < spacedim; ++i)
524 for (unsigned int j = 0; j < dim; ++j)
525 for (unsigned int l = 0; l < dim; ++l)
526 for (unsigned int m = 0; m < dim; ++m)
527 result[i][j][l][m] =
528 (third[0][j][l][m] *
529 data.mapping_support_points[0][i]);
530 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
531 for (unsigned int i = 0; i < spacedim; ++i)
532 for (unsigned int j = 0; j < dim; ++j)
533 for (unsigned int l = 0; l < dim; ++l)
534 for (unsigned int m = 0; m < dim; ++m)
535 result[i][j][l][m] +=
536 (third[k][j][l][m] *
537 data.mapping_support_points[k][i]);
538
539 for (unsigned int i = 0; i < spacedim; ++i)
540 for (unsigned int j = 0; j < dim; ++j)
541 for (unsigned int l = 0; l < dim; ++l)
542 for (unsigned int m = 0; m < dim; ++m)
543 jacobian_2nd_derivatives[point][i][j][l][m] =
544 result[i][j][l][m];
545 }
546 }
547 }
548 }
549
557 template <int dim, int spacedim>
558 void
559 maybe_update_jacobian_pushed_forward_2nd_derivatives(
560 const CellSimilarity::Similarity cell_similarity,
561 const typename QProjector<dim>::DataSetDescriptor data_set,
562 const typename ::MappingFE<dim, spacedim>::InternalData &data,
563 std::vector<Tensor<4, spacedim>>
564 &jacobian_pushed_forward_2nd_derivatives,
565 const unsigned int n_q_points)
566 {
567 const UpdateFlags update_flags = data.update_each;
569 {
570 AssertIndexRange(n_q_points,
571 jacobian_pushed_forward_2nd_derivatives.size() +
572 1);
573
574 if (cell_similarity != CellSimilarity::translation)
575 {
576 double tmp[spacedim][spacedim][spacedim][spacedim];
577 for (unsigned int point = 0; point < n_q_points; ++point)
578 {
579 const Tensor<3, dim> *third =
580 &data.third_derivative(point + data_set, 0);
581 double result[spacedim][dim][dim][dim];
582 for (unsigned int i = 0; i < spacedim; ++i)
583 for (unsigned int j = 0; j < dim; ++j)
584 for (unsigned int l = 0; l < dim; ++l)
585 for (unsigned int m = 0; m < dim; ++m)
586 result[i][j][l][m] =
587 (third[0][j][l][m] *
588 data.mapping_support_points[0][i]);
589 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
590 for (unsigned int i = 0; i < spacedim; ++i)
591 for (unsigned int j = 0; j < dim; ++j)
592 for (unsigned int l = 0; l < dim; ++l)
593 for (unsigned int m = 0; m < dim; ++m)
594 result[i][j][l][m] +=
595 (third[k][j][l][m] *
596 data.mapping_support_points[k][i]);
597
598 // push forward the j-coordinate
599 for (unsigned int i = 0; i < spacedim; ++i)
600 for (unsigned int j = 0; j < spacedim; ++j)
601 for (unsigned int l = 0; l < dim; ++l)
602 for (unsigned int m = 0; m < dim; ++m)
603 {
604 jacobian_pushed_forward_2nd_derivatives
605 [point][i][j][l][m] =
606 result[i][0][l][m] *
607 data.covariant[point][j][0];
608 for (unsigned int jr = 1; jr < dim; ++jr)
609 jacobian_pushed_forward_2nd_derivatives[point]
610 [i][j][l]
611 [m] +=
612 result[i][jr][l][m] *
613 data.covariant[point][j][jr];
614 }
615
616 // push forward the l-coordinate
617 for (unsigned int i = 0; i < spacedim; ++i)
618 for (unsigned int j = 0; j < spacedim; ++j)
619 for (unsigned int l = 0; l < spacedim; ++l)
620 for (unsigned int m = 0; m < dim; ++m)
621 {
622 tmp[i][j][l][m] =
623 jacobian_pushed_forward_2nd_derivatives[point]
624 [i][j][0]
625 [m] *
626 data.covariant[point][l][0];
627 for (unsigned int lr = 1; lr < dim; ++lr)
628 tmp[i][j][l][m] +=
629 jacobian_pushed_forward_2nd_derivatives
630 [point][i][j][lr][m] *
631 data.covariant[point][l][lr];
632 }
633
634 // push forward the m-coordinate
635 for (unsigned int i = 0; i < spacedim; ++i)
636 for (unsigned int j = 0; j < spacedim; ++j)
637 for (unsigned int l = 0; l < spacedim; ++l)
638 for (unsigned int m = 0; m < spacedim; ++m)
639 {
640 jacobian_pushed_forward_2nd_derivatives
641 [point][i][j][l][m] =
642 tmp[i][j][l][0] * data.covariant[point][m][0];
643 for (unsigned int mr = 1; mr < dim; ++mr)
644 jacobian_pushed_forward_2nd_derivatives[point]
645 [i][j][l]
646 [m] +=
647 tmp[i][j][l][mr] *
648 data.covariant[point][m][mr];
649 }
650 }
651 }
652 }
653 }
654
661 template <int dim, int spacedim>
662 void
663 maybe_update_jacobian_3rd_derivatives(
664 const CellSimilarity::Similarity cell_similarity,
665 const typename QProjector<dim>::DataSetDescriptor data_set,
666 const typename ::MappingFE<dim, spacedim>::InternalData &data,
667 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives,
668 const unsigned int n_q_points)
669 {
670 const UpdateFlags update_flags = data.update_each;
671 if (update_flags & update_jacobian_3rd_derivatives)
672 {
673 AssertIndexRange(n_q_points, jacobian_3rd_derivatives.size() + 1);
674
675 if (cell_similarity != CellSimilarity::translation)
676 {
677 for (unsigned int point = 0; point < n_q_points; ++point)
678 {
679 const Tensor<4, dim> *fourth =
680 &data.fourth_derivative(point + data_set, 0);
681 double result[spacedim][dim][dim][dim][dim];
682 for (unsigned int i = 0; i < spacedim; ++i)
683 for (unsigned int j = 0; j < dim; ++j)
684 for (unsigned int l = 0; l < dim; ++l)
685 for (unsigned int m = 0; m < dim; ++m)
686 for (unsigned int n = 0; n < dim; ++n)
687 result[i][j][l][m][n] =
688 (fourth[0][j][l][m][n] *
689 data.mapping_support_points[0][i]);
690 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
691 for (unsigned int i = 0; i < spacedim; ++i)
692 for (unsigned int j = 0; j < dim; ++j)
693 for (unsigned int l = 0; l < dim; ++l)
694 for (unsigned int m = 0; m < dim; ++m)
695 for (unsigned int n = 0; n < dim; ++n)
696 result[i][j][l][m][n] +=
697 (fourth[k][j][l][m][n] *
698 data.mapping_support_points[k][i]);
699
700 for (unsigned int i = 0; i < spacedim; ++i)
701 for (unsigned int j = 0; j < dim; ++j)
702 for (unsigned int l = 0; l < dim; ++l)
703 for (unsigned int m = 0; m < dim; ++m)
704 for (unsigned int n = 0; n < dim; ++n)
705 jacobian_3rd_derivatives[point][i][j][l][m][n] =
706 result[i][j][l][m][n];
707 }
708 }
709 }
710 }
711
719 template <int dim, int spacedim>
720 void
721 maybe_update_jacobian_pushed_forward_3rd_derivatives(
722 const CellSimilarity::Similarity cell_similarity,
723 const typename QProjector<dim>::DataSetDescriptor data_set,
724 const typename ::MappingFE<dim, spacedim>::InternalData &data,
725 std::vector<Tensor<5, spacedim>>
726 &jacobian_pushed_forward_3rd_derivatives,
727 const unsigned int n_q_points)
728 {
729 const UpdateFlags update_flags = data.update_each;
731 {
732 AssertIndexRange(n_q_points,
733 jacobian_pushed_forward_3rd_derivatives.size() +
734 1);
735
736 if (cell_similarity != CellSimilarity::translation)
737 {
738 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
739 for (unsigned int point = 0; point < n_q_points; ++point)
740 {
741 const Tensor<4, dim> *fourth =
742 &data.fourth_derivative(point + data_set, 0);
743 double result[spacedim][dim][dim][dim][dim];
744 for (unsigned int i = 0; i < spacedim; ++i)
745 for (unsigned int j = 0; j < dim; ++j)
746 for (unsigned int l = 0; l < dim; ++l)
747 for (unsigned int m = 0; m < dim; ++m)
748 for (unsigned int n = 0; n < dim; ++n)
749 result[i][j][l][m][n] =
750 (fourth[0][j][l][m][n] *
751 data.mapping_support_points[0][i]);
752 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
753 for (unsigned int i = 0; i < spacedim; ++i)
754 for (unsigned int j = 0; j < dim; ++j)
755 for (unsigned int l = 0; l < dim; ++l)
756 for (unsigned int m = 0; m < dim; ++m)
757 for (unsigned int n = 0; n < dim; ++n)
758 result[i][j][l][m][n] +=
759 (fourth[k][j][l][m][n] *
760 data.mapping_support_points[k][i]);
761
762 // push-forward the j-coordinate
763 for (unsigned int i = 0; i < spacedim; ++i)
764 for (unsigned int j = 0; j < spacedim; ++j)
765 for (unsigned int l = 0; l < dim; ++l)
766 for (unsigned int m = 0; m < dim; ++m)
767 for (unsigned int n = 0; n < dim; ++n)
768 {
769 tmp[i][j][l][m][n] =
770 result[i][0][l][m][n] *
771 data.covariant[point][j][0];
772 for (unsigned int jr = 1; jr < dim; ++jr)
773 tmp[i][j][l][m][n] +=
774 result[i][jr][l][m][n] *
775 data.covariant[point][j][jr];
776 }
777
778 // push-forward the l-coordinate
779 for (unsigned int i = 0; i < spacedim; ++i)
780 for (unsigned int j = 0; j < spacedim; ++j)
781 for (unsigned int l = 0; l < spacedim; ++l)
782 for (unsigned int m = 0; m < dim; ++m)
783 for (unsigned int n = 0; n < dim; ++n)
784 {
785 jacobian_pushed_forward_3rd_derivatives
786 [point][i][j][l][m][n] =
787 tmp[i][j][0][m][n] *
788 data.covariant[point][l][0];
789 for (unsigned int lr = 1; lr < dim; ++lr)
790 jacobian_pushed_forward_3rd_derivatives
791 [point][i][j][l][m][n] +=
792 tmp[i][j][lr][m][n] *
793 data.covariant[point][l][lr];
794 }
795
796 // push-forward the m-coordinate
797 for (unsigned int i = 0; i < spacedim; ++i)
798 for (unsigned int j = 0; j < spacedim; ++j)
799 for (unsigned int l = 0; l < spacedim; ++l)
800 for (unsigned int m = 0; m < spacedim; ++m)
801 for (unsigned int n = 0; n < dim; ++n)
802 {
803 tmp[i][j][l][m][n] =
804 jacobian_pushed_forward_3rd_derivatives
805 [point][i][j][l][0][n] *
806 data.covariant[point][m][0];
807 for (unsigned int mr = 1; mr < dim; ++mr)
808 tmp[i][j][l][m][n] +=
809 jacobian_pushed_forward_3rd_derivatives
810 [point][i][j][l][mr][n] *
811 data.covariant[point][m][mr];
812 }
813
814 // push-forward the n-coordinate
815 for (unsigned int i = 0; i < spacedim; ++i)
816 for (unsigned int j = 0; j < spacedim; ++j)
817 for (unsigned int l = 0; l < spacedim; ++l)
818 for (unsigned int m = 0; m < spacedim; ++m)
819 for (unsigned int n = 0; n < spacedim; ++n)
820 {
821 jacobian_pushed_forward_3rd_derivatives
822 [point][i][j][l][m][n] =
823 tmp[i][j][l][m][0] *
824 data.covariant[point][n][0];
825 for (unsigned int nr = 1; nr < dim; ++nr)
826 jacobian_pushed_forward_3rd_derivatives
827 [point][i][j][l][m][n] +=
828 tmp[i][j][l][m][nr] *
829 data.covariant[point][n][nr];
830 }
831 }
832 }
833 }
834 }
835 } // namespace
836 } // namespace MappingFEImplementation
837} // namespace internal
838
839
840
841template <int dim, int spacedim>
843 : fe(fe.clone())
844 , polynomial_degree(fe.tensor_degree())
845{
847 ExcMessage("It only makes sense to create polynomial mappings "
848 "with a polynomial degree greater or equal to one."));
849 Assert(fe.n_components() == 1, ExcNotImplemented());
850
851 Assert(fe.has_support_points(), ExcNotImplemented());
852
853 const auto &mapping_support_points = fe.get_unit_support_points();
854
855 const auto reference_cell = fe.reference_cell();
856
857 const unsigned int n_points = mapping_support_points.size();
858 const unsigned int n_shape_functions = reference_cell.n_vertices();
859
861 Table<2, double>(n_points, n_shape_functions);
862
863 for (unsigned int point = 0; point < n_points; ++point)
864 for (unsigned int i = 0; i < n_shape_functions; ++i)
866 reference_cell.d_linear_shape_function(mapping_support_points[point],
867 i);
868}
869
870
871
872template <int dim, int spacedim>
874 : fe(mapping.fe->clone())
875 , polynomial_degree(mapping.polynomial_degree)
876 , mapping_support_point_weights(mapping.mapping_support_point_weights)
877{}
878
879
880
881template <int dim, int spacedim>
882std::unique_ptr<Mapping<dim, spacedim>>
884{
885 return std::make_unique<MappingFE<dim, spacedim>>(*this);
886}
887
888
889
890template <int dim, int spacedim>
891unsigned int
893{
894 return polynomial_degree;
895}
896
897
898
899template <int dim, int spacedim>
903 const Point<dim> &p) const
904{
905 const std::vector<Point<spacedim>> support_points =
906 this->compute_mapping_support_points(cell);
907
908 Point<spacedim> mapped_point;
909
910 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
911 mapped_point += support_points[i] * this->fe->shape_value(i, p);
912
913 return mapped_point;
914}
915
916
917
918template <int dim, int spacedim>
922 const Point<spacedim> &p) const
923{
924 const std::vector<Point<spacedim>> support_points =
925 this->compute_mapping_support_points(cell);
926
927 const double eps = 1.e-12 * cell->diameter();
928 const unsigned int loop_limit = 10;
929
930 Point<dim> p_unit;
931
932 unsigned int loop = 0;
933
934 // This loop solves the following problem:
935 // grad_F^T residual + (grad_F^T grad_F + grad_F^T hess_F^T dp) dp = 0
936 // where the term
937 // (grad_F^T hess_F dp) is approximated by (-hess_F * residual)
938 // This is basically a second order approximation of Newton method, where the
939 // Jacobian is corrected with a higher order term coming from the hessian.
940 do
941 {
942 Point<spacedim> mapped_point;
943
944 // Transpose of the gradient map
947
948 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
949 {
950 mapped_point += support_points[i] * this->fe->shape_value(i, p_unit);
951 const auto grad_F_i = this->fe->shape_grad(i, p_unit);
952 const auto hessian_F_i = this->fe->shape_grad_grad(i, p_unit);
953 for (unsigned int j = 0; j < dim; ++j)
954 {
955 grad_FT[j] += grad_F_i[j] * support_points[i];
956 for (unsigned int l = 0; l < dim; ++l)
957 hess_FT[j][l] += hessian_F_i[j][l] * support_points[i];
958 }
959 }
960
961 // Residual
962 const auto residual = p - mapped_point;
963 // Project the residual on the reference coordinate system
964 // to compute the error, and to filter components orthogonal to the
965 // manifold, and compute a 2nd order correction of the metric tensor
966 const auto grad_FT_residual = apply_transformation(grad_FT, residual);
967
968 // Do not invert nor compute the metric if not necessary.
969 if (grad_FT_residual.norm() <= eps)
970 break;
971
972 // Now compute the (corrected) metric tensor
973 Tensor<2, dim> corrected_metric_tensor;
974 for (unsigned int j = 0; j < dim; ++j)
975 for (unsigned int l = 0; l < dim; ++l)
976 corrected_metric_tensor[j][l] =
977 -grad_FT[j] * grad_FT[l] + hess_FT[j][l] * residual;
978
979 // And compute the update
980 const auto g_inverse = invert(corrected_metric_tensor);
981 p_unit -= Point<dim>(g_inverse * grad_FT_residual);
982
983 ++loop;
984 }
985 while (loop < loop_limit);
986
987 // Here we check that in the last execution of while the first
988 // condition was already wrong, meaning the residual was below
989 // eps. Only if the first condition failed, loop will have been
990 // increased and tested, and thus have reached the limit.
991 AssertThrow(loop < loop_limit,
993
994 return p_unit;
995}
996
997
998
999template <int dim, int spacedim>
1002{
1003 // add flags if the respective quantities are necessary to compute
1004 // what we need. note that some flags appear in both the conditions
1005 // and in subsequent set operations. this leads to some circular
1006 // logic. the only way to treat this is to iterate. since there are
1007 // 5 if-clauses in the loop, it will take at most 5 iterations to
1008 // converge. do them:
1009 UpdateFlags out = in;
1010 for (unsigned int i = 0; i < 5; ++i)
1011 {
1012 // The following is a little incorrect:
1013 // If not applied on a face,
1014 // update_boundary_forms does not
1015 // make sense. On the other hand,
1016 // it is necessary on a
1017 // face. Currently,
1018 // update_boundary_forms is simply
1019 // ignored for the interior of a
1020 // cell.
1022 out |= update_boundary_forms;
1023
1028
1029 if (out &
1034
1035 // The contravariant transformation is used in the Piola
1036 // transformation, which requires the determinant of the Jacobi
1037 // matrix of the transformation. Because we have no way of
1038 // knowing here whether the finite element wants to use the
1039 // contravariant or the Piola transforms, we add the JxW values
1040 // to the list of flags to be updated for each cell.
1043
1044 // the same is true when computing normal vectors: they require
1045 // the determinant of the Jacobian
1046 if (out & update_normal_vectors)
1048 }
1049
1050 return out;
1051}
1052
1053
1054
1055template <int dim, int spacedim>
1056std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1058 const Quadrature<dim> &q) const
1059{
1060 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1061 std::make_unique<InternalData>(*this->fe);
1062 data_ptr->reinit(this->requires_update_flags(update_flags), q);
1063 return data_ptr;
1064}
1065
1066
1067
1068template <int dim, int spacedim>
1069std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1071 const UpdateFlags update_flags,
1072 const hp::QCollection<dim - 1> &quadrature) const
1073{
1074 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1075 std::make_unique<InternalData>(*this->fe);
1076 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1077 data.initialize_face(this->requires_update_flags(update_flags),
1079 this->fe->reference_cell(), quadrature),
1080 quadrature.max_n_quadrature_points());
1081
1082 return data_ptr;
1083}
1084
1085
1086
1087template <int dim, int spacedim>
1088std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1090 const UpdateFlags update_flags,
1091 const Quadrature<dim - 1> &quadrature) const
1092{
1093 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1094 std::make_unique<InternalData>(*this->fe);
1095 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1096 data.initialize_face(this->requires_update_flags(update_flags),
1098 this->fe->reference_cell(), quadrature),
1099 quadrature.size());
1100
1101 return data_ptr;
1102}
1103
1104
1105
1106template <int dim, int spacedim>
1110 const CellSimilarity::Similarity cell_similarity,
1111 const Quadrature<dim> &quadrature,
1112 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1114 &output_data) const
1115{
1116 // ensure that the following static_cast is really correct:
1117 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1119 const InternalData &data = static_cast<const InternalData &>(internal_data);
1120
1121 const unsigned int n_q_points = quadrature.size();
1122
1123 // recompute the support points of the transformation of this
1124 // cell. we tried to be clever here in an earlier version of the
1125 // library by checking whether the cell is the same as the one we
1126 // had visited last, but it turns out to be difficult to determine
1127 // that because a cell for the purposes of a mapping is
1128 // characterized not just by its (triangulation, level, index)
1129 // triple, but also by the locations of its vertices, the manifold
1130 // object attached to the cell and all of its bounding faces/edges,
1131 // etc. to reliably test that the "cell" we are on is, therefore,
1132 // not easily done
1133 data.mapping_support_points = this->compute_mapping_support_points(cell);
1134 data.cell_of_current_support_points = cell;
1135
1136 // if the order of the mapping is greater than 1, then do not reuse any cell
1137 // similarity information. This is necessary because the cell similarity
1138 // value is computed with just cell vertices and does not take into account
1139 // cell curvature.
1140 const CellSimilarity::Similarity computed_cell_similarity =
1141 (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
1142
1143 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1145 data,
1146 output_data.quadrature_points,
1147 n_q_points);
1148
1149 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1150 computed_cell_similarity,
1152 data,
1153 n_q_points);
1154
1155 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1156 computed_cell_similarity,
1158 data,
1159 output_data.jacobian_grads,
1160 n_q_points);
1161
1162 internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
1163 dim,
1164 spacedim>(computed_cell_similarity,
1166 data,
1168 n_q_points);
1169
1170 internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
1171 dim,
1172 spacedim>(computed_cell_similarity,
1174 data,
1175 output_data.jacobian_2nd_derivatives,
1176 n_q_points);
1177
1178 internal::MappingFEImplementation::
1179 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1180 computed_cell_similarity,
1182 data,
1184 n_q_points);
1185
1186 internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
1187 dim,
1188 spacedim>(computed_cell_similarity,
1190 data,
1191 output_data.jacobian_3rd_derivatives,
1192 n_q_points);
1193
1194 internal::MappingFEImplementation::
1195 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1196 computed_cell_similarity,
1198 data,
1200 n_q_points);
1201
1202 const UpdateFlags update_flags = data.update_each;
1203 const std::vector<double> &weights = quadrature.get_weights();
1204
1205 // Multiply quadrature weights by absolute value of Jacobian determinants or
1206 // the area element g=sqrt(DX^t DX) in case of codim > 0
1207
1208 if (update_flags & (update_normal_vectors | update_JxW_values))
1209 {
1210 AssertDimension(output_data.JxW_values.size(), n_q_points);
1211
1212 Assert(!(update_flags & update_normal_vectors) ||
1213 (output_data.normal_vectors.size() == n_q_points),
1214 ExcDimensionMismatch(output_data.normal_vectors.size(),
1215 n_q_points));
1216
1217
1218 if (computed_cell_similarity != CellSimilarity::translation)
1219 for (unsigned int point = 0; point < n_q_points; ++point)
1220 {
1221 if (dim == spacedim)
1222 {
1223 const double det = data.contravariant[point].determinant();
1224
1225 // check for distorted cells.
1226
1227 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1228 // 1e12 in 2d. might want to find a finer
1229 // (dimension-independent) criterion
1230 Assert(det >
1231 1e-12 * Utilities::fixed_power<dim>(
1232 cell->diameter() / std::sqrt(double(dim))),
1234 cell->center(), det, point)));
1235
1236 output_data.JxW_values[point] = weights[point] * det;
1237 }
1238 // if dim==spacedim, then there is no cell normal to
1239 // compute. since this is for FEValues (and not FEFaceValues),
1240 // there are also no face normals to compute
1241 else // codim>0 case
1242 {
1243 Tensor<1, spacedim> DX_t[dim];
1244 for (unsigned int i = 0; i < spacedim; ++i)
1245 for (unsigned int j = 0; j < dim; ++j)
1246 DX_t[j][i] = data.contravariant[point][i][j];
1247
1248 Tensor<2, dim> G; // First fundamental form
1249 for (unsigned int i = 0; i < dim; ++i)
1250 for (unsigned int j = 0; j < dim; ++j)
1251 G[i][j] = DX_t[i] * DX_t[j];
1252
1253 output_data.JxW_values[point] =
1254 std::sqrt(determinant(G)) * weights[point];
1255
1256 if (computed_cell_similarity ==
1258 {
1259 // we only need to flip the normal
1260 if (update_flags & update_normal_vectors)
1261 output_data.normal_vectors[point] *= -1.;
1262 }
1263 else
1264 {
1265 if (update_flags & update_normal_vectors)
1266 {
1267 Assert(spacedim == dim + 1,
1268 ExcMessage(
1269 "There is no (unique) cell normal for " +
1271 "-dimensional cells in " +
1272 Utilities::int_to_string(spacedim) +
1273 "-dimensional space. This only works if the "
1274 "space dimension is one greater than the "
1275 "dimensionality of the mesh cells."));
1276
1277 if (dim == 1)
1278 output_data.normal_vectors[point] =
1279 cross_product_2d(-DX_t[0]);
1280 else // dim == 2
1281 output_data.normal_vectors[point] =
1282 cross_product_3d(DX_t[0], DX_t[1]);
1283
1284 output_data.normal_vectors[point] /=
1285 output_data.normal_vectors[point].norm();
1286
1287 if (cell->direction_flag() == false)
1288 output_data.normal_vectors[point] *= -1.;
1289 }
1290 }
1291 } // codim>0 case
1292 }
1293 }
1294
1295
1296
1297 // copy values from InternalData to vector given by reference
1298 if (update_flags & update_jacobians)
1299 {
1300 AssertDimension(output_data.jacobians.size(), n_q_points);
1301 if (computed_cell_similarity != CellSimilarity::translation)
1302 for (unsigned int point = 0; point < n_q_points; ++point)
1303 output_data.jacobians[point] = data.contravariant[point];
1304 }
1305
1306 // copy values from InternalData to vector given by reference
1307 if (update_flags & update_inverse_jacobians)
1308 {
1309 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1310 if (computed_cell_similarity != CellSimilarity::translation)
1311 for (unsigned int point = 0; point < n_q_points; ++point)
1312 output_data.inverse_jacobians[point] =
1313 data.covariant[point].transpose();
1314 }
1315
1316 return computed_cell_similarity;
1317}
1318
1319
1320
1321namespace internal
1322{
1323 namespace MappingFEImplementation
1324 {
1325 namespace
1326 {
1337 template <int dim, int spacedim>
1338 void
1339 maybe_compute_face_data(
1340 const ::MappingFE<dim, spacedim> &mapping,
1341 const typename ::Triangulation<dim, spacedim>::cell_iterator
1342 &cell,
1343 const unsigned int face_no,
1344 const unsigned int subface_no,
1345 const unsigned int n_q_points,
1346 const typename QProjector<dim>::DataSetDescriptor data_set,
1347 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1349 &output_data)
1350 {
1351 const UpdateFlags update_flags = data.update_each;
1352
1353 if (update_flags &
1356 {
1357 if (update_flags & update_boundary_forms)
1358 AssertIndexRange(n_q_points,
1359 output_data.boundary_forms.size() + 1);
1360 if (update_flags & update_normal_vectors)
1361 AssertIndexRange(n_q_points,
1362 output_data.normal_vectors.size() + 1);
1363 if (update_flags & update_JxW_values)
1364 AssertIndexRange(n_q_points, output_data.JxW_values.size() + 1);
1365
1366 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1367
1368 // first compute some common data that is used for evaluating
1369 // all of the flags below
1370
1371 // map the unit tangentials to the real cell. checking for
1372 // d!=dim-1 eliminates compiler warnings regarding unsigned int
1373 // expressions < 0.
1374 for (unsigned int d = 0; d != dim - 1; ++d)
1375 {
1376 Assert(face_no + cell->n_faces() * d <
1377 data.unit_tangentials.size(),
1379 Assert(
1380 data.aux[d].size() <=
1381 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1383
1384 mapping.transform(
1386 data.unit_tangentials[face_no + cell->n_faces() * d]),
1388 data,
1389 make_array_view(data.aux[d]));
1390 }
1391
1392 if (update_flags & update_boundary_forms)
1393 {
1394 // if dim==spacedim, we can use the unit tangentials to
1395 // compute the boundary form by simply taking the cross
1396 // product
1397 if (dim == spacedim)
1398 {
1399 for (unsigned int i = 0; i < n_q_points; ++i)
1400 switch (dim)
1401 {
1402 case 1:
1403 // in 1d, we don't have access to any of the
1404 // data.aux fields (because it has only dim-1
1405 // components), but we can still compute the
1406 // boundary form by simply looking at the number
1407 // of the face
1408 output_data.boundary_forms[i][0] =
1409 (face_no == 0 ? -1 : +1);
1410 break;
1411 case 2:
1412 output_data.boundary_forms[i] =
1413 cross_product_2d(data.aux[0][i]);
1414 break;
1415 case 3:
1416 output_data.boundary_forms[i] =
1417 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1418 break;
1419 default:
1421 }
1422 }
1423 else //(dim < spacedim)
1424 {
1425 // in the codim-one case, the boundary form results from
1426 // the cross product of all the face tangential vectors
1427 // and the cell normal vector
1428 //
1429 // to compute the cell normal, use the same method used in
1430 // fill_fe_values for cells above
1431 AssertIndexRange(n_q_points, data.contravariant.size() + 1);
1432
1433 for (unsigned int point = 0; point < n_q_points; ++point)
1434 {
1435 if (dim == 1)
1436 {
1437 // J is a tangent vector
1438 output_data.boundary_forms[point] =
1439 data.contravariant[point].transpose()[0];
1440 output_data.boundary_forms[point] /=
1441 (face_no == 0 ? -1. : +1.) *
1442 output_data.boundary_forms[point].norm();
1443 }
1444
1445 if (dim == 2)
1446 {
1448 data.contravariant[point].transpose();
1449
1450 Tensor<1, spacedim> cell_normal =
1451 cross_product_3d(DX_t[0], DX_t[1]);
1452 cell_normal /= cell_normal.norm();
1453
1454 // then compute the face normal from the face
1455 // tangent and the cell normal:
1456 output_data.boundary_forms[point] =
1457 cross_product_3d(data.aux[0][point], cell_normal);
1458 }
1459 }
1460 }
1461 }
1462
1463 if (update_flags & update_JxW_values)
1464 for (unsigned int i = 0; i < n_q_points; ++i)
1465 {
1466 output_data.JxW_values[i] =
1467 output_data.boundary_forms[i].norm() *
1468 data.quadrature_weights[i + data_set];
1469
1470 if (subface_no != numbers::invalid_unsigned_int)
1471 {
1472 if (dim == 2)
1473 {
1474 const double area_ratio =
1475 1. / cell->reference_cell()
1476 .face_reference_cell(face_no)
1477 .n_isotropic_children();
1478 output_data.JxW_values[i] *= area_ratio;
1479 }
1480 else
1482 }
1483 }
1484
1485 if (update_flags & update_normal_vectors)
1486 for (unsigned int i = 0; i < n_q_points; ++i)
1487 output_data.normal_vectors[i] =
1488 Point<spacedim>(output_data.boundary_forms[i] /
1489 output_data.boundary_forms[i].norm());
1490
1491 if (update_flags & update_jacobians)
1492 for (unsigned int point = 0; point < n_q_points; ++point)
1493 output_data.jacobians[point] = data.contravariant[point];
1494
1495 if (update_flags & update_inverse_jacobians)
1496 for (unsigned int point = 0; point < n_q_points; ++point)
1497 output_data.inverse_jacobians[point] =
1498 data.covariant[point].transpose();
1499 }
1500 }
1501
1502
1509 template <int dim, int spacedim>
1510 void
1512 const ::MappingFE<dim, spacedim> &mapping,
1513 const typename ::Triangulation<dim, spacedim>::cell_iterator
1514 &cell,
1515 const unsigned int face_no,
1516 const unsigned int subface_no,
1517 const typename QProjector<dim>::DataSetDescriptor data_set,
1518 const Quadrature<dim - 1> &quadrature,
1519 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1521 &output_data)
1522 {
1523 const unsigned int n_q_points = quadrature.size();
1524
1525 maybe_compute_q_points<dim, spacedim>(data_set,
1526 data,
1527 output_data.quadrature_points,
1528 n_q_points);
1529 maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
1530 data_set,
1531 data,
1532 n_q_points);
1533 maybe_update_jacobian_grads<dim, spacedim>(CellSimilarity::none,
1534 data_set,
1535 data,
1536 output_data.jacobian_grads,
1537 n_q_points);
1538 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1540 data_set,
1541 data,
1543 n_q_points);
1544 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1546 data_set,
1547 data,
1548 output_data.jacobian_2nd_derivatives,
1549 n_q_points);
1550 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1552 data_set,
1553 data,
1555 n_q_points);
1556 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1558 data_set,
1559 data,
1560 output_data.jacobian_3rd_derivatives,
1561 n_q_points);
1562 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1564 data_set,
1565 data,
1567 n_q_points);
1568
1570 cell,
1571 face_no,
1572 subface_no,
1573 n_q_points,
1574 data_set,
1575 data,
1576 output_data);
1577 }
1578 } // namespace
1579 } // namespace MappingFEImplementation
1580} // namespace internal
1581
1582
1583
1584template <int dim, int spacedim>
1585void
1588 const unsigned int face_no,
1589 const hp::QCollection<dim - 1> &quadrature,
1590 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1592 &output_data) const
1593{
1594 // ensure that the following cast is really correct:
1595 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1597 const InternalData &data = static_cast<const InternalData &>(internal_data);
1598
1599 // if necessary, recompute the support points of the transformation of this
1600 // cell (note that we need to first check the triangulation pointer, since
1601 // otherwise the second test might trigger an exception if the
1602 // triangulations are not the same)
1603 if ((data.mapping_support_points.empty()) ||
1604 (&cell->get_triangulation() !=
1605 &data.cell_of_current_support_points->get_triangulation()) ||
1606 (cell != data.cell_of_current_support_points))
1607 {
1608 data.mapping_support_points = this->compute_mapping_support_points(cell);
1609 data.cell_of_current_support_points = cell;
1610 }
1611
1612 internal::MappingFEImplementation::do_fill_fe_face_values(
1613 *this,
1614 cell,
1615 face_no,
1617 QProjector<dim>::DataSetDescriptor::face(this->fe->reference_cell(),
1618 face_no,
1619 cell->combined_face_orientation(
1620 face_no),
1621 quadrature),
1622 quadrature[quadrature.size() == 1 ? 0 : face_no],
1623 data,
1624 output_data);
1625}
1626
1627
1628
1629template <int dim, int spacedim>
1630void
1633 const unsigned int face_no,
1634 const unsigned int subface_no,
1635 const Quadrature<dim - 1> &quadrature,
1636 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1638 &output_data) const
1639{
1640 // ensure that the following cast is really correct:
1641 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1643 const InternalData &data = static_cast<const InternalData &>(internal_data);
1644
1645 // if necessary, recompute the support points of the transformation of this
1646 // cell (note that we need to first check the triangulation pointer, since
1647 // otherwise the second test might trigger an exception if the
1648 // triangulations are not the same)
1649 if ((data.mapping_support_points.empty()) ||
1650 (&cell->get_triangulation() !=
1651 &data.cell_of_current_support_points->get_triangulation()) ||
1652 (cell != data.cell_of_current_support_points))
1653 {
1654 data.mapping_support_points = this->compute_mapping_support_points(cell);
1655 data.cell_of_current_support_points = cell;
1656 }
1657
1658 internal::MappingFEImplementation::do_fill_fe_face_values(
1659 *this,
1660 cell,
1661 face_no,
1662 subface_no,
1663 QProjector<dim>::DataSetDescriptor::subface(this->fe->reference_cell(),
1664 face_no,
1665 subface_no,
1666 cell->combined_face_orientation(
1667 face_no),
1668 quadrature.size(),
1669 cell->subface_case(face_no)),
1670 quadrature,
1671 data,
1672 output_data);
1673}
1674
1675
1676
1677namespace internal
1678{
1679 namespace MappingFEImplementation
1680 {
1681 namespace
1682 {
1683 template <int dim, int spacedim, int rank>
1684 void
1685 transform_fields(
1686 const ArrayView<const Tensor<rank, dim>> &input,
1687 const MappingKind mapping_kind,
1688 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1689 const ArrayView<Tensor<rank, spacedim>> &output)
1690 {
1691 // In the case of wedges and pyramids, faces might have different
1692 // numbers of quadrature points on each face with the result
1693 // that input and output have different sizes, since input has
1694 // the correct size but the size of output is the maximum of
1695 // all possible sizes.
1696 AssertIndexRange(input.size(), output.size() + 1);
1697
1698 Assert(
1699 (dynamic_cast<
1700 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1701 &mapping_data) != nullptr),
1703 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1704 static_cast<
1705 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1706 mapping_data);
1707
1708 switch (mapping_kind)
1709 {
1711 {
1712 Assert(
1715 "update_contravariant_transformation"));
1716
1717 for (unsigned int i = 0; i < input.size(); ++i)
1718 output[i] =
1719 apply_transformation(data.contravariant[i], input[i]);
1720
1721 return;
1722 }
1723
1724 case mapping_piola:
1725 {
1726 Assert(
1729 "update_contravariant_transformation"));
1730 Assert(
1731 data.update_each & update_volume_elements,
1733 "update_volume_elements"));
1734 Assert(rank == 1, ExcMessage("Only for rank 1"));
1735 if (rank != 1)
1736 return;
1737
1738 for (unsigned int i = 0; i < input.size(); ++i)
1739 {
1740 output[i] =
1741 apply_transformation(data.contravariant[i], input[i]);
1742 output[i] /= data.volume_elements[i];
1743 }
1744 return;
1745 }
1746 // We still allow this operation as in the
1747 // reference cell Derivatives are Tensor
1748 // rather than DerivativeForm
1749 case mapping_covariant:
1750 {
1751 Assert(
1754 "update_covariant_transformation"));
1755
1756 for (unsigned int i = 0; i < input.size(); ++i)
1757 output[i] = apply_transformation(data.covariant[i], input[i]);
1758
1759 return;
1760 }
1761
1762 default:
1764 }
1765 }
1766
1767
1768 template <int dim, int spacedim, int rank>
1769 void
1771 const ArrayView<const Tensor<rank, dim>> &input,
1772 const MappingKind mapping_kind,
1773 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1774 const ArrayView<Tensor<rank, spacedim>> &output)
1775 {
1776 AssertDimension(input.size(), output.size());
1777 Assert(
1778 (dynamic_cast<
1779 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1780 &mapping_data) != nullptr),
1782 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1783 static_cast<
1784 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1785 mapping_data);
1786
1787 switch (mapping_kind)
1788 {
1790 {
1791 Assert(
1794 "update_covariant_transformation"));
1795 Assert(
1798 "update_contravariant_transformation"));
1799 Assert(rank == 2, ExcMessage("Only for rank 2"));
1800
1801 for (unsigned int i = 0; i < output.size(); ++i)
1802 {
1804 apply_transformation(data.contravariant[i],
1805 transpose(input[i]));
1806 output[i] =
1807 apply_transformation(data.covariant[i], A.transpose());
1808 }
1809
1810 return;
1811 }
1812
1814 {
1815 Assert(
1818 "update_covariant_transformation"));
1819 Assert(rank == 2, ExcMessage("Only for rank 2"));
1820
1821 for (unsigned int i = 0; i < output.size(); ++i)
1822 {
1824 apply_transformation(data.covariant[i],
1825 transpose(input[i]));
1826 output[i] =
1827 apply_transformation(data.covariant[i], A.transpose());
1828 }
1829
1830 return;
1831 }
1832
1834 {
1835 Assert(
1838 "update_covariant_transformation"));
1839 Assert(
1842 "update_contravariant_transformation"));
1843 Assert(
1844 data.update_each & update_volume_elements,
1846 "update_volume_elements"));
1847 Assert(rank == 2, ExcMessage("Only for rank 2"));
1848
1849 for (unsigned int i = 0; i < output.size(); ++i)
1850 {
1852 apply_transformation(data.covariant[i], input[i]);
1853 const Tensor<2, spacedim> T =
1854 apply_transformation(data.contravariant[i],
1855 A.transpose());
1856
1857 output[i] = transpose(T);
1858 output[i] /= data.volume_elements[i];
1859 }
1860
1861 return;
1862 }
1863
1864 default:
1866 }
1867 }
1868
1869
1870
1871 template <int dim, int spacedim>
1872 void
1874 const ArrayView<const Tensor<3, dim>> &input,
1875 const MappingKind mapping_kind,
1876 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1877 const ArrayView<Tensor<3, spacedim>> &output)
1878 {
1879 AssertDimension(input.size(), output.size());
1880 Assert(
1881 (dynamic_cast<
1882 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1883 &mapping_data) != nullptr),
1885 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1886 static_cast<
1887 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1888 mapping_data);
1889
1890 switch (mapping_kind)
1891 {
1893 {
1894 Assert(
1897 "update_covariant_transformation"));
1898 Assert(
1901 "update_contravariant_transformation"));
1902
1903 for (unsigned int q = 0; q < output.size(); ++q)
1904 output[q] =
1906 data.contravariant[q],
1907 input[q]);
1908
1909 return;
1910 }
1911
1913 {
1914 Assert(
1917 "update_covariant_transformation"));
1918
1919 for (unsigned int q = 0; q < output.size(); ++q)
1920 output[q] =
1922 input[q]);
1923
1924 return;
1925 }
1926
1928 {
1929 Assert(
1932 "update_covariant_transformation"));
1933 Assert(
1936 "update_contravariant_transformation"));
1937 Assert(
1938 data.update_each & update_volume_elements,
1940 "update_volume_elements"));
1941
1942 for (unsigned int q = 0; q < output.size(); ++q)
1943 output[q] =
1945 data.contravariant[q],
1946 data.volume_elements[q],
1947 input[q]);
1948
1949 return;
1950 }
1951
1952 default:
1954 }
1955 }
1956
1957
1958
1959 template <int dim, int spacedim, int rank>
1960 void
1963 const MappingKind mapping_kind,
1964 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1966 {
1967 AssertDimension(input.size(), output.size());
1968 Assert(
1969 (dynamic_cast<
1970 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1971 &mapping_data) != nullptr),
1973 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1974 static_cast<
1975 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1976 mapping_data);
1977
1978 switch (mapping_kind)
1979 {
1980 case mapping_covariant:
1981 {
1982 Assert(
1985 "update_covariant_transformation"));
1986
1987 for (unsigned int i = 0; i < output.size(); ++i)
1988 output[i] = apply_transformation(data.covariant[i], input[i]);
1989
1990 return;
1991 }
1992 default:
1994 }
1995 }
1996 } // namespace
1997 } // namespace MappingFEImplementation
1998} // namespace internal
1999
2000
2001
2002template <int dim, int spacedim>
2003void
2005 const ArrayView<const Tensor<1, dim>> &input,
2006 const MappingKind mapping_kind,
2007 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2008 const ArrayView<Tensor<1, spacedim>> &output) const
2009{
2010 internal::MappingFEImplementation::transform_fields(input,
2011 mapping_kind,
2012 mapping_data,
2013 output);
2014}
2015
2016
2017
2018template <int dim, int spacedim>
2019void
2021 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2022 const MappingKind mapping_kind,
2023 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2024 const ArrayView<Tensor<2, spacedim>> &output) const
2025{
2026 internal::MappingFEImplementation::transform_differential_forms(input,
2027 mapping_kind,
2028 mapping_data,
2029 output);
2030}
2031
2032
2033
2034template <int dim, int spacedim>
2035void
2037 const ArrayView<const Tensor<2, dim>> &input,
2038 const MappingKind mapping_kind,
2039 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2040 const ArrayView<Tensor<2, spacedim>> &output) const
2041{
2042 switch (mapping_kind)
2043 {
2045 internal::MappingFEImplementation::transform_fields(input,
2046 mapping_kind,
2047 mapping_data,
2048 output);
2049 return;
2050
2054 internal::MappingFEImplementation::transform_gradients(input,
2055 mapping_kind,
2056 mapping_data,
2057 output);
2058 return;
2059 default:
2061 }
2062}
2063
2064
2065
2066template <int dim, int spacedim>
2067void
2069 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2070 const MappingKind mapping_kind,
2071 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2072 const ArrayView<Tensor<3, spacedim>> &output) const
2073{
2074 AssertDimension(input.size(), output.size());
2075 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2077 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2078
2079 switch (mapping_kind)
2080 {
2082 {
2085 "update_covariant_transformation"));
2086
2087 for (unsigned int q = 0; q < output.size(); ++q)
2088 for (unsigned int i = 0; i < spacedim; ++i)
2089 for (unsigned int j = 0; j < spacedim; ++j)
2090 {
2091 double tmp[dim];
2092 for (unsigned int K = 0; K < dim; ++K)
2093 {
2094 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
2095 for (unsigned int J = 1; J < dim; ++J)
2096 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
2097 }
2098 for (unsigned int k = 0; k < spacedim; ++k)
2099 {
2100 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
2101 for (unsigned int K = 1; K < dim; ++K)
2102 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
2103 }
2104 }
2105 return;
2106 }
2107
2108 default:
2110 }
2111}
2112
2113
2114
2115template <int dim, int spacedim>
2116void
2118 const ArrayView<const Tensor<3, dim>> &input,
2119 const MappingKind mapping_kind,
2120 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2121 const ArrayView<Tensor<3, spacedim>> &output) const
2122{
2123 switch (mapping_kind)
2124 {
2128 internal::MappingFEImplementation::transform_hessians(input,
2129 mapping_kind,
2130 mapping_data,
2131 output);
2132 return;
2133 default:
2135 }
2136}
2137
2138
2139
2140namespace
2141{
2142 template <int spacedim>
2143 bool
2144 check_all_manifold_ids_identical(
2146 {
2147 return true;
2148 }
2149
2150
2151
2152 template <int spacedim>
2153 bool
2154 check_all_manifold_ids_identical(
2156 {
2157 const auto m_id = cell->manifold_id();
2158
2159 for (const auto f : cell->face_indices())
2160 if (m_id != cell->face(f)->manifold_id())
2161 return false;
2162
2163 return true;
2164 }
2165
2166
2167
2168 template <int spacedim>
2169 bool
2170 check_all_manifold_ids_identical(
2172 {
2173 const auto m_id = cell->manifold_id();
2174
2175 for (const auto f : cell->face_indices())
2176 if (m_id != cell->face(f)->manifold_id())
2177 return false;
2178
2179 for (const auto l : cell->line_indices())
2180 if (m_id != cell->line(l)->manifold_id())
2181 return false;
2182
2183 return true;
2184 }
2185} // namespace
2186
2187
2188
2189template <int dim, int spacedim>
2190std::vector<Point<spacedim>>
2192 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2193{
2194 std::vector<Point<spacedim>> mapping_support_points(
2195 fe->get_unit_support_points().size());
2196
2197 std::vector<Point<spacedim>> vertices(cell->n_vertices());
2198 for (const unsigned int i : cell->vertex_indices())
2199 vertices[i] = cell->vertex(i);
2200
2201 if (check_all_manifold_ids_identical(cell))
2202 {
2203 // version 1) all geometric entities have the same manifold:
2204
2205 cell->get_manifold().get_new_points(vertices,
2206 mapping_support_point_weights,
2207 mapping_support_points);
2208 }
2209 else
2210 {
2211 // version 1) geometric entities have different manifold
2212
2213 // helper function to compute mapped points on subentities
2214 // note[PM]: this function currently only uses the vertices
2215 // of cells to create new points on subentities; however,
2216 // one should use all bounding points to create new points
2217 // as in the case of MappingQ.
2218 const auto process = [&](const auto &manifold,
2219 const auto &indices,
2220 const unsigned n_points) {
2221 if ((indices.size() == 0) || (n_points == 0))
2222 return;
2223
2224 const unsigned int n_shape_functions =
2225 this->fe->reference_cell().n_vertices();
2226
2227 Table<2, double> mapping_support_point_weights_local(n_points,
2228 n_shape_functions);
2229 std::vector<Point<spacedim>> mapping_support_points_local(n_points);
2230
2231 for (unsigned int p = 0; p < n_points; ++p)
2232 for (unsigned int i = 0; i < n_shape_functions; ++i)
2233 mapping_support_point_weights_local(p, i) =
2234 mapping_support_point_weights(
2235 indices[p + (indices.size() - n_points)], i);
2236
2237 manifold.get_new_points(vertices,
2238 mapping_support_point_weights_local,
2239 mapping_support_points_local);
2240
2241 for (unsigned int p = 0; p < n_points; ++p)
2242 mapping_support_points[indices[p + (indices.size() - n_points)]] =
2243 mapping_support_points_local[p];
2244 };
2245
2246 // create dummy DoFHandler to extract indices on subobjects
2247 const auto &fe = *this->fe;
2249 GridGenerator::reference_cell(tria, fe.reference_cell());
2250 DoFHandler<dim, spacedim> dof_handler(tria);
2251 dof_handler.distribute_dofs(fe);
2252 const auto &cell_ref = dof_handler.begin_active();
2253
2254 std::vector<types::global_dof_index> indices;
2255
2256 // add vertices
2257 for (const unsigned int i : cell->vertex_indices())
2258 mapping_support_points[i] = cell->vertex(i);
2259
2260 // process and add line support points
2261 for (unsigned int l = 0; l < cell_ref->n_lines(); ++l)
2262 {
2263 const auto accessor = cell_ref->line(l);
2264 indices.resize(fe.n_dofs_per_line() + 2 * fe.n_dofs_per_vertex());
2265 accessor->get_dof_indices(indices);
2266 process(cell->line(l)->get_manifold(), indices, fe.n_dofs_per_line());
2267 }
2268
2269 // process and add face support points
2270 if constexpr (dim >= 3)
2271 {
2272 for (unsigned int f = 0; f < cell_ref->n_faces(); ++f)
2273 {
2274 const auto accessor = cell_ref->face(f);
2275 indices.resize(fe.n_dofs_per_face());
2276 accessor->get_dof_indices(indices);
2277 process(cell->face(f)->get_manifold(),
2278 indices,
2279 fe.n_dofs_per_quad());
2280 }
2281 }
2282
2283 // process and add volume support points
2284 if constexpr (dim >= 2)
2285 {
2286 indices.resize(fe.n_dofs_per_cell());
2287 cell_ref->get_dof_indices(indices);
2288 process(cell->get_manifold(),
2289 indices,
2290 (dim == 2) ? fe.n_dofs_per_quad() : fe.n_dofs_per_hex());
2291 }
2292 }
2293
2294 return mapping_support_points;
2295}
2296
2297
2298
2299template <int dim, int spacedim>
2302 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2303{
2304 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
2305}
2306
2307
2308
2309template <int dim, int spacedim>
2310bool
2312 const ReferenceCell &reference_cell) const
2313{
2314 Assert(dim == reference_cell.get_dimension(),
2315 ExcMessage("The dimension of your mapping (" +
2317 ") and the reference cell cell_type (" +
2318 Utilities::to_string(reference_cell.get_dimension()) +
2319 " ) do not agree."));
2320
2321 return fe->reference_cell() == reference_cell;
2322}
2323
2324
2325
2326//--------------------------- Explicit instantiations -----------------------
2327#include "fe/mapping_fe.inst"
2328
2329
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:61
InternalData(const FiniteElement< dim, spacedim > &fe)
Definition mapping_fe.cc:50
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
Definition mapping_fe.cc:82
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:750
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)
Tensor< 3, spacedim, Number > apply_contravariant_hessian(const DerivativeForm< 1, dim, spacedim, Number > &covariant, const DerivativeForm< 1, dim, spacedim, Number > &contravariant, const Tensor< 3, dim, Number > &input)
Tensor< 3, spacedim, Number > apply_piola_hessian(const DerivativeForm< 1, dim, spacedim, Number > &covariant, const DerivativeForm< 1, dim, spacedim, Number > &contravariant, const Number &volume_element, const Tensor< 3, dim, Number > &input)
Tensor< 3, spacedim, Number > apply_covariant_hessian(const DerivativeForm< 1, dim, spacedim, Number > &covariant, const Tensor< 3, dim, Number > &input)
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 > &)