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