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