deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mapping_q_internal.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_mapping_q_internal_h
16#define dealii_mapping_q_internal_h
17
18#include <deal.II/base/config.h>
19
23#include <deal.II/base/point.h>
25#include <deal.II/base/table.h>
26#include <deal.II/base/tensor.h>
27
28#include <deal.II/fe/fe_tools.h>
32
34
41
42#include <array>
43#include <limits>
44#include <numeric>
45
46
48
49namespace internal
50{
56 namespace MappingQ1
57 {
58 // These are left as templates on the spatial dimension (even though dim
59 // == spacedim must be true for them to make sense) because templates are
60 // expanded before the compiler eliminates code due to the 'if (dim ==
61 // spacedim)' statement (see the body of the general
62 // transform_real_to_unit_cell).
63 template <int spacedim>
64 inline Point<1>
67 &vertices,
68 const Point<spacedim> &p)
69 {
70 Assert(spacedim == 1, ExcInternalError());
71 return Point<1>((p[0] - vertices[0][0]) /
72 (vertices[1][0] - vertices[0][0]));
73 }
74
75
76
77 template <int spacedim>
78 inline Point<2>
81 &vertices,
82 const Point<spacedim> &p)
83 {
84 Assert(spacedim == 2, ExcInternalError());
85
86 // For accuracy reasons, we do all arithmetic in extended precision
87 // (long double). This has a noticeable effect on the hit rate for
88 // borderline cases and thus makes the algorithm more robust.
89 const long double x = p[0];
90 const long double y = p[1];
91
92 const long double x0 = vertices[0][0];
93 const long double x1 = vertices[1][0];
94 const long double x2 = vertices[2][0];
95 const long double x3 = vertices[3][0];
96
97 const long double y0 = vertices[0][1];
98 const long double y1 = vertices[1][1];
99 const long double y2 = vertices[2][1];
100 const long double y3 = vertices[3][1];
101
102 const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3);
103 const long double b = -(x0 - x1 - x2 + x3) * y + (x - 2 * x1 + x3) * y0 -
104 (x - 2 * x0 + x2) * y1 - (x - x1) * y2 +
105 (x - x0) * y3;
106 const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1;
107
108 const long double discriminant = b * b - 4 * a * c;
109 // exit if the point is not in the cell (this is the only case where the
110 // discriminant is negative)
112 discriminant > 0.0,
114
115 long double eta1;
116 long double eta2;
117 const long double sqrt_discriminant = std::sqrt(discriminant);
118 // special case #1: if a is near-zero to make the discriminant exactly
119 // equal b, then use the linear formula
120 if (b != 0.0 && std::abs(b) == sqrt_discriminant)
121 {
122 eta1 = -c / b;
123 eta2 = -c / b;
124 }
125 // special case #2: a is zero for parallelograms and very small for
126 // near-parallelograms:
127 else if (std::abs(a) < 1e-8 * std::abs(b))
128 {
129 // if both a and c are very small then the root should be near
130 // zero: this first case will capture that
131 eta1 = 2 * c / (-b - sqrt_discriminant);
132 eta2 = 2 * c / (-b + sqrt_discriminant);
133 }
134 // finally, use the plain version:
135 else
136 {
137 eta1 = (-b - sqrt_discriminant) / (2 * a);
138 eta2 = (-b + sqrt_discriminant) / (2 * a);
139 }
140 // pick the one closer to the center of the cell.
141 const long double eta =
142 (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
143
144 /*
145 * There are two ways to compute xi from eta, but either one may have a
146 * zero denominator.
147 */
148 const long double subexpr0 = -eta * x2 + x0 * (eta - 1);
149 const long double xi_denominator0 = eta * x3 - x1 * (eta - 1) + subexpr0;
150 const long double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
152
153 if (std::abs(xi_denominator0) > 1e-10 * max_x)
154 {
155 const double xi = (x + subexpr0) / xi_denominator0;
156 return {xi, static_cast<double>(eta)};
157 }
158 else
159 {
160 const long double max_y =
163 const long double subexpr1 = -eta * y2 + y0 * (eta - 1);
164 const long double xi_denominator1 =
165 eta * y3 - y1 * (eta - 1) + subexpr1;
166 if (std::abs(xi_denominator1) > 1e-10 * max_y)
167 {
168 const double xi = (subexpr1 + y) / xi_denominator1;
169 return {xi, static_cast<double>(eta)};
170 }
171 else // give up and try Newton iteration
172 {
174 false,
175 (typename Mapping<spacedim,
176 spacedim>::ExcTransformationFailed()));
177 }
178 }
179 // bogus return to placate compiler. It should not be possible to get
180 // here.
182 return {std::numeric_limits<double>::quiet_NaN(),
183 std::numeric_limits<double>::quiet_NaN()};
184 }
185
186
187
188 template <int spacedim>
189 inline Point<3>
192 & /*vertices*/,
193 const Point<spacedim> & /*p*/)
194 {
195 // It should not be possible to get here
197 return {std::numeric_limits<double>::quiet_NaN(),
198 std::numeric_limits<double>::quiet_NaN(),
199 std::numeric_limits<double>::quiet_NaN()};
200 }
201 } // namespace MappingQ1
202
203
204
210 namespace MappingQImplementation
211 {
216 template <int dim>
217 std::vector<Point<dim>>
218 unit_support_points(const std::vector<Point<1>> &line_support_points,
219 const std::vector<unsigned int> &renumbering)
220 {
221 AssertDimension(Utilities::pow(line_support_points.size(), dim),
222 renumbering.size());
223 std::vector<Point<dim>> points(renumbering.size());
224 const unsigned int n1 = line_support_points.size();
225 for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2)
226 for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1)
227 for (unsigned int q0 = 0; q0 < n1; ++q0, ++q)
228 {
229 points[renumbering[q]][0] = line_support_points[q0][0];
230 if (dim > 1)
231 points[renumbering[q]][1] = line_support_points[q1][0];
232 if (dim > 2)
233 points[renumbering[q]][2] = line_support_points[q2][0];
234 }
235 return points;
236 }
237
238
239
247 inline ::Table<2, double>
248 compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
249 {
251
252 // we are asked to compute weights for interior support points, but
253 // there are no interior points if degree==1
254 if (polynomial_degree == 1)
255 return loqvs;
256
257 const unsigned int M = polynomial_degree - 1;
258 const unsigned int n_inner_2d = M * M;
259 const unsigned int n_outer_2d = 4 + 4 * M;
260
261 // set the weights of transfinite interpolation
262 loqvs.reinit(n_inner_2d, n_outer_2d);
263 QGaussLobatto<2> gl(polynomial_degree + 1);
264 for (unsigned int i = 0; i < M; ++i)
265 for (unsigned int j = 0; j < M; ++j)
266 {
267 const Point<2> &p =
268 gl.point((i + 1) * (polynomial_degree + 1) + (j + 1));
269 const unsigned int index_table = i * M + j;
270 for (unsigned int v = 0; v < 4; ++v)
271 loqvs(index_table, v) =
273 loqvs(index_table, 4 + i) = 1. - p[0];
274 loqvs(index_table, 4 + i + M) = p[0];
275 loqvs(index_table, 4 + j + 2 * M) = 1. - p[1];
276 loqvs(index_table, 4 + j + 3 * M) = p[1];
277 }
278
279 // the sum of weights of the points at the outer rim should be one.
280 // check this
281 for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
282 Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
283 loqvs[unit_point].end(),
284 0.) -
285 1) < 1e-13 * polynomial_degree,
287
288 return loqvs;
289 }
290
291
292
299 inline ::Table<2, double>
300 compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
301 {
303
304 // we are asked to compute weights for interior support points, but
305 // there are no interior points if degree==1
306 if (polynomial_degree == 1)
307 return lohvs;
308
309 const unsigned int M = polynomial_degree - 1;
310
311 const unsigned int n_inner = Utilities::fixed_power<3>(M);
312 const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
313
314 // set the weights of transfinite interpolation
315 lohvs.reinit(n_inner, n_outer);
316 QGaussLobatto<3> gl(polynomial_degree + 1);
317 for (unsigned int i = 0; i < M; ++i)
318 for (unsigned int j = 0; j < M; ++j)
319 for (unsigned int k = 0; k < M; ++k)
320 {
321 const Point<3> &p = gl.point((i + 1) * (M + 2) * (M + 2) +
322 (j + 1) * (M + 2) + (k + 1));
323 const unsigned int index_table = i * M * M + j * M + k;
324
325 // vertices
326 for (unsigned int v = 0; v < 8; ++v)
327 lohvs(index_table, v) =
329
330 // lines
331 {
332 constexpr std::array<unsigned int, 4> line_coordinates_y(
333 {{0, 1, 4, 5}});
334 const Point<2> py(p[0], p[2]);
335 for (unsigned int l = 0; l < 4; ++l)
336 lohvs(index_table, 8 + line_coordinates_y[l] * M + j) =
338 }
339
340 {
341 constexpr std::array<unsigned int, 4> line_coordinates_x(
342 {{2, 3, 6, 7}});
343 const Point<2> px(p[1], p[2]);
344 for (unsigned int l = 0; l < 4; ++l)
345 lohvs(index_table, 8 + line_coordinates_x[l] * M + k) =
347 }
348
349 {
350 constexpr std::array<unsigned int, 4> line_coordinates_z(
351 {{8, 9, 10, 11}});
352 const Point<2> pz(p[0], p[1]);
353 for (unsigned int l = 0; l < 4; ++l)
354 lohvs(index_table, 8 + line_coordinates_z[l] * M + i) =
356 }
357
358 // quads
359 lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
360 1. - p[0];
361 lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0];
362 lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) =
363 1. - p[1];
364 lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1];
365 lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) =
366 1. - p[2];
367 lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
368 }
369
370 // the sum of weights of the points at the outer rim should be one.
371 // check this
372 for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
373 Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
374 lohvs[unit_point].end(),
375 0.) -
376 1) < 1e-13 * polynomial_degree,
378
379 return lohvs;
380 }
381
382
383
388 inline std::vector<::Table<2, double>>
390 const unsigned int polynomial_degree,
391 const unsigned int dim)
392 {
393 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
394 std::vector<::Table<2, double>> output(dim);
395 if (polynomial_degree <= 1)
396 return output;
397
398 // fill the 1d interior weights
399 QGaussLobatto<1> quadrature(polynomial_degree + 1);
400 output[0].reinit(polynomial_degree - 1,
402 for (unsigned int q = 0; q < polynomial_degree - 1; ++q)
403 for (const unsigned int i : GeometryInfo<1>::vertex_indices())
404 output[0](q, i) =
405 GeometryInfo<1>::d_linear_shape_function(quadrature.point(q + 1),
406 i);
407
408 if (dim > 1)
409 output[1] = compute_support_point_weights_on_quad(polynomial_degree);
410
411 if (dim > 2)
412 output[2] = compute_support_point_weights_on_hex(polynomial_degree);
413
414 return output;
415 }
416
417
418
422 template <int dim>
423 inline ::Table<2, double>
424 compute_support_point_weights_cell(const unsigned int polynomial_degree)
425 {
426 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
427 if (polynomial_degree <= 1)
428 return ::Table<2, double>();
429
430 QGaussLobatto<dim> quadrature(polynomial_degree + 1);
431 const std::vector<unsigned int> h2l =
432 FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
433
434 ::Table<2, double> output(quadrature.size() -
437 for (unsigned int q = 0; q < output.size(0); ++q)
438 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
440 quadrature.point(h2l[q + GeometryInfo<dim>::vertices_per_cell]), i);
441
442 return output;
443 }
444
445
446
454 template <int dim, int spacedim>
455 inline Point<spacedim>
457 const typename ::MappingQ<dim, spacedim>::InternalData &data)
458 {
459 AssertDimension(data.shape_values.size(),
460 data.mapping_support_points.size());
461
462 // use now the InternalData to compute the point in real space.
464 for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i)
465 p_real += data.mapping_support_points[i] * data.shape(0, i);
466
467 return p_real;
468 }
469
470
471
476 template <int dim, int spacedim, typename Number>
477 inline Point<dim, Number>
481 const ArrayView<const Point<spacedim>> &points,
482 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
483 const std::vector<unsigned int> &renumber,
484 const bool print_iterations_to_deallog = false)
485 {
487 deallog << "Start MappingQ::do_transform_real_to_unit_cell for real "
488 << "point [ " << p << " ] " << std::endl;
489
490 AssertDimension(points.size(),
491 Utilities::pow(polynomials_1d.size(), dim));
492
493 // Newton iteration to solve
494 // f(x)=p(x)-p=0
495 // where we are looking for 'x' and p(x) is the forward transformation
496 // from unit to real cell. We solve this using a Newton iteration
497 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
498 // The start value is set to be the linear approximation to the cell
499
500 // The shape values and derivatives of the mapping at this point are
501 // previously computed.
502
505 polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber);
506
508
509 // early out if we already have our point in all SIMD lanes, i.e.,
510 // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable
511 // this step also for VectorizedArray where we do not have operator <,
512 // compare the result to zero which is defined for SIMD types
513 if (std::max(Number(0.),
514 f.norm_square() - 1e-24 * p_real.second[0].norm_square()) ==
515 Number(0.))
516 return p_unit;
517
518 // we need to compare the position of the computed p(x) against the
519 // given point 'p'. We will terminate the iteration and return 'x' if
520 // they are less than eps apart. The question is how to choose eps --
521 // or, put maybe more generally: in which norm we want these 'p' and
522 // 'p(x)' to be eps apart.
523 //
524 // the question is difficult since we may have to deal with very
525 // elongated cells where we may achieve 1e-12*h for the distance of
526 // these two points in the 'long' direction, but achieving this
527 // tolerance in the 'short' direction of the cell may not be possible
528 //
529 // what we do instead is then to terminate iterations if
530 // \| p(x) - p \|_A < eps
531 // where the A-norm is somehow induced by the transformation of the
532 // cell. in particular, we want to measure distances relative to the
533 // sizes of the cell in its principal directions.
534 //
535 // to define what exactly A should be, note that to first order we have
536 // the following (assuming that x* is the solution of the problem, i.e.,
537 // p(x*)=p):
538 // p(x) - p = p(x) - p(x*)
539 // = -grad p(x) * (x*-x) + higher order terms
540 // This suggest to measure with a norm that corresponds to
541 // A = {[grad p(x)]^T [grad p(x)]}^{-1}
542 // because then
543 // \| p(x) - p \|_A \approx \| x - x* \|
544 // Consequently, we will try to enforce that
545 // \| p(x) - p \|_A = \| f \| <= eps
546 //
547 // Note that using this norm is a bit dangerous since the norm changes
548 // in every iteration (A isn't fixed by depending on xk). However, if
549 // the cell is not too deformed (it may be stretched, but not twisted)
550 // then the mapping is almost linear and A is indeed constant or
551 // nearly so.
552 const double eps = 1.e-11;
553 const unsigned int newton_iteration_limit = 20;
554
556 invalid_point[0] = std::numeric_limits<double>::lowest();
557 bool tried_project_to_unit_cell = false;
558
559 unsigned int newton_iteration = 0;
560 Number f_weighted_norm_square = 1.;
561 Number last_f_weighted_norm_square = 1.;
562
563 do
564 {
566 deallog << "Newton iteration " << newton_iteration
567 << " for unit point guess " << p_unit << std::endl;
568
569 // f'(x)
571 for (unsigned int d = 0; d < spacedim; ++d)
572 for (unsigned int e = 0; e < dim; ++e)
573 df[d][e] = p_real.second[e][d];
574
575 // Check determinant(df) > 0 on all SIMD lanes. The
576 // condition here is unreadable (but really corresponds to
577 // asking whether det(df) > 0 for all elements of the
578 // vector) because VectorizedArray does not have a member
579 // that can be used to check that all vector elements are
580 // positive. But VectorizedArray has a std::min() function,
581 // and operator==().
582 if (!(std::min(determinant(df),
583 Number(std::numeric_limits<double>::min())) ==
584 Number(std::numeric_limits<double>::min())))
585 {
586 // We allow to enter this function with an initial guess
587 // outside the unit cell. We might have run into invalid
588 // Jacobians due to that, so we should at least try once to go
589 // back to the unit cell and go on with the Newton iteration
590 // from there. Since the outside case is unlikely, we can
591 // afford spending the extra effort at this place.
592 if (tried_project_to_unit_cell == false)
593 {
596 polynomials_1d,
597 points,
598 p_unit,
599 polynomials_1d.size() == 2,
600 renumber);
601 f = p_real.first - p;
605 continue;
606 }
607 else
608 return invalid_point;
609 }
610
611 // Solve [f'(x)]d=f(x)
613 const Tensor<1, spacedim, Number> delta = df_inverse * f;
615
617 deallog << " delta=" << delta << std::endl;
618
619 // do a line search
620 double step_length = 1.0;
621 do
622 {
623 // update of p_unit. The spacedim-th component of transformed
624 // point is simply ignored in codimension one case. When this
625 // component is not zero, then we are projecting the point to
626 // the surface or curve identified by the cell.
628 for (unsigned int i = 0; i < dim; ++i)
629 p_unit_trial[i] -= step_length * delta[i];
630
631 // shape values and derivatives at new p_unit point
632 const auto p_real_trial =
634 polynomials_1d,
635 points,
637 polynomials_1d.size() == 2,
638 renumber);
640 p_real_trial.first - p;
641 f_weighted_norm_square = (df_inverse * f_trial).norm_square();
642
644 {
645 deallog << " step_length=" << step_length << std::endl;
646 if (step_length == 1.0)
647 deallog << " ||f || =" << f.norm() << std::endl;
648 deallog << " ||f*|| =" << f_trial.norm() << std::endl
649 << " ||f*||_A ="
650 << std::sqrt(f_weighted_norm_square) << std::endl;
651 }
652
653 // See if we are making progress with the current step length
654 // and if not, reduce it by a factor of two and try again.
655 //
656 // Strictly speaking, we should probably use the same norm as we
657 // use for the outer algorithm. In practice, line search is just
658 // a crutch to find a "reasonable" step length, and so using the
659 // l2 norm is probably just fine.
660 //
661 // check f_trial.norm() < f.norm() in SIMD form. This is a bit
662 // more complicated because some SIMD lanes might not be doing
663 // any progress any more as they have already reached roundoff
664 // accuracy. We define that as the case of updates less than
665 // 1e-6. The tolerance might seem coarse but since we are
666 // dealing with a Newton iteration of a polynomial function we
667 // either converge quadratically or not any more. Thus, our
668 // selection is to terminate if either the norm of f is
669 // decreasing or that threshold of 1e-6 is reached.
670 if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) *
672 Number(0.)) ==
673 Number(0.))
674 {
677 f = f_trial;
678 break;
679 }
680 else if (step_length > 0.05)
681 step_length *= 0.5;
682 else
683 break;
684 }
685 while (true);
686
687 // In case we terminated the line search due to the step becoming
688 // too small, we give the iteration another try with the
689 // projection of the initial guess to the unit cell before we give
690 // up, just like for the negative determinant case.
691 if (step_length <= 0.05 && tried_project_to_unit_cell == false)
692 {
695 polynomials_1d,
696 points,
697 p_unit,
698 polynomials_1d.size() == 2,
699 renumber);
700 f = p_real.first - p;
704 continue;
705 }
706 else if (step_length <= 0.05)
707 return invalid_point;
708
711 return invalid_point;
712 }
713 // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the
714 // weighted norm is less than 1e-6 and the improvement against the
715 // previous step was less than a factor of 10 (in that regime, we
716 // either have quadratic convergence or roundoff errors due to a bad
717 // mapping)
718 while (
719 !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) *
721 std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) *
722 100.,
723 Number(0.)) ==
724 Number(0.)));
725
727 deallog << "Iteration converged for p_unit = [ " << p_unit
728 << " ] and iteration error "
729 << std::sqrt(f_weighted_norm_square) << std::endl;
730
731 return p_unit;
732 }
733
734
735
739 template <int dim>
740 inline Point<dim>
742 const Point<dim + 1> &p,
744 const ArrayView<const Point<dim + 1>> &points,
745 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
746 const std::vector<unsigned int> &renumber)
747 {
748 const unsigned int spacedim = dim + 1;
749
750 AssertDimension(points.size(),
751 Utilities::pow(polynomials_1d.size(), dim));
752
754
755 const double eps = 1.e-12;
756 const unsigned int loop_limit = 10;
757
758 unsigned int loop = 0;
759 double f_weighted_norm_square = 1.;
760
761 while (f_weighted_norm_square > eps * eps && loop++ < loop_limit)
762 {
763 const auto p_real =
765 polynomials_1d,
766 points,
767 p_unit,
768 polynomials_1d.size() == 2,
769 renumber);
772
774 polynomials_1d, points, p_unit, renumber);
775 Point<dim> f;
777 for (unsigned int j = 0; j < dim; ++j)
778 {
779 f[j] = DF[j] * p_minus_F;
780 for (unsigned int l = 0; l < dim; ++l)
781 df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F;
782 }
783
784 // Solve [df(x)]d=f(x)
785 const Tensor<1, dim> d =
786 invert(df) * static_cast<const Tensor<1, dim> &>(f);
788 p_unit -= d;
789 }
790
791
792 // Here we check that in the last execution of while the first
793 // condition was already wrong, meaning the residual was below
794 // eps. Only if the first condition failed, loop will have been
795 // increased and tested, and thus have reached the limit.
796 AssertThrow(loop < loop_limit,
798
799 return p_unit;
800 }
801
802
803
821 template <int dim, int spacedim>
823 {
824 public:
828 static constexpr unsigned int n_functions =
829 (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
830
844 const std::vector<Point<dim>> &unit_support_points)
847 1. / real_support_points[0].distance(real_support_points[1]))
848 , is_affine(true)
849 {
851
852 // For the bi-/trilinear approximation, we cannot build a quadratic
853 // polynomial due to a lack of points (interpolation matrix would
854 // get singular). Similarly, it is not entirely clear how to gather
855 // enough information for the case dim < spacedim.
856 //
857 // In both cases we require the vector real_support_points to
858 // contain the vertex positions and fall back to an affine
859 // approximation:
860 Assert(dim == spacedim || real_support_points.size() ==
864 {
865 const auto affine = GridTools::affine_cell_approximation<dim>(
868 affine.first.covariant_form().transpose();
869
870 // The code for evaluation assumes an additional transformation of
871 // the form (x - normalization_shift) * normalization_length --
872 // account for this in the definition of the coefficients.
873 coefficients[0] =
875 for (unsigned int d = 0; d < spacedim; ++d)
876 for (unsigned int e = 0; e < dim; ++e)
877 coefficients[1 + d][e] =
878 A_inv[e][d] * (1.0 / normalization_length);
879 is_affine = true;
880 return;
881 }
882
884 std::array<double, n_functions> shape_values;
885 for (unsigned int q = 0; q < unit_support_points.size(); ++q)
886 {
887 // Evaluate quadratic shape functions in point, with the
888 // normalization applied in order to avoid roundoff issues with
889 // scaling far away from 1.
890 shape_values[0] = 1.;
894 for (unsigned int d = 0; d < spacedim; ++d)
895 shape_values[1 + d] = p_scaled[d];
896 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
897 for (unsigned int e = 0; e <= d; ++e, ++c)
898 shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
899
900 // Build lower diagonal of least squares matrix and rhs, the
901 // essential part being that we construct the matrix with the
902 // real points and the right hand side by comparing to the
903 // reference point positions which sets up an inverse
904 // interpolation.
905 for (unsigned int i = 0; i < n_functions; ++i)
906 for (unsigned int j = 0; j < n_functions; ++j)
907 matrix[i][j] += shape_values[i] * shape_values[j];
908 for (unsigned int i = 0; i < n_functions; ++i)
909 coefficients[i] += shape_values[i] * unit_support_points[q];
910 }
911
912 // Factorize the matrix A = L * L^T in-place with the
913 // Cholesky-Banachiewicz algorithm. The implementation is similar to
914 // FullMatrix::cholesky() but re-implemented to avoid memory
915 // allocations and some unnecessary divisions which we can do here as
916 // we only need to solve with dim right hand sides.
917 for (unsigned int i = 0; i < n_functions; ++i)
918 {
919 double Lij_sum = 0;
920 for (unsigned int j = 0; j < i; ++j)
921 {
922 double Lik_Ljk_sum = 0;
923 for (unsigned int k = 0; k < j; ++k)
924 Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
925 matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
926 Lij_sum += matrix[i][j] * matrix[i][j];
927 }
928 AssertThrow(matrix[i][i] - Lij_sum >= 0,
929 ExcMessage("Matrix of normal equations not positive "
930 "definite"));
931
932 // Store the inverse in the diagonal since that is the quantity
933 // needed later in the factorization as well as the forward and
934 // backward substitution, minimizing the number of divisions.
935 matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
936 }
937
938 // Solve lower triangular part, L * y = rhs.
939 for (unsigned int i = 0; i < n_functions; ++i)
940 {
941 Point<dim> sum = coefficients[i];
942 for (unsigned int j = 0; j < i; ++j)
943 sum -= matrix[i][j] * coefficients[j];
944 coefficients[i] = sum * matrix[i][i];
945 }
946
947 // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
948 for (unsigned int i = n_functions; i > 0;)
949 {
950 --i;
951 Point<dim> sum = coefficients[i];
952 for (unsigned int j = i + 1; j < n_functions; ++j)
953 sum -= matrix[j][i] * coefficients[j];
954 coefficients[i] = sum * matrix[i][i];
955 }
956
957 // Check whether the approximation is indeed affine, allowing to
958 // skip the quadratic terms.
959 is_affine = true;
960 for (unsigned int i = dim + 1; i < n_functions; ++i)
961 if (coefficients[i].norm_square() > 1e-20)
962 {
963 is_affine = false;
964 break;
965 }
966 }
967
972 default;
973
977 template <typename Number>
980 {
982 for (unsigned int d = 0; d < dim; ++d)
983 result[d] = coefficients[0][d];
984
985 // Apply the normalization to ensure a good conditioning. Since Number
986 // might be a vectorized array whereas the normalization is a point of
987 // doubles, we cannot use the overload of operator- and must instead
988 // loop over the components of the point.
990 for (unsigned int d = 0; d < spacedim; ++d)
992
993 for (unsigned int d = 0; d < spacedim; ++d)
994 result += coefficients[1 + d] * p_scaled[d];
995
996 if (!is_affine)
997 {
999 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
1000 for (unsigned int e = 0; e <= d; ++e, ++c)
1001 result +=
1002 coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
1003
1004 // Check if the quadratic approximation ends up considerably
1005 // farther outside the unit cell on some or all SIMD lanes than
1006 // the affine approximation - in that case, we switch those
1007 // components back to the affine approximation. Note that the
1008 // quadratic approximation will grow more quickly away from the
1009 // unit cell. We make the selection for each SIMD lane with a
1010 // ternary operation.
1011 const Number distance_to_unit_cell = result.distance_square(
1013 const Number affine_distance_to_unit_cell =
1014 result_affine.distance_square(
1016 for (unsigned int d = 0; d < dim; ++d)
1018 distance_to_unit_cell,
1020 result_affine[d],
1021 result[d]);
1022 }
1023 return result;
1024 }
1025
1026 private:
1035
1040
1044 std::array<Point<dim>, n_functions> coefficients;
1045
1052 };
1053
1054
1055
1061 template <int dim, int spacedim>
1062 inline void
1064 const CellSimilarity::Similarity cell_similarity,
1065 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1066 std::vector<Point<spacedim>> &quadrature_points,
1067 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1068 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians,
1069 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1070 {
1071 const UpdateFlags update_flags = data.update_each;
1072
1073 using VectorizedArrayType =
1074 typename ::MappingQ<dim,
1075 spacedim>::InternalData::VectorizedArrayType;
1076 const unsigned int n_shape_values = data.n_shape_functions;
1077 const unsigned int n_q_points = data.shape_info.n_q_points;
1078 constexpr unsigned int n_lanes = VectorizedArrayType::size();
1079 constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1080 constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1081
1082 // Since MappingQ::InternalData does not have separate arrays for the
1083 // covariant and contravariant transformations, but uses the arrays in
1084 // the `MappingRelatedData`, it can happen that vectors do not have the
1085 // right size
1086 if (update_flags & update_contravariant_transformation)
1087 jacobians.resize(n_q_points);
1088 if (update_flags & update_covariant_transformation)
1089 inverse_jacobians.resize(n_q_points);
1090
1094 ((cell_similarity != CellSimilarity::translation) &&
1095 (update_flags & update_contravariant_transformation) ?
1098 ((cell_similarity != CellSimilarity::translation) &&
1099 (update_flags & update_jacobian_grads) ?
1102
1103 Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1106 n_q_points == quadrature_points.size(),
1107 ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1109 data.n_shape_functions > 0,
1112 n_q_points == jacobian_grads.size(),
1113 ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1114
1115 // shortcut in case we have an identity interpolation and only request
1116 // the quadrature points
1118 data.shape_info.element_type ==
1120 {
1121 for (unsigned int q = 0; q < n_q_points; ++q)
1122 quadrature_points[q] =
1123 data.mapping_support_points[data.shape_info
1124 .lexicographic_numbering[q]];
1125 return;
1126 }
1127
1129
1130 // prepare arrays
1132 {
1133 eval.set_data_pointers(&data.scratch, n_comp);
1134
1135 // make sure to initialize on all lanes also when some are unused in
1136 // the code below
1137 for (unsigned int i = 0; i < n_shape_values * n_comp; ++i)
1138 eval.begin_dof_values()[i] = VectorizedArrayType();
1139
1140 const std::vector<unsigned int> &renumber_to_lexicographic =
1141 data.shape_info.lexicographic_numbering;
1142 for (unsigned int i = 0; i < n_shape_values; ++i)
1143 for (unsigned int d = 0; d < spacedim; ++d)
1144 {
1145 const unsigned int in_comp = d % n_lanes;
1146 const unsigned int out_comp = d / n_lanes;
1147 eval
1148 .begin_dof_values()[out_comp * n_shape_values + i][in_comp] =
1149 data.mapping_support_points[renumber_to_lexicographic[i]][d];
1150 }
1151
1152 // do the actual tensorized evaluation
1154 n_comp, evaluation_flag, eval.begin_dof_values(), eval);
1155 }
1156
1157 // do the postprocessing
1159 {
1160 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1161 for (unsigned int i = 0; i < n_q_points; ++i)
1162 for (unsigned int in_comp = 0;
1163 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1164 ++in_comp)
1165 quadrature_points[i][out_comp * n_lanes + in_comp] =
1166 eval.begin_values()[out_comp * n_q_points + i][in_comp];
1167 }
1168
1170 {
1171 // We need to reinterpret the data after evaluate has been applied.
1172 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1173 for (unsigned int point = 0; point < n_q_points; ++point)
1174 for (unsigned int in_comp = 0;
1175 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1176 ++in_comp)
1177 for (unsigned int j = 0; j < dim; ++j)
1178 {
1179 jacobians[point][out_comp * n_lanes + in_comp][j] =
1180 eval.begin_gradients()[(out_comp * n_q_points + point) *
1181 dim +
1182 j][in_comp];
1183 }
1184 }
1185
1186 if (update_flags & update_volume_elements)
1187 if (cell_similarity != CellSimilarity::translation)
1188 for (unsigned int point = 0; point < n_q_points; ++point)
1189 data.volume_elements[point] = jacobians[point].determinant();
1190
1191 // copy values from InternalData to vector given by reference
1192 if (update_flags & update_covariant_transformation)
1193 {
1194 AssertDimension(jacobians.size(), n_q_points);
1195 AssertDimension(inverse_jacobians.size(), n_q_points);
1196 if (cell_similarity != CellSimilarity::translation)
1197 for (unsigned int point = 0; point < n_q_points; ++point)
1198 inverse_jacobians[point] =
1199 jacobians[point].covariant_form().transpose();
1200 }
1201
1203 {
1204 constexpr int desymmetrize_3d[6][2] = {
1205 {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1206 constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1207
1208 // We need to reinterpret the data after evaluate has been applied.
1209 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1210 for (unsigned int point = 0; point < n_q_points; ++point)
1211 for (unsigned int j = 0; j < n_hessians; ++j)
1212 for (unsigned int in_comp = 0;
1213 in_comp < n_lanes &&
1214 in_comp < spacedim - out_comp * n_lanes;
1215 ++in_comp)
1216 {
1217 const unsigned int total_number = point * n_hessians + j;
1218 const unsigned int new_point = total_number % n_q_points;
1219 const unsigned int new_hessian_comp =
1220 total_number / n_q_points;
1221 const unsigned int new_hessian_comp_i =
1222 dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1224 const unsigned int new_hessian_comp_j =
1225 dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1227 const double value =
1228 eval.begin_hessians()[(out_comp * n_q_points + point) *
1229 n_hessians +
1230 j][in_comp];
1231 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1233 value;
1234 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1236 value;
1237 }
1238 }
1239 }
1240
1241
1242
1243 template <int dim, int spacedim>
1244 inline void
1246 const CellSimilarity::Similarity cell_similarity,
1247 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1248 const ArrayView<const Point<dim>> &unit_points,
1249 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1250 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1251 std::vector<Point<spacedim>> &quadrature_points,
1252 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1253 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
1254 {
1255 const UpdateFlags update_flags = data.update_each;
1257 data.mapping_support_points);
1258
1259 const unsigned int n_points = unit_points.size();
1260 const unsigned int n_lanes = VectorizedArray<double>::size();
1261
1262 // Since MappingQ::InternalData does not have separate arrays for the
1263 // covariant and contravariant transformations, but uses the arrays in
1264 // the `MappingRelatedData`, it can happen that vectors do not have the
1265 // right size
1266 if (update_flags & update_contravariant_transformation)
1267 jacobians.resize(n_points);
1268 if (update_flags & update_covariant_transformation)
1269 inverse_jacobians.resize(n_points);
1270
1271 const bool is_translation =
1272 cell_similarity == CellSimilarity::translation;
1273
1274 const bool needs_gradient =
1275 !is_translation &&
1276 (update_flags &
1279
1280 // Use the more heavy VectorizedArray code path if there is more than
1281 // one point left to compute
1282 for (unsigned int i = 0; i < n_points; i += n_lanes)
1283 if (n_points - i > 1)
1284 {
1286 for (unsigned int j = 0; j < n_lanes; ++j)
1287 if (i + j < n_points)
1288 for (unsigned int d = 0; d < dim; ++d)
1289 p_vec[d][j] = unit_points[i + j][d];
1290 else
1291 for (unsigned int d = 0; d < dim; ++d)
1292 p_vec[d][j] = unit_points[i][d];
1293
1296 derivative;
1297 if (needs_gradient)
1298 {
1299 const auto result =
1301 polynomials_1d,
1303 p_vec,
1304 polynomials_1d.size() == 2,
1305 renumber_lexicographic_to_hierarchic);
1306
1307 value = result.first;
1308
1309 for (unsigned int d = 0; d < spacedim; ++d)
1310 for (unsigned int e = 0; e < dim; ++e)
1311 derivative[d][e] = result.second[e][d];
1312 }
1313 else
1315 polynomials_1d,
1317 p_vec,
1318 polynomials_1d.size() == 2,
1319 renumber_lexicographic_to_hierarchic);
1320
1321 if (update_flags & update_quadrature_points)
1322 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1323 for (unsigned int d = 0; d < spacedim; ++d)
1324 quadrature_points[i + j][d] = value[d][j];
1325
1326 if (is_translation)
1327 continue;
1328
1329 if (update_flags & update_contravariant_transformation)
1330 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1331 for (unsigned int d = 0; d < spacedim; ++d)
1332 for (unsigned int e = 0; e < dim; ++e)
1333 jacobians[i + j][d][e] = derivative[d][e][j];
1334
1335 if (update_flags & update_volume_elements)
1336 {
1337 const auto determinant = derivative.determinant();
1338 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1339 data.volume_elements[i + j] = determinant[j];
1340 }
1341
1342 if (update_flags & update_covariant_transformation)
1343 {
1344 const auto covariant = derivative.covariant_form();
1345 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1346 for (unsigned int d = 0; d < dim; ++d)
1347 for (unsigned int e = 0; e < spacedim; ++e)
1348 inverse_jacobians[i + j][d][e] = covariant[e][d][j];
1349 }
1350 }
1351 else
1352 {
1355 if (needs_gradient)
1356 {
1357 const auto result =
1359 polynomials_1d,
1361 unit_points[i],
1362 polynomials_1d.size() == 2,
1363 renumber_lexicographic_to_hierarchic);
1364
1365 value = result.first;
1366
1367 for (unsigned int d = 0; d < spacedim; ++d)
1368 for (unsigned int e = 0; e < dim; ++e)
1369 derivative[d][e] = result.second[e][d];
1370 }
1371 else
1373 polynomials_1d,
1375 unit_points[i],
1376 polynomials_1d.size() == 2,
1377 renumber_lexicographic_to_hierarchic);
1378
1379 if (update_flags & update_quadrature_points)
1380 quadrature_points[i] = value;
1381
1382 if (is_translation)
1383 continue;
1384
1385 if (dim == spacedim && update_flags & update_volume_elements)
1386 data.volume_elements[i] = derivative.determinant();
1387
1388 if (update_flags & update_contravariant_transformation)
1389 jacobians[i] = derivative;
1390
1391 if (update_flags & update_covariant_transformation)
1392 inverse_jacobians[i] = jacobians[i].covariant_form().transpose();
1393 }
1394 }
1395
1396
1397
1404 template <int dim, int spacedim>
1405 inline void
1407 const CellSimilarity::Similarity cell_similarity,
1408 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1409 const ArrayView<const Point<dim>> &unit_points,
1410 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1411 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1412 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1413 {
1414 if (data.update_each & update_jacobian_grads)
1415 {
1417 data.mapping_support_points);
1418 const unsigned int n_q_points = jacobian_grads.size();
1419
1420 if (cell_similarity != CellSimilarity::translation)
1421 for (unsigned int point = 0; point < n_q_points; ++point)
1422 {
1425 polynomials_1d,
1427 unit_points[point],
1428 renumber_lexicographic_to_hierarchic);
1429
1430 for (unsigned int i = 0; i < spacedim; ++i)
1431 for (unsigned int j = 0; j < dim; ++j)
1432 for (unsigned int l = 0; l < dim; ++l)
1433 jacobian_grads[point][i][j][l] = second[j][l][i];
1434 }
1435 }
1436 }
1437
1438
1439
1446 template <int dim, int spacedim>
1447 inline void
1449 const CellSimilarity::Similarity cell_similarity,
1450 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1451 const ArrayView<const Point<dim>> &unit_points,
1452 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1453 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1454 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1455 {
1457 {
1459 data.mapping_support_points);
1460 const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
1461
1462 if (cell_similarity != CellSimilarity::translation)
1463 {
1464 double tmp[spacedim][spacedim][spacedim];
1465 for (unsigned int point = 0; point < n_q_points; ++point)
1466 {
1469 polynomials_1d,
1471 unit_points[point],
1472 renumber_lexicographic_to_hierarchic);
1473 const DerivativeForm<1, dim, spacedim> covariant =
1474 data.output_data->inverse_jacobians[point].transpose();
1475
1476 // first push forward the j-components
1477 for (unsigned int i = 0; i < spacedim; ++i)
1478 for (unsigned int j = 0; j < spacedim; ++j)
1479 for (unsigned int l = 0; l < dim; ++l)
1480 {
1481 tmp[i][j][l] = second[0][l][i] * covariant[j][0];
1482 for (unsigned int jr = 1; jr < dim; ++jr)
1483 {
1484 tmp[i][j][l] +=
1485 second[jr][l][i] * covariant[j][jr];
1486 }
1487 }
1488
1489 // now, pushing forward the l-components
1490 for (unsigned int i = 0; i < spacedim; ++i)
1491 for (unsigned int j = 0; j < spacedim; ++j)
1492 for (unsigned int l = 0; l < spacedim; ++l)
1493 {
1494 jacobian_pushed_forward_grads[point][i][j][l] =
1495 tmp[i][j][0] * covariant[l][0];
1496 for (unsigned int lr = 1; lr < dim; ++lr)
1497 {
1498 jacobian_pushed_forward_grads[point][i][j][l] +=
1499 tmp[i][j][lr] * covariant[l][lr];
1500 }
1501 }
1502 }
1503 }
1504 }
1505 }
1506
1507
1508
1509 template <int dim, int spacedim, int length_tensor>
1512 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1513 {
1514 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1516 for (unsigned int i = 0; i < spacedim; ++i)
1517 {
1518 if (dim == 1)
1519 result[i][0][0][0] = compressed[0][i];
1520 else if (dim >= 2)
1521 {
1522 for (unsigned int d = 0; d < 2; ++d)
1523 for (unsigned int e = 0; e < 2; ++e)
1524 for (unsigned int f = 0; f < 2; ++f)
1525 result[i][d][e][f] = compressed[d + e + f][i];
1526 if (dim == 3)
1527 {
1529
1530 // the derivatives are ordered in rows by increasing z
1531 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1532 // j walks through the columns
1533 for (unsigned int d = 0; d < 2; ++d)
1534 for (unsigned int e = 0; e < 2; ++e)
1535 {
1536 result[i][d][e][2] = compressed[4 + d + e][i];
1537 result[i][d][2][e] = compressed[4 + d + e][i];
1538 result[i][2][d][e] = compressed[4 + d + e][i];
1539 }
1540 for (unsigned int d = 0; d < 2; ++d)
1541 {
1542 result[i][d][2][2] = compressed[7 + d][i];
1543 result[i][2][d][2] = compressed[7 + d][i];
1544 result[i][2][2][d] = compressed[7 + d][i];
1545 }
1546 result[i][2][2][2] = compressed[9][i];
1547 }
1548 }
1549 }
1550 return result;
1551 }
1552
1553
1554
1561 template <int dim, int spacedim>
1562 inline void
1564 const CellSimilarity::Similarity cell_similarity,
1565 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1566 const ArrayView<const Point<dim>> &unit_points,
1567 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1568 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1569 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1570 {
1571 if (data.update_each & update_jacobian_2nd_derivatives)
1572 {
1574 data.mapping_support_points);
1575 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1576
1577 if (cell_similarity != CellSimilarity::translation)
1578 {
1579 for (unsigned int point = 0; point < n_q_points; ++point)
1580 {
1581 jacobian_2nd_derivatives[point] = expand_3rd_derivatives<dim>(
1582 internal::evaluate_tensor_product_higher_derivatives<3>(
1583 polynomials_1d,
1585 unit_points[point],
1586 renumber_lexicographic_to_hierarchic));
1587 }
1588 }
1589 }
1590 }
1591
1592
1593
1601 template <int dim, int spacedim>
1602 inline void
1604 const CellSimilarity::Similarity cell_similarity,
1605 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1606 const ArrayView<const Point<dim>> &unit_points,
1607 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1608 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1609 std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
1610 {
1612 {
1614 data.mapping_support_points);
1615 const unsigned int n_q_points =
1616 jacobian_pushed_forward_2nd_derivatives.size();
1617
1618 if (cell_similarity != CellSimilarity::translation)
1619 {
1622 for (unsigned int point = 0; point < n_q_points; ++point)
1623 {
1626 internal::evaluate_tensor_product_higher_derivatives<3>(
1627 polynomials_1d,
1629 unit_points[point],
1630 renumber_lexicographic_to_hierarchic));
1631 const DerivativeForm<1, dim, spacedim> covariant =
1632 data.output_data->inverse_jacobians[point].transpose();
1633
1634 // push forward the j-coordinate
1635 for (unsigned int i = 0; i < spacedim; ++i)
1636 for (unsigned int j = 0; j < spacedim; ++j)
1637 for (unsigned int l = 0; l < dim; ++l)
1638 for (unsigned int m = 0; m < dim; ++m)
1639 {
1640 tmp[i][j][l][m] =
1641 third[i][0][l][m] * covariant[j][0];
1642 for (unsigned int jr = 1; jr < dim; ++jr)
1643 tmp[i][j][l][m] +=
1644 third[i][jr][l][m] * covariant[j][jr];
1645 }
1646
1647 // push forward the l-coordinate
1648 for (unsigned int i = 0; i < spacedim; ++i)
1649 for (unsigned int j = 0; j < spacedim; ++j)
1650 for (unsigned int l = 0; l < spacedim; ++l)
1651 for (unsigned int m = 0; m < dim; ++m)
1652 {
1653 tmp2[i][j][l][m] =
1654 tmp[i][j][0][m] * covariant[l][0];
1655 for (unsigned int lr = 1; lr < dim; ++lr)
1656 tmp2[i][j][l][m] +=
1657 tmp[i][j][lr][m] * covariant[l][lr];
1658 }
1659
1660 // push forward the m-coordinate
1661 for (unsigned int i = 0; i < spacedim; ++i)
1662 for (unsigned int j = 0; j < spacedim; ++j)
1663 for (unsigned int l = 0; l < spacedim; ++l)
1664 for (unsigned int m = 0; m < spacedim; ++m)
1665 {
1666 jacobian_pushed_forward_2nd_derivatives
1667 [point][i][j][l][m] =
1668 tmp2[i][j][l][0] * covariant[m][0];
1669 for (unsigned int mr = 1; mr < dim; ++mr)
1670 jacobian_pushed_forward_2nd_derivatives[point][i]
1671 [j][l]
1672 [m] +=
1673 tmp2[i][j][l][mr] * covariant[m][mr];
1674 }
1675 }
1676 }
1677 }
1678 }
1679
1680
1681
1682 template <int dim, int spacedim, int length_tensor>
1685 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1686 {
1687 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1689 for (unsigned int i = 0; i < spacedim; ++i)
1690 {
1691 if (dim == 1)
1692 result[i][0][0][0][0] = compressed[0][i];
1693 else if (dim >= 2)
1694 {
1695 for (unsigned int d = 0; d < 2; ++d)
1696 for (unsigned int e = 0; e < 2; ++e)
1697 for (unsigned int f = 0; f < 2; ++f)
1698 for (unsigned int g = 0; g < 2; ++g)
1699 result[i][d][e][f][g] = compressed[d + e + f + g][i];
1700 if (dim == 3)
1701 {
1703
1704 // the derivatives are ordered in rows by increasing z
1705 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1706 // j walks through the columns
1707 for (unsigned int d = 0; d < 2; ++d)
1708 for (unsigned int e = 0; e < 2; ++e)
1709 for (unsigned int f = 0; f < 2; ++f)
1710 {
1711 result[i][d][e][f][2] = compressed[5 + d + e + f][i];
1712 result[i][d][e][2][f] = compressed[5 + d + e + f][i];
1713 result[i][d][2][e][f] = compressed[5 + d + e + f][i];
1714 result[i][2][d][e][f] = compressed[5 + d + e + f][i];
1715 }
1716 for (unsigned int d = 0; d < 2; ++d)
1717 for (unsigned int e = 0; e < 2; ++e)
1718 {
1719 result[i][d][e][2][2] = compressed[9 + d + e][i];
1720 result[i][d][2][e][2] = compressed[9 + d + e][i];
1721 result[i][d][2][2][e] = compressed[9 + d + e][i];
1722 result[i][2][d][e][2] = compressed[9 + d + e][i];
1723 result[i][2][d][2][e] = compressed[9 + d + e][i];
1724 result[i][2][2][d][e] = compressed[9 + d + e][i];
1725 }
1726 for (unsigned int d = 0; d < 2; ++d)
1727 {
1728 result[i][d][2][2][2] = compressed[12 + d][i];
1729 result[i][2][d][2][2] = compressed[12 + d][i];
1730 result[i][2][2][d][2] = compressed[12 + d][i];
1731 result[i][2][2][2][d] = compressed[12 + d][i];
1732 }
1733 result[i][2][2][2][2] = compressed[14][i];
1734 }
1735 }
1736 }
1737 return result;
1738 }
1739
1740
1741
1748 template <int dim, int spacedim>
1749 inline void
1751 const CellSimilarity::Similarity cell_similarity,
1752 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1753 const ArrayView<const Point<dim>> &unit_points,
1754 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1755 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1756 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1757 {
1758 if (data.update_each & update_jacobian_3rd_derivatives)
1759 {
1761 data.mapping_support_points);
1762 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1763
1764 if (cell_similarity != CellSimilarity::translation)
1765 {
1766 for (unsigned int point = 0; point < n_q_points; ++point)
1767 {
1768 jacobian_3rd_derivatives[point] = expand_4th_derivatives<dim>(
1769 internal::evaluate_tensor_product_higher_derivatives<4>(
1770 polynomials_1d,
1772 unit_points[point],
1773 renumber_lexicographic_to_hierarchic));
1774 }
1775 }
1776 }
1777 }
1778
1779
1780
1788 template <int dim, int spacedim>
1789 inline void
1791 const CellSimilarity::Similarity cell_similarity,
1792 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1793 const ArrayView<const Point<dim>> &unit_points,
1794 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1795 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1796 std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
1797 {
1799 {
1801 data.mapping_support_points);
1802 const unsigned int n_q_points =
1803 jacobian_pushed_forward_3rd_derivatives.size();
1804
1805 if (cell_similarity != CellSimilarity::translation)
1806 {
1807 ::
1809 tmp;
1811 tmp2;
1812 for (unsigned int point = 0; point < n_q_points; ++point)
1813 {
1816 internal::evaluate_tensor_product_higher_derivatives<4>(
1817 polynomials_1d,
1819 unit_points[point],
1820 renumber_lexicographic_to_hierarchic));
1821
1822 const DerivativeForm<1, dim, spacedim> covariant =
1823 data.output_data->inverse_jacobians[point].transpose();
1824
1825 // push-forward the j-coordinate
1826 for (unsigned int i = 0; i < spacedim; ++i)
1827 for (unsigned int j = 0; j < spacedim; ++j)
1828 for (unsigned int l = 0; l < dim; ++l)
1829 for (unsigned int m = 0; m < dim; ++m)
1830 for (unsigned int n = 0; n < dim; ++n)
1831 {
1832 tmp[i][j][l][m][n] =
1833 fourth[i][0][l][m][n] * covariant[j][0];
1834 for (unsigned int jr = 1; jr < dim; ++jr)
1835 tmp[i][j][l][m][n] +=
1836 fourth[i][jr][l][m][n] * covariant[j][jr];
1837 }
1838
1839 // push-forward the l-coordinate
1840 for (unsigned int i = 0; i < spacedim; ++i)
1841 for (unsigned int j = 0; j < spacedim; ++j)
1842 for (unsigned int l = 0; l < spacedim; ++l)
1843 for (unsigned int m = 0; m < dim; ++m)
1844 for (unsigned int n = 0; n < dim; ++n)
1845 {
1846 tmp2[i][j][l][m][n] =
1847 tmp[i][j][0][m][n] * covariant[l][0];
1848 for (unsigned int lr = 1; lr < dim; ++lr)
1849 tmp2[i][j][l][m][n] +=
1850 tmp[i][j][lr][m][n] * covariant[l][lr];
1851 }
1852
1853 // push-forward the m-coordinate
1854 for (unsigned int i = 0; i < spacedim; ++i)
1855 for (unsigned int j = 0; j < spacedim; ++j)
1856 for (unsigned int l = 0; l < spacedim; ++l)
1857 for (unsigned int m = 0; m < spacedim; ++m)
1858 for (unsigned int n = 0; n < dim; ++n)
1859 {
1860 tmp[i][j][l][m][n] =
1861 tmp2[i][j][l][0][n] * covariant[m][0];
1862 for (unsigned int mr = 1; mr < dim; ++mr)
1863 tmp[i][j][l][m][n] +=
1864 tmp2[i][j][l][mr][n] * covariant[m][mr];
1865 }
1866
1867 // push-forward the n-coordinate
1868 for (unsigned int i = 0; i < spacedim; ++i)
1869 for (unsigned int j = 0; j < spacedim; ++j)
1870 for (unsigned int l = 0; l < spacedim; ++l)
1871 for (unsigned int m = 0; m < spacedim; ++m)
1872 for (unsigned int n = 0; n < spacedim; ++n)
1873 {
1874 jacobian_pushed_forward_3rd_derivatives
1875 [point][i][j][l][m][n] =
1876 tmp[i][j][l][m][0] * covariant[n][0];
1877 for (unsigned int nr = 1; nr < dim; ++nr)
1878 jacobian_pushed_forward_3rd_derivatives[point]
1879 [i][j][l]
1880 [m][n] +=
1881 tmp[i][j][l][m][nr] * covariant[n][nr];
1882 }
1883 }
1884 }
1885 }
1886 }
1887
1888
1889
1899 template <int dim, int spacedim>
1900 inline void
1903 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
1904 const unsigned int face_no,
1905 const unsigned int subface_no,
1906 const unsigned int n_q_points,
1907 const std::vector<double> &weights,
1908 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1910 &output_data)
1911 {
1912 const UpdateFlags update_flags = data.update_each;
1913
1914 if (update_flags &
1916 {
1917 if (update_flags & update_boundary_forms)
1918 AssertDimension(output_data.boundary_forms.size(), n_q_points);
1919 if (update_flags & update_normal_vectors)
1920 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1921 if (update_flags & update_JxW_values)
1922 AssertDimension(output_data.JxW_values.size(), n_q_points);
1923
1924 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1925
1926 // first compute some common data that is used for evaluating
1927 // all of the flags below
1928
1929 // map the unit tangentials to the real cell. checking for d!=dim-1
1930 // eliminates compiler warnings regarding unsigned int expressions <
1931 // 0.
1932 for (unsigned int d = 0; d != dim - 1; ++d)
1933 {
1934 const unsigned int vector_index =
1936 Assert(vector_index < data.unit_tangentials.size(),
1938 Assert(data.aux[d].size() <=
1939 data.unit_tangentials[vector_index].size(),
1941 mapping.transform(make_array_view(
1942 data.unit_tangentials[vector_index]),
1944 data,
1945 make_array_view(data.aux[d]));
1946 }
1947
1948 if (update_flags & update_boundary_forms)
1949 {
1950 // if dim==spacedim, we can use the unit tangentials to compute
1951 // the boundary form by simply taking the cross product
1952 if (dim == spacedim)
1953 {
1954 for (unsigned int i = 0; i < n_q_points; ++i)
1955 switch (dim)
1956 {
1957 case 1:
1958 // in 1d, we don't have access to any of the
1959 // data.aux fields (because it has only dim-1
1960 // components), but we can still compute the
1961 // boundary form by simply looking at the number of
1962 // the face
1963 output_data.boundary_forms[i][0] =
1964 (face_no == 0 ? -1 : +1);
1965 break;
1966 case 2:
1967 output_data.boundary_forms[i] =
1968 cross_product_2d(data.aux[0][i]);
1969 break;
1970 case 3:
1971 output_data.boundary_forms[i] =
1972 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1973 break;
1974 default:
1976 }
1977 }
1978 else //(dim < spacedim)
1979 {
1980 // in the codim-one case, the boundary form results from the
1981 // cross product of all the face tangential vectors and the
1982 // cell normal vector
1983 //
1984 // to compute the cell normal, use the same method used in
1985 // fill_fe_values for cells above
1986 AssertDimension(data.output_data->jacobians.size(),
1987 n_q_points);
1988
1989 for (unsigned int point = 0; point < n_q_points; ++point)
1990 {
1991 const DerivativeForm<1, dim, spacedim> contravariant =
1992 data.output_data->jacobians[point];
1993 if (dim == 1)
1994 {
1995 // J is a tangent vector
1996 output_data.boundary_forms[point] =
1997 contravariant.transpose()[0];
1998 output_data.boundary_forms[point] /=
1999 (face_no == 0 ? -1. : +1.) *
2000 output_data.boundary_forms[point].norm();
2001 }
2002
2003 if (dim == 2)
2004 {
2006 contravariant.transpose();
2007
2009 cross_product_3d(DX_t[0], DX_t[1]);
2011
2012 // then compute the face normal from the face
2013 // tangent and the cell normal:
2014 output_data.boundary_forms[point] =
2015 cross_product_3d(data.aux[0][point], cell_normal);
2016 }
2017 }
2018 }
2019 }
2020
2021 if (update_flags & update_JxW_values)
2022 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
2023 {
2024 output_data.JxW_values[i] =
2025 output_data.boundary_forms[i].norm() * weights[i];
2026
2028 {
2030 cell->subface_case(face_no), subface_no);
2031 output_data.JxW_values[i] *= area_ratio;
2032 }
2033 }
2034
2035 if (update_flags & update_normal_vectors)
2036 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
2037 output_data.normal_vectors[i] =
2038 Point<spacedim>(output_data.boundary_forms[i] /
2039 output_data.boundary_forms[i].norm());
2040 }
2041 }
2042
2043
2044
2051 template <int dim, int spacedim>
2052 inline void
2055 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
2056 const unsigned int face_no,
2057 const unsigned int subface_no,
2059 const Quadrature<dim - 1> &quadrature,
2060 const typename ::MappingQ<dim, spacedim>::InternalData &data,
2061 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
2062 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
2064 &output_data)
2065 {
2066 const ArrayView<const Point<dim>> quadrature_points(
2067 &data.quadrature_points[data_set], quadrature.size());
2068 if (dim > 1 && data.tensor_product_quadrature)
2069 {
2072 data,
2073 output_data.quadrature_points,
2074 output_data.jacobians,
2075 output_data.inverse_jacobians,
2076 output_data.jacobian_grads);
2077 }
2078 else
2079 {
2083 data,
2084 quadrature_points,
2085 polynomials_1d,
2086 renumber_lexicographic_to_hierarchic,
2087 output_data.quadrature_points,
2088 output_data.jacobians,
2089 output_data.inverse_jacobians);
2092 data,
2093 quadrature_points,
2094 polynomials_1d,
2095 renumber_lexicographic_to_hierarchic,
2096 output_data.jacobian_grads);
2097 }
2100 data,
2101 quadrature_points,
2102 polynomials_1d,
2103 renumber_lexicographic_to_hierarchic,
2104 output_data.jacobian_pushed_forward_grads);
2107 data,
2108 quadrature_points,
2109 polynomials_1d,
2110 renumber_lexicographic_to_hierarchic,
2111 output_data.jacobian_2nd_derivatives);
2114 data,
2115 quadrature_points,
2116 polynomials_1d,
2117 renumber_lexicographic_to_hierarchic,
2121 data,
2122 quadrature_points,
2123 polynomials_1d,
2124 renumber_lexicographic_to_hierarchic,
2125 output_data.jacobian_3rd_derivatives);
2128 data,
2129 quadrature_points,
2130 polynomials_1d,
2131 renumber_lexicographic_to_hierarchic,
2133
2135 cell,
2136 face_no,
2137 subface_no,
2138 quadrature.size(),
2139 quadrature.get_weights(),
2140 data,
2141 output_data);
2142 }
2143
2144
2145
2149 template <int dim, int spacedim, int rank>
2150 inline void
2152 const ArrayView<const Tensor<rank, dim>> &input,
2153 const MappingKind mapping_kind,
2154 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2155 const ArrayView<Tensor<rank, spacedim>> &output)
2156 {
2157 AssertDimension(input.size(), output.size());
2158 Assert((dynamic_cast<
2159 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2160 &mapping_data) != nullptr),
2162 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2163 static_cast<
2164 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2165 mapping_data);
2166
2167 switch (mapping_kind)
2168 {
2170 {
2173 "update_contravariant_transformation"));
2174
2175 for (unsigned int i = 0; i < output.size(); ++i)
2176 output[i] = apply_transformation(data.output_data->jacobians[i],
2177 input[i]);
2178
2179 return;
2180 }
2181
2182 case mapping_piola:
2183 {
2186 "update_contravariant_transformation"));
2187 Assert(data.update_each & update_volume_elements,
2189 "update_volume_elements"));
2190 Assert(rank == 1, ExcMessage("Only for rank 1"));
2191 if (rank != 1)
2192 return;
2193
2194 for (unsigned int i = 0; i < output.size(); ++i)
2195 {
2196 output[i] =
2197 apply_transformation(data.output_data->jacobians[i],
2198 input[i]);
2199 output[i] /= data.volume_elements[i];
2200 }
2201 return;
2202 }
2203 // We still allow this operation as in the
2204 // reference cell Derivatives are Tensor
2205 // rather than DerivativeForm
2206 case mapping_covariant:
2207 {
2210 "update_covariant_transformation"));
2211
2212 for (unsigned int i = 0; i < output.size(); ++i)
2213 output[i] = apply_transformation(
2214 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2215
2216 return;
2217 }
2218
2219 default:
2221 }
2222 }
2223
2224
2225
2229 template <int dim, int spacedim, int rank>
2230 inline void
2232 const ArrayView<const Tensor<rank, dim>> &input,
2233 const MappingKind mapping_kind,
2234 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2235 const ArrayView<Tensor<rank, spacedim>> &output)
2236 {
2237 AssertDimension(input.size(), output.size());
2238 Assert((dynamic_cast<
2239 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2240 &mapping_data) != nullptr),
2242 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2243 static_cast<
2244 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2245 mapping_data);
2246
2247 switch (mapping_kind)
2248 {
2250 {
2253 "update_covariant_transformation"));
2256 "update_contravariant_transformation"));
2257 Assert(rank == 2, ExcMessage("Only for rank 2"));
2258
2259 for (unsigned int i = 0; i < output.size(); ++i)
2260 {
2262 apply_transformation(data.output_data->jacobians[i],
2263 transpose(input[i]));
2264 output[i] = apply_transformation(
2265 data.output_data->inverse_jacobians[i].transpose(),
2266 A.transpose());
2267 }
2268
2269 return;
2270 }
2271
2273 {
2276 "update_covariant_transformation"));
2277 Assert(rank == 2, ExcMessage("Only for rank 2"));
2278
2279 for (unsigned int i = 0; i < output.size(); ++i)
2280 {
2281 const DerivativeForm<1, dim, spacedim> covariant =
2282 data.output_data->inverse_jacobians[i].transpose();
2284 apply_transformation(covariant, transpose(input[i]));
2285 output[i] = apply_transformation(covariant, A.transpose());
2286 }
2287
2288 return;
2289 }
2290
2292 {
2295 "update_covariant_transformation"));
2298 "update_contravariant_transformation"));
2299 Assert(data.update_each & update_volume_elements,
2301 "update_volume_elements"));
2302 Assert(rank == 2, ExcMessage("Only for rank 2"));
2303
2304 for (unsigned int i = 0; i < output.size(); ++i)
2305 {
2306 const DerivativeForm<1, dim, spacedim> covariant =
2307 data.output_data->inverse_jacobians[i].transpose();
2309 apply_transformation(covariant, input[i]);
2310 const Tensor<2, spacedim> T =
2311 apply_transformation(data.output_data->jacobians[i],
2312 A.transpose());
2313
2314 output[i] = transpose(T);
2315 output[i] /= data.volume_elements[i];
2316 }
2317
2318 return;
2319 }
2320
2321 default:
2323 }
2324 }
2325
2326
2327
2331 template <int dim, int spacedim>
2332 inline void
2334 const ArrayView<const Tensor<3, dim>> &input,
2335 const MappingKind mapping_kind,
2336 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2337 const ArrayView<Tensor<3, spacedim>> &output)
2338 {
2339 AssertDimension(input.size(), output.size());
2340 Assert((dynamic_cast<
2341 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2342 &mapping_data) != nullptr),
2344 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2345 static_cast<
2346 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2347 mapping_data);
2348
2349 switch (mapping_kind)
2350 {
2352 {
2355 "update_covariant_transformation"));
2358 "update_contravariant_transformation"));
2359
2360 for (unsigned int q = 0; q < output.size(); ++q)
2361 {
2362 const DerivativeForm<1, dim, spacedim> covariant =
2363 data.output_data->inverse_jacobians[q].transpose();
2364 const DerivativeForm<1, dim, spacedim> contravariant =
2365 data.output_data->jacobians[q];
2366
2367 for (unsigned int i = 0; i < spacedim; ++i)
2368 {
2369 double tmp1[dim][dim];
2370 for (unsigned int J = 0; J < dim; ++J)
2371 for (unsigned int K = 0; K < dim; ++K)
2372 {
2373 tmp1[J][K] =
2374 contravariant[i][0] * input[q][0][J][K];
2375 for (unsigned int I = 1; I < dim; ++I)
2376 tmp1[J][K] +=
2377 contravariant[i][I] * input[q][I][J][K];
2378 }
2379 for (unsigned int j = 0; j < spacedim; ++j)
2380 {
2381 double tmp2[dim];
2382 for (unsigned int K = 0; K < dim; ++K)
2383 {
2384 tmp2[K] = covariant[j][0] * tmp1[0][K];
2385 for (unsigned int J = 1; J < dim; ++J)
2386 tmp2[K] += covariant[j][J] * tmp1[J][K];
2387 }
2388 for (unsigned int k = 0; k < spacedim; ++k)
2389 {
2390 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2391 for (unsigned int K = 1; K < dim; ++K)
2392 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2393 }
2394 }
2395 }
2396 }
2397 return;
2398 }
2399
2401 {
2404 "update_covariant_transformation"));
2405
2406 for (unsigned int q = 0; q < output.size(); ++q)
2407 {
2408 const DerivativeForm<1, dim, spacedim> covariant =
2409 data.output_data->inverse_jacobians[q].transpose();
2410
2411 for (unsigned int i = 0; i < spacedim; ++i)
2412 {
2413 double tmp1[dim][dim];
2414 for (unsigned int J = 0; J < dim; ++J)
2415 for (unsigned int K = 0; K < dim; ++K)
2416 {
2417 tmp1[J][K] = covariant[i][0] * input[q][0][J][K];
2418 for (unsigned int I = 1; I < dim; ++I)
2419 tmp1[J][K] += covariant[i][I] * input[q][I][J][K];
2420 }
2421 for (unsigned int j = 0; j < spacedim; ++j)
2422 {
2423 double tmp2[dim];
2424 for (unsigned int K = 0; K < dim; ++K)
2425 {
2426 tmp2[K] = covariant[j][0] * tmp1[0][K];
2427 for (unsigned int J = 1; J < dim; ++J)
2428 tmp2[K] += covariant[j][J] * tmp1[J][K];
2429 }
2430 for (unsigned int k = 0; k < spacedim; ++k)
2431 {
2432 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2433 for (unsigned int K = 1; K < dim; ++K)
2434 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2435 }
2436 }
2437 }
2438 }
2439
2440 return;
2441 }
2442
2444 {
2447 "update_covariant_transformation"));
2450 "update_contravariant_transformation"));
2451 Assert(data.update_each & update_volume_elements,
2453 "update_volume_elements"));
2454
2455 for (unsigned int q = 0; q < output.size(); ++q)
2456 {
2457 const DerivativeForm<1, dim, spacedim> covariant =
2458 data.output_data->inverse_jacobians[q].transpose();
2459 const DerivativeForm<1, dim, spacedim> contravariant =
2460 data.output_data->jacobians[q];
2461 for (unsigned int i = 0; i < spacedim; ++i)
2462 {
2463 double factor[dim];
2464 for (unsigned int I = 0; I < dim; ++I)
2465 factor[I] =
2466 contravariant[i][I] * (1. / data.volume_elements[q]);
2467 double tmp1[dim][dim];
2468 for (unsigned int J = 0; J < dim; ++J)
2469 for (unsigned int K = 0; K < dim; ++K)
2470 {
2471 tmp1[J][K] = factor[0] * input[q][0][J][K];
2472 for (unsigned int I = 1; I < dim; ++I)
2473 tmp1[J][K] += factor[I] * input[q][I][J][K];
2474 }
2475 for (unsigned int j = 0; j < spacedim; ++j)
2476 {
2477 double tmp2[dim];
2478 for (unsigned int K = 0; K < dim; ++K)
2479 {
2480 tmp2[K] = covariant[j][0] * tmp1[0][K];
2481 for (unsigned int J = 1; J < dim; ++J)
2482 tmp2[K] += covariant[j][J] * tmp1[J][K];
2483 }
2484 for (unsigned int k = 0; k < spacedim; ++k)
2485 {
2486 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2487 for (unsigned int K = 1; K < dim; ++K)
2488 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2489 }
2490 }
2491 }
2492 }
2493
2494 return;
2495 }
2496
2497 default:
2499 }
2500 }
2501
2502
2503
2508 template <int dim, int spacedim, int rank>
2509 inline void
2512 const MappingKind mapping_kind,
2513 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2515 {
2516 AssertDimension(input.size(), output.size());
2517 Assert((dynamic_cast<
2518 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2519 &mapping_data) != nullptr),
2521 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2522 static_cast<
2523 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2524 mapping_data);
2525
2526 switch (mapping_kind)
2527 {
2528 case mapping_covariant:
2529 {
2532 "update_covariant_transformation"));
2533
2534 for (unsigned int i = 0; i < output.size(); ++i)
2535 output[i] = apply_transformation(
2536 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2537
2538 return;
2539 }
2540 default:
2542 }
2543 }
2544 } // namespace MappingQImplementation
2545} // namespace internal
2546
2548
2549#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
Abstract base class for mapping classes.
Definition mapping.h:320
numbers::NumberTraits< Number >::real_type norm() const
friend class Tensor
Definition tensor.h:865
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
InverseQuadraticApproximation(const InverseQuadraticApproximation &)=default
Point< dim, Number > compute(const Point< spacedim, Number > &p) const
InverseQuadraticApproximation(const ArrayView< const Point< spacedim > > &real_support_points, const std::vector< Point< dim > > &unit_support_points)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Point< 2 > second
Definition grid_out.cc:4630
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_volume_elements
Determinant of the Jacobian.
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_jacobian_3rd_derivatives
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_quadrature_points
Transformed quadrature points.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
MappingKind
Definition mapping.h:81
@ mapping_piola
Definition mapping.h:116
@ mapping_covariant_gradient
Definition mapping.h:102
@ mapping_covariant
Definition mapping.h:91
@ mapping_contravariant
Definition mapping.h:96
@ mapping_contravariant_hessian
Definition mapping.h:158
@ mapping_covariant_hessian
Definition mapping.h:152
@ mapping_contravariant_gradient
Definition mapping.h:108
@ mapping_piola_gradient
Definition mapping.h:122
@ mapping_piola_hessian
Definition mapping.h:164
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
LogStream deallog
Definition logstream.cc:38
std::vector< index_type > data
Definition mpi.cc:740
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 4, spacedim > > &jacobian_pushed_forward_2nd_derivatives)
Point< dim, Number > do_transform_real_to_unit_cell_internal(const Point< spacedim, Number > &p, const Point< dim, Number > &initial_p_unit, const ArrayView< const Point< spacedim > > &points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber, const bool print_iterations_to_deallog=false)
DerivativeForm< 3, dim, spacedim > expand_3rd_derivatives(const Tensor< 1, length_tensor, Tensor< 1, spacedim > > &compressed)
inline ::Table< 2, double > compute_support_point_weights_cell(const unsigned int polynomial_degree)
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
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)
inline ::Table< 2, double > compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
void transform_fields(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 4, dim, spacedim > > &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 3, spacedim > > &jacobian_pushed_forward_grads)
inline ::Table< 2, double > compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
void maybe_update_q_points_Jacobians_generic(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Point< spacedim > > &quadrature_points, std::vector< DerivativeForm< 1, dim, spacedim > > &jacobians, std::vector< DerivativeForm< 1, spacedim, dim > > &inverse_jacobians)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
void transform_gradients(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void maybe_update_q_points_Jacobians_and_grads_tensor(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Point< spacedim > > &quadrature_points, std::vector< DerivativeForm< 1, dim, spacedim > > &jacobians, std::vector< DerivativeForm< 1, spacedim, dim > > &inverse_jacobians, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
void transform_hessians(const ArrayView< const Tensor< 3, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim > > &output)
void maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 5, spacedim > > &jacobian_pushed_forward_3rd_derivatives)
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)
DerivativeForm< 4, dim, spacedim > expand_4th_derivatives(const Tensor< 1, length_tensor, Tensor< 1, spacedim > > &compressed)
Point< spacedim > compute_mapped_location_of_point(const typename ::MappingQ< dim, spacedim >::InternalData &data)
Point< dim > do_transform_real_to_unit_cell_internal_codim1(const Point< dim+1 > &p, const Point< dim > &initial_p_unit, const ArrayView< const Point< dim+1 > > &points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber)
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)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 3, dim, spacedim > > &jacobian_2nd_derivatives)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)