Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+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
30#include <deal.II/grid/tria.h>
32
34
35#include <boost/container/small_vector.hpp>
36
37#include <algorithm>
38#include <array>
39#include <cmath>
40#include <memory>
41#include <numeric>
42
43
45
46
47template <int dim, int spacedim>
50 : fe(fe)
51 , polynomial_degree(fe.tensor_degree())
52 , n_shape_functions(fe.n_dofs_per_cell())
53{}
54
55
56
57template <int dim, int spacedim>
58std::size_t
75
76
77template <int dim, int spacedim>
78void
80 const UpdateFlags update_flags,
81 const Quadrature<dim> &q,
82 const unsigned int n_original_q_points)
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_original_q_points);
92
93 if (this->update_each & update_contravariant_transformation)
94 contravariant.resize(n_original_q_points);
95
96 if (this->update_each & update_volume_elements)
97 volume_elements.resize(n_original_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 initialize(update_flags, q, n_original_q_points);
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 unit_tangential_vectors<dim>(i, 0));
161 if (dim > 2)
162 {
163 unit_tangentials[n_faces + i].resize(n_original_q_points);
164 std::fill(
165 unit_tangentials[n_faces + i].begin(),
166 unit_tangentials[n_faces + i].end(),
167 reference_cell.template unit_tangential_vectors<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 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1063 data.initialize(this->requires_update_flags(update_flags), q, q.size());
1064
1065 return data_ptr;
1066}
1067
1068
1069
1070template <int dim, int spacedim>
1071std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1073 const UpdateFlags update_flags,
1074 const hp::QCollection<dim - 1> &quadrature) const
1075{
1076 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1077 std::make_unique<InternalData>(*this->fe);
1078 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1079 data.initialize_face(this->requires_update_flags(update_flags),
1081 this->fe->reference_cell(), quadrature),
1082 quadrature.max_n_quadrature_points());
1083
1084 return data_ptr;
1085}
1086
1087
1088
1089template <int dim, int spacedim>
1090std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1092 const UpdateFlags update_flags,
1093 const Quadrature<dim - 1> &quadrature) const
1094{
1095 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1096 std::make_unique<InternalData>(*this->fe);
1097 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1098 data.initialize_face(this->requires_update_flags(update_flags),
1100 this->fe->reference_cell(), quadrature),
1101 quadrature.size());
1102
1103 return data_ptr;
1104}
1105
1106
1107
1108template <int dim, int spacedim>
1112 const CellSimilarity::Similarity cell_similarity,
1113 const Quadrature<dim> &quadrature,
1114 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1116 &output_data) const
1117{
1118 // ensure that the following static_cast is really correct:
1119 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1121 const InternalData &data = static_cast<const InternalData &>(internal_data);
1122
1123 const unsigned int n_q_points = quadrature.size();
1124
1125 // recompute the support points of the transformation of this
1126 // cell. we tried to be clever here in an earlier version of the
1127 // library by checking whether the cell is the same as the one we
1128 // had visited last, but it turns out to be difficult to determine
1129 // that because a cell for the purposes of a mapping is
1130 // characterized not just by its (triangulation, level, index)
1131 // triple, but also by the locations of its vertices, the manifold
1132 // object attached to the cell and all of its bounding faces/edges,
1133 // etc. to reliably test that the "cell" we are on is, therefore,
1134 // not easily done
1135 data.mapping_support_points = this->compute_mapping_support_points(cell);
1137
1138 // if the order of the mapping is greater than 1, then do not reuse any cell
1139 // similarity information. This is necessary because the cell similarity
1140 // value is computed with just cell vertices and does not take into account
1141 // cell curvature.
1142 const CellSimilarity::Similarity computed_cell_similarity =
1143 (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
1144
1145 internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1147 data,
1148 output_data.quadrature_points,
1149 n_q_points);
1150
1151 internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1152 computed_cell_similarity,
1154 data,
1155 n_q_points);
1156
1157 internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1158 computed_cell_similarity,
1160 data,
1161 output_data.jacobian_grads,
1162 n_q_points);
1163
1164 internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
1165 dim,
1166 spacedim>(computed_cell_similarity,
1168 data,
1170 n_q_points);
1171
1172 internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
1173 dim,
1174 spacedim>(computed_cell_similarity,
1176 data,
1177 output_data.jacobian_2nd_derivatives,
1178 n_q_points);
1179
1180 internal::MappingFEImplementation::
1181 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1182 computed_cell_similarity,
1184 data,
1186 n_q_points);
1187
1188 internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
1189 dim,
1190 spacedim>(computed_cell_similarity,
1192 data,
1193 output_data.jacobian_3rd_derivatives,
1194 n_q_points);
1195
1196 internal::MappingFEImplementation::
1197 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1198 computed_cell_similarity,
1200 data,
1202 n_q_points);
1203
1204 const UpdateFlags update_flags = data.update_each;
1205 const std::vector<double> &weights = quadrature.get_weights();
1206
1207 // Multiply quadrature weights by absolute value of Jacobian determinants or
1208 // the area element g=sqrt(DX^t DX) in case of codim > 0
1209
1210 if (update_flags & (update_normal_vectors | update_JxW_values))
1211 {
1212 AssertDimension(output_data.JxW_values.size(), n_q_points);
1213
1214 Assert(!(update_flags & update_normal_vectors) ||
1215 (output_data.normal_vectors.size() == n_q_points),
1216 ExcDimensionMismatch(output_data.normal_vectors.size(),
1217 n_q_points));
1218
1219
1220 if (computed_cell_similarity != CellSimilarity::translation)
1221 for (unsigned int point = 0; point < n_q_points; ++point)
1222 {
1223 if (dim == spacedim)
1224 {
1225 const double det = data.contravariant[point].determinant();
1226
1227 // check for distorted cells.
1228
1229 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1230 // 1e12 in 2d. might want to find a finer
1231 // (dimension-independent) criterion
1232 Assert(det >
1233 1e-12 * Utilities::fixed_power<dim>(
1234 cell->diameter() / std::sqrt(double(dim))),
1236 cell->center(), det, point)));
1237
1238 output_data.JxW_values[point] = weights[point] * det;
1239 }
1240 // if dim==spacedim, then there is no cell normal to
1241 // compute. since this is for FEValues (and not FEFaceValues),
1242 // there are also no face normals to compute
1243 else // codim>0 case
1244 {
1245 Tensor<1, spacedim> DX_t[dim];
1246 for (unsigned int i = 0; i < spacedim; ++i)
1247 for (unsigned int j = 0; j < dim; ++j)
1248 DX_t[j][i] = data.contravariant[point][i][j];
1249
1250 Tensor<2, dim> G; // First fundamental form
1251 for (unsigned int i = 0; i < dim; ++i)
1252 for (unsigned int j = 0; j < dim; ++j)
1253 G[i][j] = DX_t[i] * DX_t[j];
1254
1255 output_data.JxW_values[point] =
1256 std::sqrt(determinant(G)) * weights[point];
1257
1258 if (computed_cell_similarity ==
1260 {
1261 // we only need to flip the normal
1262 if (update_flags & update_normal_vectors)
1263 output_data.normal_vectors[point] *= -1.;
1264 }
1265 else
1266 {
1267 if (update_flags & update_normal_vectors)
1268 {
1269 Assert(spacedim == dim + 1,
1270 ExcMessage(
1271 "There is no (unique) cell normal for " +
1273 "-dimensional cells in " +
1274 Utilities::int_to_string(spacedim) +
1275 "-dimensional space. This only works if the "
1276 "space dimension is one greater than the "
1277 "dimensionality of the mesh cells."));
1278
1279 if (dim == 1)
1280 output_data.normal_vectors[point] =
1281 cross_product_2d(-DX_t[0]);
1282 else // dim == 2
1283 output_data.normal_vectors[point] =
1284 cross_product_3d(DX_t[0], DX_t[1]);
1285
1286 output_data.normal_vectors[point] /=
1287 output_data.normal_vectors[point].norm();
1288
1289 if (cell->direction_flag() == false)
1290 output_data.normal_vectors[point] *= -1.;
1291 }
1292 }
1293 } // codim>0 case
1294 }
1295 }
1296
1297
1298
1299 // copy values from InternalData to vector given by reference
1300 if (update_flags & update_jacobians)
1301 {
1302 AssertDimension(output_data.jacobians.size(), n_q_points);
1303 if (computed_cell_similarity != CellSimilarity::translation)
1304 for (unsigned int point = 0; point < n_q_points; ++point)
1305 output_data.jacobians[point] = data.contravariant[point];
1306 }
1307
1308 // copy values from InternalData to vector given by reference
1309 if (update_flags & update_inverse_jacobians)
1310 {
1311 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1312 if (computed_cell_similarity != CellSimilarity::translation)
1313 for (unsigned int point = 0; point < n_q_points; ++point)
1314 output_data.inverse_jacobians[point] =
1315 data.covariant[point].transpose();
1316 }
1317
1318 return computed_cell_similarity;
1319}
1320
1321
1322
1323namespace internal
1324{
1325 namespace MappingFEImplementation
1326 {
1327 namespace
1328 {
1339 template <int dim, int spacedim>
1340 void
1341 maybe_compute_face_data(
1342 const ::MappingFE<dim, spacedim> &mapping,
1343 const typename ::Triangulation<dim, spacedim>::cell_iterator
1344 &cell,
1345 const unsigned int face_no,
1346 const unsigned int subface_no,
1347 const unsigned int n_q_points,
1348 const typename QProjector<dim>::DataSetDescriptor data_set,
1349 const typename ::MappingFE<dim, spacedim>::InternalData &data,
1351 &output_data)
1352 {
1353 const UpdateFlags update_flags = data.update_each;
1354
1355 if (update_flags &
1358 {
1359 if (update_flags & update_boundary_forms)
1360 AssertIndexRange(n_q_points,
1361 output_data.boundary_forms.size() + 1);
1362 if (update_flags & update_normal_vectors)
1363 AssertIndexRange(n_q_points,
1364 output_data.normal_vectors.size() + 1);
1365 if (update_flags & update_JxW_values)
1366 AssertIndexRange(n_q_points, output_data.JxW_values.size() + 1);
1367
1368 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1369
1370 // first compute some common data that is used for evaluating
1371 // all of the flags below
1372
1373 // map the unit tangentials to the real cell. checking for
1374 // d!=dim-1 eliminates compiler warnings regarding unsigned int
1375 // expressions < 0.
1376 for (unsigned int d = 0; d != dim - 1; ++d)
1377 {
1378 Assert(face_no + cell->n_faces() * d <
1379 data.unit_tangentials.size(),
1381 Assert(
1382 data.aux[d].size() <=
1383 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1385
1386 mapping.transform(
1388 data.unit_tangentials[face_no + cell->n_faces() * d]),
1390 data,
1391 make_array_view(data.aux[d]));
1392 }
1393
1394 if (update_flags & update_boundary_forms)
1395 {
1396 // if dim==spacedim, we can use the unit tangentials to
1397 // compute the boundary form by simply taking the cross
1398 // product
1399 if (dim == spacedim)
1400 {
1401 for (unsigned int i = 0; i < n_q_points; ++i)
1402 switch (dim)
1403 {
1404 case 1:
1405 // in 1d, we don't have access to any of the
1406 // data.aux fields (because it has only dim-1
1407 // components), but we can still compute the
1408 // boundary form by simply looking at the number
1409 // of the face
1410 output_data.boundary_forms[i][0] =
1411 (face_no == 0 ? -1 : +1);
1412 break;
1413 case 2:
1414 output_data.boundary_forms[i] =
1415 cross_product_2d(data.aux[0][i]);
1416 break;
1417 case 3:
1418 output_data.boundary_forms[i] =
1419 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1420 break;
1421 default:
1423 }
1424 }
1425 else //(dim < spacedim)
1426 {
1427 // in the codim-one case, the boundary form results from
1428 // the cross product of all the face tangential vectors
1429 // and the cell normal vector
1430 //
1431 // to compute the cell normal, use the same method used in
1432 // fill_fe_values for cells above
1433 AssertIndexRange(n_q_points, data.contravariant.size() + 1);
1434
1435 for (unsigned int point = 0; point < n_q_points; ++point)
1436 {
1437 if (dim == 1)
1438 {
1439 // J is a tangent vector
1440 output_data.boundary_forms[point] =
1441 data.contravariant[point].transpose()[0];
1442 output_data.boundary_forms[point] /=
1443 (face_no == 0 ? -1. : +1.) *
1444 output_data.boundary_forms[point].norm();
1445 }
1446
1447 if (dim == 2)
1448 {
1450 data.contravariant[point].transpose();
1451
1452 Tensor<1, spacedim> cell_normal =
1453 cross_product_3d(DX_t[0], DX_t[1]);
1454 cell_normal /= cell_normal.norm();
1455
1456 // then compute the face normal from the face
1457 // tangent and the cell normal:
1458 output_data.boundary_forms[point] =
1459 cross_product_3d(data.aux[0][point], cell_normal);
1460 }
1461 }
1462 }
1463 }
1464
1465 if (update_flags & update_JxW_values)
1466 for (unsigned int i = 0; i < n_q_points; ++i)
1467 {
1468 output_data.JxW_values[i] =
1469 output_data.boundary_forms[i].norm() *
1470 data.quadrature_weights[i + data_set];
1471
1472 if (subface_no != numbers::invalid_unsigned_int)
1473 {
1474#if false
1475 const double area_ratio =
1477 cell->subface_case(face_no), subface_no);
1478 output_data.JxW_values[i] *= area_ratio;
1479#else
1481#endif
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() !=
1606 (cell != data.cell_of_current_support_points))
1607 {
1608 data.mapping_support_points = this->compute_mapping_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->face_orientation(face_no),
1620 cell->face_flip(face_no),
1621 cell->face_rotation(face_no),
1622 quadrature),
1623 quadrature[quadrature.size() == 1 ? 0 : face_no],
1624 data,
1625 output_data);
1626}
1627
1628
1629
1630template <int dim, int spacedim>
1631void
1634 const unsigned int face_no,
1635 const unsigned int subface_no,
1636 const Quadrature<dim - 1> &quadrature,
1637 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1639 &output_data) const
1640{
1641 // ensure that the following cast is really correct:
1642 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1644 const InternalData &data = static_cast<const InternalData &>(internal_data);
1645
1646 // if necessary, recompute the support points of the transformation of this
1647 // cell (note that we need to first check the triangulation pointer, since
1648 // otherwise the second test might trigger an exception if the
1649 // triangulations are not the same)
1650 if ((data.mapping_support_points.empty()) ||
1651 (&cell->get_triangulation() !=
1653 (cell != data.cell_of_current_support_points))
1654 {
1655 data.mapping_support_points = this->compute_mapping_support_points(cell);
1657 }
1658
1659 internal::MappingFEImplementation::do_fill_fe_face_values(
1660 *this,
1661 cell,
1662 face_no,
1663 subface_no,
1664 QProjector<dim>::DataSetDescriptor::subface(this->fe->reference_cell(),
1665 face_no,
1666 subface_no,
1667 cell->face_orientation(face_no),
1668 cell->face_flip(face_no),
1669 cell->face_rotation(face_no),
1670 quadrature.size(),
1671 cell->subface_case(face_no)),
1672 quadrature,
1673 data,
1674 output_data);
1675}
1676
1677
1678
1679namespace internal
1680{
1681 namespace MappingFEImplementation
1682 {
1683 namespace
1684 {
1685 template <int dim, int spacedim, int rank>
1686 void
1687 transform_fields(
1688 const ArrayView<const Tensor<rank, dim>> &input,
1689 const MappingKind mapping_kind,
1690 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1691 const ArrayView<Tensor<rank, spacedim>> &output)
1692 {
1693 // In the case of wedges and pyramids, faces might have different
1694 // numbers of quadrature points on each face with the result
1695 // that input and output have different sizes, since input has
1696 // the correct size but the size of output is the maximum of
1697 // all possible sizes.
1698 AssertIndexRange(input.size(), output.size() + 1);
1699
1700 Assert(
1701 (dynamic_cast<
1702 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1703 &mapping_data) != nullptr),
1705 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1706 static_cast<
1707 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1708 mapping_data);
1709
1710 switch (mapping_kind)
1711 {
1713 {
1714 Assert(
1715 data.update_each & update_contravariant_transformation,
1717 "update_contravariant_transformation"));
1718
1719 for (unsigned int i = 0; i < input.size(); ++i)
1720 output[i] =
1721 apply_transformation(data.contravariant[i], input[i]);
1722
1723 return;
1724 }
1725
1726 case mapping_piola:
1727 {
1728 Assert(
1729 data.update_each & update_contravariant_transformation,
1731 "update_contravariant_transformation"));
1732 Assert(
1733 data.update_each & update_volume_elements,
1735 "update_volume_elements"));
1736 Assert(rank == 1, ExcMessage("Only for rank 1"));
1737 if (rank != 1)
1738 return;
1739
1740 for (unsigned int i = 0; i < input.size(); ++i)
1741 {
1742 output[i] =
1743 apply_transformation(data.contravariant[i], input[i]);
1744 output[i] /= data.volume_elements[i];
1745 }
1746 return;
1747 }
1748 // We still allow this operation as in the
1749 // reference cell Derivatives are Tensor
1750 // rather than DerivativeForm
1751 case mapping_covariant:
1752 {
1753 Assert(
1754 data.update_each & update_contravariant_transformation,
1756 "update_covariant_transformation"));
1757
1758 for (unsigned int i = 0; i < input.size(); ++i)
1759 output[i] = apply_transformation(data.covariant[i], input[i]);
1760
1761 return;
1762 }
1763
1764 default:
1766 }
1767 }
1768
1769
1770 template <int dim, int spacedim, int rank>
1771 void
1773 const ArrayView<const Tensor<rank, dim>> &input,
1774 const MappingKind mapping_kind,
1775 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1776 const ArrayView<Tensor<rank, spacedim>> &output)
1777 {
1778 AssertDimension(input.size(), output.size());
1779 Assert(
1780 (dynamic_cast<
1781 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1782 &mapping_data) != nullptr),
1784 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1785 static_cast<
1786 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1787 mapping_data);
1788
1789 switch (mapping_kind)
1790 {
1792 {
1793 Assert(
1794 data.update_each & update_covariant_transformation,
1796 "update_covariant_transformation"));
1797 Assert(
1798 data.update_each & update_contravariant_transformation,
1800 "update_contravariant_transformation"));
1801 Assert(rank == 2, ExcMessage("Only for rank 2"));
1802
1803 for (unsigned int i = 0; i < output.size(); ++i)
1804 {
1806 apply_transformation(data.contravariant[i],
1807 transpose(input[i]));
1808 output[i] =
1809 apply_transformation(data.covariant[i], A.transpose());
1810 }
1811
1812 return;
1813 }
1814
1816 {
1817 Assert(
1818 data.update_each & update_covariant_transformation,
1820 "update_covariant_transformation"));
1821 Assert(rank == 2, ExcMessage("Only for rank 2"));
1822
1823 for (unsigned int i = 0; i < output.size(); ++i)
1824 {
1826 apply_transformation(data.covariant[i],
1827 transpose(input[i]));
1828 output[i] =
1829 apply_transformation(data.covariant[i], A.transpose());
1830 }
1831
1832 return;
1833 }
1834
1836 {
1837 Assert(
1838 data.update_each & update_covariant_transformation,
1840 "update_covariant_transformation"));
1841 Assert(
1842 data.update_each & update_contravariant_transformation,
1844 "update_contravariant_transformation"));
1845 Assert(
1846 data.update_each & update_volume_elements,
1848 "update_volume_elements"));
1849 Assert(rank == 2, ExcMessage("Only for rank 2"));
1850
1851 for (unsigned int i = 0; i < output.size(); ++i)
1852 {
1854 apply_transformation(data.covariant[i], input[i]);
1855 const Tensor<2, spacedim> T =
1856 apply_transformation(data.contravariant[i],
1857 A.transpose());
1858
1859 output[i] = transpose(T);
1860 output[i] /= data.volume_elements[i];
1861 }
1862
1863 return;
1864 }
1865
1866 default:
1868 }
1869 }
1870
1871
1872
1873 template <int dim, int spacedim>
1874 void
1876 const ArrayView<const Tensor<3, dim>> &input,
1877 const MappingKind mapping_kind,
1878 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1879 const ArrayView<Tensor<3, spacedim>> &output)
1880 {
1881 AssertDimension(input.size(), output.size());
1882 Assert(
1883 (dynamic_cast<
1884 const typename ::MappingFE<dim, spacedim>::InternalData *>(
1885 &mapping_data) != nullptr),
1887 const typename ::MappingFE<dim, spacedim>::InternalData &data =
1888 static_cast<
1889 const typename ::MappingFE<dim, spacedim>::InternalData &>(
1890 mapping_data);
1891
1892 switch (mapping_kind)
1893 {
1895 {
1896 Assert(
1897 data.update_each & update_covariant_transformation,
1899 "update_covariant_transformation"));
1900 Assert(
1901 data.update_each & update_contravariant_transformation,
1903 "update_contravariant_transformation"));
1904
1905 for (unsigned int q = 0; q < output.size(); ++q)
1906 for (unsigned int i = 0; i < spacedim; ++i)
1907 {
1908 double tmp1[dim][dim];
1909 for (unsigned int J = 0; J < dim; ++J)
1910 for (unsigned int K = 0; K < dim; ++K)
1911 {
1912 tmp1[J][K] =
1913 data.contravariant[q][i][0] * input[q][0][J][K];
1914 for (unsigned int I = 1; I < dim; ++I)
1915 tmp1[J][K] +=
1916 data.contravariant[q][i][I] * input[q][I][J][K];
1917 }
1918 for (unsigned int j = 0; j < spacedim; ++j)
1919 {
1920 double tmp2[dim];
1921 for (unsigned int K = 0; K < dim; ++K)
1922 {
1923 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1924 for (unsigned int J = 1; J < dim; ++J)
1925 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1926 }
1927 for (unsigned int k = 0; k < spacedim; ++k)
1928 {
1929 output[q][i][j][k] =
1930 data.covariant[q][k][0] * tmp2[0];
1931 for (unsigned int K = 1; K < dim; ++K)
1932 output[q][i][j][k] +=
1933 data.covariant[q][k][K] * tmp2[K];
1934 }
1935 }
1936 }
1937 return;
1938 }
1939
1941 {
1942 Assert(
1943 data.update_each & update_covariant_transformation,
1945 "update_covariant_transformation"));
1946
1947 for (unsigned int q = 0; q < output.size(); ++q)
1948 for (unsigned int i = 0; i < spacedim; ++i)
1949 {
1950 double tmp1[dim][dim];
1951 for (unsigned int J = 0; J < dim; ++J)
1952 for (unsigned int K = 0; K < dim; ++K)
1953 {
1954 tmp1[J][K] =
1955 data.covariant[q][i][0] * input[q][0][J][K];
1956 for (unsigned int I = 1; I < dim; ++I)
1957 tmp1[J][K] +=
1958 data.covariant[q][i][I] * input[q][I][J][K];
1959 }
1960 for (unsigned int j = 0; j < spacedim; ++j)
1961 {
1962 double tmp2[dim];
1963 for (unsigned int K = 0; K < dim; ++K)
1964 {
1965 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1966 for (unsigned int J = 1; J < dim; ++J)
1967 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1968 }
1969 for (unsigned int k = 0; k < spacedim; ++k)
1970 {
1971 output[q][i][j][k] =
1972 data.covariant[q][k][0] * tmp2[0];
1973 for (unsigned int K = 1; K < dim; ++K)
1974 output[q][i][j][k] +=
1975 data.covariant[q][k][K] * tmp2[K];
1976 }
1977 }
1978 }
1979
1980 return;
1981 }
1982
1984 {
1985 Assert(
1986 data.update_each & update_covariant_transformation,
1988 "update_covariant_transformation"));
1989 Assert(
1990 data.update_each & update_contravariant_transformation,
1992 "update_contravariant_transformation"));
1993 Assert(
1994 data.update_each & update_volume_elements,
1996 "update_volume_elements"));
1997
1998 for (unsigned int q = 0; q < output.size(); ++q)
1999 for (unsigned int i = 0; i < spacedim; ++i)
2000 {
2001 double factor[dim];
2002 for (unsigned int I = 0; I < dim; ++I)
2003 factor[I] =
2004 data.contravariant[q][i][I] / data.volume_elements[q];
2005 double tmp1[dim][dim];
2006 for (unsigned int J = 0; J < dim; ++J)
2007 for (unsigned int K = 0; K < dim; ++K)
2008 {
2009 tmp1[J][K] = factor[0] * input[q][0][J][K];
2010 for (unsigned int I = 1; I < dim; ++I)
2011 tmp1[J][K] += factor[I] * input[q][I][J][K];
2012 }
2013 for (unsigned int j = 0; j < spacedim; ++j)
2014 {
2015 double tmp2[dim];
2016 for (unsigned int K = 0; K < dim; ++K)
2017 {
2018 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2019 for (unsigned int J = 1; J < dim; ++J)
2020 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2021 }
2022 for (unsigned int k = 0; k < spacedim; ++k)
2023 {
2024 output[q][i][j][k] =
2025 data.covariant[q][k][0] * tmp2[0];
2026 for (unsigned int K = 1; K < dim; ++K)
2027 output[q][i][j][k] +=
2028 data.covariant[q][k][K] * tmp2[K];
2029 }
2030 }
2031 }
2032
2033 return;
2034 }
2035
2036 default:
2038 }
2039 }
2040
2041
2042
2043 template <int dim, int spacedim, int rank>
2044 void
2047 const MappingKind mapping_kind,
2048 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2050 {
2051 AssertDimension(input.size(), output.size());
2052 Assert(
2053 (dynamic_cast<
2054 const typename ::MappingFE<dim, spacedim>::InternalData *>(
2055 &mapping_data) != nullptr),
2057 const typename ::MappingFE<dim, spacedim>::InternalData &data =
2058 static_cast<
2059 const typename ::MappingFE<dim, spacedim>::InternalData &>(
2060 mapping_data);
2061
2062 switch (mapping_kind)
2063 {
2064 case mapping_covariant:
2065 {
2066 Assert(
2067 data.update_each & update_contravariant_transformation,
2069 "update_covariant_transformation"));
2070
2071 for (unsigned int i = 0; i < output.size(); ++i)
2072 output[i] = apply_transformation(data.covariant[i], input[i]);
2073
2074 return;
2075 }
2076 default:
2078 }
2079 }
2080 } // namespace
2081 } // namespace MappingFEImplementation
2082} // namespace internal
2083
2084
2085
2086template <int dim, int spacedim>
2087void
2089 const ArrayView<const Tensor<1, dim>> &input,
2090 const MappingKind mapping_kind,
2091 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2092 const ArrayView<Tensor<1, spacedim>> &output) const
2093{
2094 internal::MappingFEImplementation::transform_fields(input,
2095 mapping_kind,
2096 mapping_data,
2097 output);
2098}
2099
2100
2101
2102template <int dim, int spacedim>
2103void
2105 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2106 const MappingKind mapping_kind,
2107 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2108 const ArrayView<Tensor<2, spacedim>> &output) const
2109{
2110 internal::MappingFEImplementation::transform_differential_forms(input,
2111 mapping_kind,
2112 mapping_data,
2113 output);
2114}
2115
2116
2117
2118template <int dim, int spacedim>
2119void
2121 const ArrayView<const Tensor<2, dim>> &input,
2122 const MappingKind mapping_kind,
2123 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2124 const ArrayView<Tensor<2, spacedim>> &output) const
2125{
2126 switch (mapping_kind)
2127 {
2129 internal::MappingFEImplementation::transform_fields(input,
2130 mapping_kind,
2131 mapping_data,
2132 output);
2133 return;
2134
2138 internal::MappingFEImplementation::transform_gradients(input,
2139 mapping_kind,
2140 mapping_data,
2141 output);
2142 return;
2143 default:
2145 }
2146}
2147
2148
2149
2150template <int dim, int spacedim>
2151void
2153 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2154 const MappingKind mapping_kind,
2155 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2156 const ArrayView<Tensor<3, spacedim>> &output) const
2157{
2158 AssertDimension(input.size(), output.size());
2159 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2161 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2162
2163 switch (mapping_kind)
2164 {
2166 {
2169 "update_covariant_transformation"));
2170
2171 for (unsigned int q = 0; q < output.size(); ++q)
2172 for (unsigned int i = 0; i < spacedim; ++i)
2173 for (unsigned int j = 0; j < spacedim; ++j)
2174 {
2175 double tmp[dim];
2176 for (unsigned int K = 0; K < dim; ++K)
2177 {
2178 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
2179 for (unsigned int J = 1; J < dim; ++J)
2180 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
2181 }
2182 for (unsigned int k = 0; k < spacedim; ++k)
2183 {
2184 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
2185 for (unsigned int K = 1; K < dim; ++K)
2186 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
2187 }
2188 }
2189 return;
2190 }
2191
2192 default:
2194 }
2195}
2196
2197
2198
2199template <int dim, int spacedim>
2200void
2202 const ArrayView<const Tensor<3, dim>> &input,
2203 const MappingKind mapping_kind,
2204 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2205 const ArrayView<Tensor<3, spacedim>> &output) const
2206{
2207 switch (mapping_kind)
2208 {
2212 internal::MappingFEImplementation::transform_hessians(input,
2213 mapping_kind,
2214 mapping_data,
2215 output);
2216 return;
2217 default:
2219 }
2220}
2221
2222
2223
2224namespace
2225{
2226 template <int spacedim>
2227 bool
2228 check_all_manifold_ids_identical(
2230 {
2231 return true;
2232 }
2233
2234
2235
2236 template <int spacedim>
2237 bool
2238 check_all_manifold_ids_identical(
2240 {
2241 const auto m_id = cell->manifold_id();
2242
2243 for (const auto f : cell->face_indices())
2244 if (m_id != cell->face(f)->manifold_id())
2245 return false;
2246
2247 return true;
2248 }
2249
2250
2251
2252 template <int spacedim>
2253 bool
2254 check_all_manifold_ids_identical(
2256 {
2257 const auto m_id = cell->manifold_id();
2258
2259 for (const auto f : cell->face_indices())
2260 if (m_id != cell->face(f)->manifold_id())
2261 return false;
2262
2263 for (const auto l : cell->line_indices())
2264 if (m_id != cell->line(l)->manifold_id())
2265 return false;
2266
2267 return true;
2268 }
2269} // namespace
2270
2271
2272
2273template <int dim, int spacedim>
2274std::vector<Point<spacedim>>
2276 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2277{
2278 Assert(
2279 check_all_manifold_ids_identical(cell),
2280 ExcMessage(
2281 "All entities of a cell need to have the same manifold id as the cell has."));
2282
2283 std::vector<Point<spacedim>> vertices(cell->n_vertices());
2284
2285 for (const unsigned int i : cell->vertex_indices())
2286 vertices[i] = cell->vertex(i);
2287
2288 std::vector<Point<spacedim>> mapping_support_points(
2289 fe->get_unit_support_points().size());
2290
2291 cell->get_manifold().get_new_points(vertices,
2292 mapping_support_point_weights,
2293 mapping_support_points);
2294
2295 return mapping_support_points;
2296}
2297
2298
2299
2300template <int dim, int spacedim>
2303 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2304{
2305 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
2306}
2307
2308
2309
2310template <int dim, int spacedim>
2311bool
2313 const ReferenceCell &reference_cell) const
2314{
2315 Assert(dim == reference_cell.get_dimension(),
2316 ExcMessage("The dimension of your mapping (" +
2318 ") and the reference cell cell_type (" +
2319 Utilities::to_string(reference_cell.get_dimension()) +
2320 " ) do not agree."));
2321
2322 return fe->reference_cell() == reference_cell;
2323}
2324
2325
2326
2327//--------------------------- Explicit instantiations -----------------------
2328#include "mapping_fe.inst"
2329
2330
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
DerivativeForm< 1, spacedim, dim, Number > transpose() 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:381
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_fe.cc:79
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition mapping_fe.h:387
void compute_shape_function_values(const std::vector< Point< dim > > &unit_points)
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition mapping_fe.h:370
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Definition mapping_fe.h:361
virtual std::size_t memory_consumption() const override
Definition mapping_fe.cc:59
InternalData(const FiniteElement< dim, spacedim > &fe)
Definition mapping_fe.cc:48
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:466
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:476
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:460
Abstract base class for mapping classes.
Definition mapping.h:316
Definition point.h:111
static DataSetDescriptor cell()
Definition qprojector.h:393
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int size() const
Definition collection.h:264
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
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)
Point< 2 > second
Definition grid_out.cc:4614
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:77
@ mapping_piola
Definition mapping.h:112
@ mapping_covariant_gradient
Definition mapping.h:98
@ mapping_covariant
Definition mapping.h:87
@ mapping_contravariant
Definition mapping.h:92
@ mapping_contravariant_hessian
Definition mapping.h:154
@ mapping_covariant_hessian
Definition mapping.h:148
@ mapping_contravariant_gradient
Definition mapping.h:104
@ mapping_piola_gradient
Definition mapping.h:118
@ mapping_piola_hessian
Definition mapping.h:160
#define DEAL_II_NOT_IMPLEMENTED()
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
if(marked_vertices.size() !=0) for(auto it
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:191
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:479
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void transform_gradients(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void transform_hessians(const ArrayView< const Tensor< 3, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim > > &output)
void maybe_compute_face_data(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int n_q_points, const std::vector< double > &weights, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int manifold_id
Definition types.h:156
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)