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