Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mapping_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
45
47
48namespace internal
49{
55 namespace MappingQ1
56 {
57 // These are left as templates on the spatial dimension (even though dim
58 // == spacedim must be true for them to make sense) because templates are
59 // expanded before the compiler eliminates code due to the 'if (dim ==
60 // spacedim)' statement (see the body of the general
61 // transform_real_to_unit_cell).
62 template <int spacedim>
63 inline Point<1>
66 &vertices,
67 const Point<spacedim> &p)
68 {
69 Assert(spacedim == 1, ExcInternalError());
70 return Point<1>((p[0] - vertices[0][0]) /
71 (vertices[1][0] - vertices[0][0]));
72 }
73
74
75
76 template <int spacedim>
77 inline Point<2>
80 &vertices,
81 const Point<spacedim> &p)
82 {
83 Assert(spacedim == 2, ExcInternalError());
84
85 // For accuracy reasons, we do all arithmetic in extended precision
86 // (long double). This has a noticeable effect on the hit rate for
87 // borderline cases and thus makes the algorithm more robust.
88 const long double x = p[0];
89 const long double y = p[1];
90
91 const long double x0 = vertices[0][0];
92 const long double x1 = vertices[1][0];
93 const long double x2 = vertices[2][0];
94 const long double x3 = vertices[3][0];
95
96 const long double y0 = vertices[0][1];
97 const long double y1 = vertices[1][1];
98 const long double y2 = vertices[2][1];
99 const long double y3 = vertices[3][1];
100
101 const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3);
102 const long double b = -(x0 - x1 - x2 + x3) * y + (x - 2 * x1 + x3) * y0 -
103 (x - 2 * x0 + x2) * y1 - (x - x1) * y2 +
104 (x - x0) * y3;
105 const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1;
106
107 const long double discriminant = b * b - 4 * a * c;
108 // exit if the point is not in the cell (this is the only case where the
109 // discriminant is negative)
111 discriminant > 0.0,
113
114 long double eta1;
115 long double eta2;
116 const long double sqrt_discriminant = std::sqrt(discriminant);
117 // special case #1: if a is near-zero to make the discriminant exactly
118 // equal b, then use the linear formula
119 if (b != 0.0 && std::abs(b) == sqrt_discriminant)
120 {
121 eta1 = -c / b;
122 eta2 = -c / b;
123 }
124 // special case #2: a is zero for parallelograms and very small for
125 // near-parallelograms:
126 else if (std::abs(a) < 1e-8 * std::abs(b))
127 {
128 // if both a and c are very small then the root should be near
129 // zero: this first case will capture that
130 eta1 = 2 * c / (-b - sqrt_discriminant);
131 eta2 = 2 * c / (-b + sqrt_discriminant);
132 }
133 // finally, use the plain version:
134 else
135 {
136 eta1 = (-b - sqrt_discriminant) / (2 * a);
137 eta2 = (-b + sqrt_discriminant) / (2 * a);
138 }
139 // pick the one closer to the center of the cell.
140 const long double eta =
141 (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
142
143 /*
144 * There are two ways to compute xi from eta, but either one may have a
145 * zero denominator.
146 */
147 const long double subexpr0 = -eta * x2 + x0 * (eta - 1);
148 const long double xi_denominator0 = eta * x3 - x1 * (eta - 1) + subexpr0;
149 const long double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
150 std::max(std::abs(x2), std::abs(x3)));
151
152 if (std::abs(xi_denominator0) > 1e-10 * max_x)
153 {
154 const double xi = (x + subexpr0) / xi_denominator0;
155 return {xi, static_cast<double>(eta)};
156 }
157 else
158 {
159 const long double max_y =
161 std::max(std::abs(y2), std::abs(y3)));
162 const long double subexpr1 = -eta * y2 + y0 * (eta - 1);
163 const long double xi_denominator1 =
164 eta * y3 - y1 * (eta - 1) + subexpr1;
165 if (std::abs(xi_denominator1) > 1e-10 * max_y)
166 {
167 const double xi = (subexpr1 + y) / xi_denominator1;
168 return {xi, static_cast<double>(eta)};
169 }
170 else // give up and try Newton iteration
171 {
173 false,
174 (typename Mapping<spacedim,
175 spacedim>::ExcTransformationFailed()));
176 }
177 }
178 // bogus return to placate compiler. It should not be possible to get
179 // here.
181 return {std::numeric_limits<double>::quiet_NaN(),
182 std::numeric_limits<double>::quiet_NaN()};
183 }
184
185
186
187 template <int spacedim>
188 inline Point<3>
191 & /*vertices*/,
192 const Point<spacedim> & /*p*/)
193 {
194 // It should not be possible to get here
196 return {std::numeric_limits<double>::quiet_NaN(),
197 std::numeric_limits<double>::quiet_NaN(),
198 std::numeric_limits<double>::quiet_NaN()};
199 }
200 } // namespace MappingQ1
201
202
203
209 namespace MappingQImplementation
210 {
215 template <int dim>
216 std::vector<Point<dim>>
217 unit_support_points(const std::vector<Point<1>> &line_support_points,
218 const std::vector<unsigned int> &renumbering)
219 {
220 AssertDimension(Utilities::pow(line_support_points.size(), dim),
221 renumbering.size());
222 std::vector<Point<dim>> points(renumbering.size());
223 const unsigned int n1 = line_support_points.size();
224 for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2)
225 for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1)
226 for (unsigned int q0 = 0; q0 < n1; ++q0, ++q)
227 {
228 points[renumbering[q]][0] = line_support_points[q0][0];
229 if (dim > 1)
230 points[renumbering[q]][1] = line_support_points[q1][0];
231 if (dim > 2)
232 points[renumbering[q]][2] = line_support_points[q2][0];
233 }
234 return points;
235 }
236
237
238
246 inline ::Table<2, double>
247 compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
248 {
249 ::Table<2, double> loqvs;
250
251 // we are asked to compute weights for interior support points, but
252 // there are no interior points if degree==1
253 if (polynomial_degree == 1)
254 return loqvs;
255
256 const unsigned int M = polynomial_degree - 1;
257 const unsigned int n_inner_2d = M * M;
258 const unsigned int n_outer_2d = 4 + 4 * M;
259
260 // set the weights of transfinite interpolation
261 loqvs.reinit(n_inner_2d, n_outer_2d);
262 QGaussLobatto<2> gl(polynomial_degree + 1);
263 for (unsigned int i = 0; i < M; ++i)
264 for (unsigned int j = 0; j < M; ++j)
265 {
266 const Point<2> &p =
267 gl.point((i + 1) * (polynomial_degree + 1) + (j + 1));
268 const unsigned int index_table = i * M + j;
269 for (unsigned int v = 0; v < 4; ++v)
270 loqvs(index_table, v) =
272 loqvs(index_table, 4 + i) = 1. - p[0];
273 loqvs(index_table, 4 + i + M) = p[0];
274 loqvs(index_table, 4 + j + 2 * M) = 1. - p[1];
275 loqvs(index_table, 4 + j + 3 * M) = p[1];
276 }
277
278 // the sum of weights of the points at the outer rim should be one.
279 // check this
280 for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
281 Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
282 loqvs[unit_point].end(),
283 0.) -
284 1) < 1e-13 * polynomial_degree,
286
287 return loqvs;
288 }
289
290
291
298 inline ::Table<2, double>
299 compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
300 {
301 ::Table<2, double> lohvs;
302
303 // we are asked to compute weights for interior support points, but
304 // there are no interior points if degree==1
305 if (polynomial_degree == 1)
306 return lohvs;
307
308 const unsigned int M = polynomial_degree - 1;
309
310 const unsigned int n_inner = Utilities::fixed_power<3>(M);
311 const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
312
313 // set the weights of transfinite interpolation
314 lohvs.reinit(n_inner, n_outer);
315 QGaussLobatto<3> gl(polynomial_degree + 1);
316 for (unsigned int i = 0; i < M; ++i)
317 for (unsigned int j = 0; j < M; ++j)
318 for (unsigned int k = 0; k < M; ++k)
319 {
320 const Point<3> &p = gl.point((i + 1) * (M + 2) * (M + 2) +
321 (j + 1) * (M + 2) + (k + 1));
322 const unsigned int index_table = i * M * M + j * M + k;
323
324 // vertices
325 for (unsigned int v = 0; v < 8; ++v)
326 lohvs(index_table, v) =
328
329 // lines
330 {
331 constexpr std::array<unsigned int, 4> line_coordinates_y(
332 {{0, 1, 4, 5}});
333 const Point<2> py(p[0], p[2]);
334 for (unsigned int l = 0; l < 4; ++l)
335 lohvs(index_table, 8 + line_coordinates_y[l] * M + j) =
337 }
338
339 {
340 constexpr std::array<unsigned int, 4> line_coordinates_x(
341 {{2, 3, 6, 7}});
342 const Point<2> px(p[1], p[2]);
343 for (unsigned int l = 0; l < 4; ++l)
344 lohvs(index_table, 8 + line_coordinates_x[l] * M + k) =
346 }
347
348 {
349 constexpr std::array<unsigned int, 4> line_coordinates_z(
350 {{8, 9, 10, 11}});
351 const Point<2> pz(p[0], p[1]);
352 for (unsigned int l = 0; l < 4; ++l)
353 lohvs(index_table, 8 + line_coordinates_z[l] * M + i) =
355 }
356
357 // quads
358 lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
359 1. - p[0];
360 lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0];
361 lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) =
362 1. - p[1];
363 lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1];
364 lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) =
365 1. - p[2];
366 lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
367 }
368
369 // the sum of weights of the points at the outer rim should be one.
370 // check this
371 for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
372 Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
373 lohvs[unit_point].end(),
374 0.) -
375 1) < 1e-13 * polynomial_degree,
377
378 return lohvs;
379 }
380
381
382
387 inline std::vector<::Table<2, double>>
389 const unsigned int polynomial_degree,
390 const unsigned int dim)
391 {
392 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
393 std::vector<::Table<2, double>> output(dim);
394 if (polynomial_degree <= 1)
395 return output;
396
397 // fill the 1d interior weights
398 QGaussLobatto<1> quadrature(polynomial_degree + 1);
399 output[0].reinit(polynomial_degree - 1,
401 for (unsigned int q = 0; q < polynomial_degree - 1; ++q)
402 for (const unsigned int i : GeometryInfo<1>::vertex_indices())
403 output[0](q, i) =
404 GeometryInfo<1>::d_linear_shape_function(quadrature.point(q + 1),
405 i);
406
407 if (dim > 1)
408 output[1] = compute_support_point_weights_on_quad(polynomial_degree);
409
410 if (dim > 2)
411 output[2] = compute_support_point_weights_on_hex(polynomial_degree);
412
413 return output;
414 }
415
416
417
421 template <int dim>
422 inline ::Table<2, double>
423 compute_support_point_weights_cell(const unsigned int polynomial_degree)
424 {
425 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
426 if (polynomial_degree <= 1)
427 return ::Table<2, double>();
428
429 QGaussLobatto<dim> quadrature(polynomial_degree + 1);
430 const std::vector<unsigned int> h2l =
432
433 ::Table<2, double> output(quadrature.size() -
436 for (unsigned int q = 0; q < output.size(0); ++q)
437 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
439 quadrature.point(h2l[q + GeometryInfo<dim>::vertices_per_cell]), i);
440
441 return output;
442 }
443
444
445
453 template <int dim, int spacedim>
454 inline Point<spacedim>
456 const typename ::MappingQ<dim, spacedim>::InternalData &data)
457 {
458 AssertDimension(data.shape_values.size(),
459 data.mapping_support_points.size());
460
461 // use now the InternalData to compute the point in real space.
462 Point<spacedim> p_real;
463 for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i)
464 p_real += data.mapping_support_points[i] * data.shape(0, i);
465
466 return p_real;
467 }
468
469
470
475 template <int dim, int spacedim, typename Number>
476 inline Point<dim, Number>
479 const Point<dim, Number> &initial_p_unit,
480 const ArrayView<const Point<spacedim>> &points,
481 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
482 const std::vector<unsigned int> &renumber,
483 const bool print_iterations_to_deallog = false)
484 {
485 if (print_iterations_to_deallog)
486 deallog << "Start MappingQ::do_transform_real_to_unit_cell for real "
487 << "point [ " << p << " ] " << std::endl;
488
489 AssertDimension(points.size(),
490 Utilities::pow(polynomials_1d.size(), dim));
491
492 // Newton iteration to solve
493 // f(x)=p(x)-p=0
494 // where we are looking for 'x' and p(x) is the forward transformation
495 // from unit to real cell. We solve this using a Newton iteration
496 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
497 // The start value is set to be the linear approximation to the cell
498
499 // The shape values and derivatives of the mapping at this point are
500 // previously computed.
501
502 Point<dim, Number> p_unit = initial_p_unit;
504 polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber);
505
506 Tensor<1, spacedim, Number> f = p_real.first - p;
507
508 // early out if we already have our point in all SIMD lanes, i.e.,
509 // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable
510 // this step also for VectorizedArray where we do not have operator <,
511 // compare the result to zero which is defined for SIMD types
512 if (std::max(Number(0.),
513 f.norm_square() - 1e-24 * p_real.second[0].norm_square()) ==
514 Number(0.))
515 return p_unit;
516
517 // we need to compare the position of the computed p(x) against the
518 // given point 'p'. We will terminate the iteration and return 'x' if
519 // they are less than eps apart. The question is how to choose eps --
520 // or, put maybe more generally: in which norm we want these 'p' and
521 // 'p(x)' to be eps apart.
522 //
523 // the question is difficult since we may have to deal with very
524 // elongated cells where we may achieve 1e-12*h for the distance of
525 // these two points in the 'long' direction, but achieving this
526 // tolerance in the 'short' direction of the cell may not be possible
527 //
528 // what we do instead is then to terminate iterations if
529 // \| p(x) - p \|_A < eps
530 // where the A-norm is somehow induced by the transformation of the
531 // cell. in particular, we want to measure distances relative to the
532 // sizes of the cell in its principal directions.
533 //
534 // to define what exactly A should be, note that to first order we have
535 // the following (assuming that x* is the solution of the problem, i.e.,
536 // p(x*)=p):
537 // p(x) - p = p(x) - p(x*)
538 // = -grad p(x) * (x*-x) + higher order terms
539 // This suggest to measure with a norm that corresponds to
540 // A = {[grad p(x)]^T [grad p(x)]}^{-1}
541 // because then
542 // \| p(x) - p \|_A \approx \| x - x* \|
543 // Consequently, we will try to enforce that
544 // \| p(x) - p \|_A = \| f \| <= eps
545 //
546 // Note that using this norm is a bit dangerous since the norm changes
547 // in every iteration (A isn't fixed by depending on xk). However, if
548 // the cell is not too deformed (it may be stretched, but not twisted)
549 // then the mapping is almost linear and A is indeed constant or
550 // nearly so.
551 const double eps = 1.e-11;
552 const unsigned int newton_iteration_limit = 20;
553
554 Point<dim, Number> invalid_point;
555 invalid_point[0] = std::numeric_limits<double>::infinity();
556 bool tried_project_to_unit_cell = false;
557
558 unsigned int newton_iteration = 0;
559 Number f_weighted_norm_square = 1.;
560 Number last_f_weighted_norm_square = 1.;
561
562 do
563 {
564 if (print_iterations_to_deallog)
565 deallog << "Newton iteration " << newton_iteration
566 << " for unit point guess " << p_unit << std::endl;
567
568 // f'(x)
570 for (unsigned int d = 0; d < spacedim; ++d)
571 for (unsigned int e = 0; e < dim; ++e)
572 df[d][e] = p_real.second[e][d];
573
574 // check determinand(df) > 0 on all SIMD lanes
575 if (!(std::min(determinant(df),
576 Number(std::numeric_limits<double>::min())) ==
577 Number(std::numeric_limits<double>::min())))
578 {
579 // We allow to enter this function with an initial guess
580 // outside the unit cell. We might have run into invalid
581 // Jacobians due to that, so we should at least try once to go
582 // back to the unit cell and go on with the Newton iteration
583 // from there. Since the outside case is unlikely, we can
584 // afford spending the extra effort at this place.
585 if (tried_project_to_unit_cell == false)
586 {
589 polynomials_1d,
590 points,
591 p_unit,
592 polynomials_1d.size() == 2,
593 renumber);
594 f = p_real.first - p;
595 f_weighted_norm_square = 1.;
596 last_f_weighted_norm_square = 1;
597 tried_project_to_unit_cell = true;
598 continue;
599 }
600 else
601 return invalid_point;
602 }
603
604 // Solve [f'(x)]d=f(x)
605 const Tensor<2, spacedim, Number> df_inverse = invert(df);
606 const Tensor<1, spacedim, Number> delta = df_inverse * f;
607 last_f_weighted_norm_square = delta.norm_square();
608
609 if (print_iterations_to_deallog)
610 deallog << " delta=" << delta << std::endl;
611
612 // do a line search
613 double step_length = 1.0;
614 do
615 {
616 // update of p_unit. The spacedim-th component of transformed
617 // point is simply ignored in codimension one case. When this
618 // component is not zero, then we are projecting the point to
619 // the surface or curve identified by the cell.
620 Point<dim, Number> p_unit_trial = p_unit;
621 for (unsigned int i = 0; i < dim; ++i)
622 p_unit_trial[i] -= step_length * delta[i];
623
624 // shape values and derivatives at new p_unit point
625 const auto p_real_trial =
627 polynomials_1d,
628 points,
629 p_unit_trial,
630 polynomials_1d.size() == 2,
631 renumber);
632 const Tensor<1, spacedim, Number> f_trial =
633 p_real_trial.first - p;
634 f_weighted_norm_square = (df_inverse * f_trial).norm_square();
635
636 if (print_iterations_to_deallog)
637 {
638 deallog << " step_length=" << step_length << std::endl;
639 if (step_length == 1.0)
640 deallog << " ||f || =" << f.norm() << std::endl;
641 deallog << " ||f*|| =" << f_trial.norm() << std::endl
642 << " ||f*||_A ="
643 << std::sqrt(f_weighted_norm_square) << std::endl;
644 }
645
646 // See if we are making progress with the current step length
647 // and if not, reduce it by a factor of two and try again.
648 //
649 // Strictly speaking, we should probably use the same norm as we
650 // use for the outer algorithm. In practice, line search is just
651 // a crutch to find a "reasonable" step length, and so using the
652 // l2 norm is probably just fine.
653 //
654 // check f_trial.norm() < f.norm() in SIMD form. This is a bit
655 // more complicated because some SIMD lanes might not be doing
656 // any progress any more as they have already reached roundoff
657 // accuracy. We define that as the case of updates less than
658 // 1e-6. The tolerance might seem coarse but since we are
659 // dealing with a Newton iteration of a polynomial function we
660 // either converge quadratically or not any more. Thus, our
661 // selection is to terminate if either the norm of f is
662 // decreasing or that threshold of 1e-6 is reached.
663 if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) *
664 std::max(f_trial.norm_square() - f.norm_square(),
665 Number(0.)) ==
666 Number(0.))
667 {
668 p_real = p_real_trial;
669 p_unit = p_unit_trial;
670 f = f_trial;
671 break;
672 }
673 else if (step_length > 0.05)
674 step_length *= 0.5;
675 else
676 break;
677 }
678 while (true);
679
680 // In case we terminated the line search due to the step becoming
681 // too small, we give the iteration another try with the
682 // projection of the initial guess to the unit cell before we give
683 // up, just like for the negative determinant case.
684 if (step_length <= 0.05 && tried_project_to_unit_cell == false)
685 {
688 polynomials_1d,
689 points,
690 p_unit,
691 polynomials_1d.size() == 2,
692 renumber);
693 f = p_real.first - p;
694 f_weighted_norm_square = 1.;
695 last_f_weighted_norm_square = 1;
696 tried_project_to_unit_cell = true;
697 continue;
698 }
699 else if (step_length <= 0.05)
700 return invalid_point;
701
702 ++newton_iteration;
703 if (newton_iteration > newton_iteration_limit)
704 return invalid_point;
705 }
706 // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the
707 // weighted norm is less than 1e-6 and the improvement against the
708 // previous step was less than a factor of 10 (in that regime, we
709 // either have quadratic convergence or roundoff errors due to a bad
710 // mapping)
711 while (
712 !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) *
713 std::max(last_f_weighted_norm_square -
714 std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) *
715 100.,
716 Number(0.)) ==
717 Number(0.)));
718
719 if (print_iterations_to_deallog)
720 deallog << "Iteration converged for p_unit = [ " << p_unit
721 << " ] and iteration error "
722 << std::sqrt(f_weighted_norm_square) << std::endl;
723
724 return p_unit;
725 }
726
727
728
732 template <int dim>
733 inline Point<dim>
735 const Point<dim + 1> &p,
736 const Point<dim> &initial_p_unit,
737 const ArrayView<const Point<dim + 1>> &points,
738 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
739 const std::vector<unsigned int> &renumber)
740 {
741 const unsigned int spacedim = dim + 1;
742
743 AssertDimension(points.size(),
744 Utilities::pow(polynomials_1d.size(), dim));
745
746 Point<dim> p_unit = initial_p_unit;
747
748 const double eps = 1.e-12;
749 const unsigned int loop_limit = 10;
750
751 unsigned int loop = 0;
752 double f_weighted_norm_square = 1.;
753
754 while (f_weighted_norm_square > eps * eps && loop++ < loop_limit)
755 {
756 const auto p_real =
758 polynomials_1d,
759 points,
760 p_unit,
761 polynomials_1d.size() == 2,
762 renumber);
763 Tensor<1, spacedim> p_minus_F = p - p_real.first;
764 const DerivativeForm<1, spacedim, dim> DF = p_real.second;
765
767 polynomials_1d, points, p_unit, renumber);
768 Point<dim> f;
770 for (unsigned int j = 0; j < dim; ++j)
771 {
772 f[j] = DF[j] * p_minus_F;
773 for (unsigned int l = 0; l < dim; ++l)
774 df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F;
775 }
776
777 // Solve [df(x)]d=f(x)
778 const Tensor<1, dim> d =
779 invert(df) * static_cast<const Tensor<1, dim> &>(f);
780 f_weighted_norm_square = d.norm_square();
781 p_unit -= d;
782 }
783
784
785 // Here we check that in the last execution of while the first
786 // condition was already wrong, meaning the residual was below
787 // eps. Only if the first condition failed, loop will have been
788 // increased and tested, and thus have reached the limit.
789 AssertThrow(loop < loop_limit,
791
792 return p_unit;
793 }
794
795
796
814 template <int dim, int spacedim>
815 class InverseQuadraticApproximation
816 {
817 public:
821 static constexpr unsigned int n_functions =
822 (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
823
836 const ArrayView<const Point<spacedim>> &real_support_points,
837 const std::vector<Point<dim>> &unit_support_points)
838 : normalization_shift(real_support_points[0])
840 1. / real_support_points[0].distance(real_support_points[1]))
841 , is_affine(true)
842 {
843 AssertDimension(real_support_points.size(), unit_support_points.size());
844
845 // For the bi-/trilinear approximation, we cannot build a quadratic
846 // polynomial due to a lack of points (interpolation matrix would
847 // get singular). Similarly, it is not entirely clear how to gather
848 // enough information for the case dim < spacedim.
849 //
850 // In both cases we require the vector real_support_points to
851 // contain the vertex positions and fall back to an affine
852 // approximation:
853 Assert(dim == spacedim || real_support_points.size() ==
856 if (real_support_points.size() == GeometryInfo<dim>::vertices_per_cell)
857 {
859 make_array_view(real_support_points));
861 affine.first.covariant_form().transpose();
862
863 // The code for evaluation assumes an additional transformation of
864 // the form (x - normalization_shift) * normalization_length --
865 // account for this in the definition of the coefficients.
866 coefficients[0] =
867 apply_transformation(A_inv, normalization_shift - affine.second);
868 for (unsigned int d = 0; d < spacedim; ++d)
869 for (unsigned int e = 0; e < dim; ++e)
870 coefficients[1 + d][e] =
871 A_inv[e][d] * (1.0 / normalization_length);
872 is_affine = true;
873 return;
874 }
875
877 std::array<double, n_functions> shape_values;
878 for (unsigned int q = 0; q < unit_support_points.size(); ++q)
879 {
880 // Evaluate quadratic shape functions in point, with the
881 // normalization applied in order to avoid roundoff issues with
882 // scaling far away from 1.
883 shape_values[0] = 1.;
884 const Tensor<1, spacedim> p_scaled =
886 (real_support_points[q] - normalization_shift);
887 for (unsigned int d = 0; d < spacedim; ++d)
888 shape_values[1 + d] = p_scaled[d];
889 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
890 for (unsigned int e = 0; e <= d; ++e, ++c)
891 shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
892
893 // Build lower diagonal of least squares matrix and rhs, the
894 // essential part being that we construct the matrix with the
895 // real points and the right hand side by comparing to the
896 // reference point positions which sets up an inverse
897 // interpolation.
898 for (unsigned int i = 0; i < n_functions; ++i)
899 for (unsigned int j = 0; j < n_functions; ++j)
900 matrix[i][j] += shape_values[i] * shape_values[j];
901 for (unsigned int i = 0; i < n_functions; ++i)
902 coefficients[i] += shape_values[i] * unit_support_points[q];
903 }
904
905 // Factorize the matrix A = L * L^T in-place with the
906 // Cholesky-Banachiewicz algorithm. The implementation is similar to
907 // FullMatrix::cholesky() but re-implemented to avoid memory
908 // allocations and some unnecessary divisions which we can do here as
909 // we only need to solve with dim right hand sides.
910 for (unsigned int i = 0; i < n_functions; ++i)
911 {
912 double Lij_sum = 0;
913 for (unsigned int j = 0; j < i; ++j)
914 {
915 double Lik_Ljk_sum = 0;
916 for (unsigned int k = 0; k < j; ++k)
917 Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
918 matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
919 Lij_sum += matrix[i][j] * matrix[i][j];
920 }
921 AssertThrow(matrix[i][i] - Lij_sum >= 0,
922 ExcMessage("Matrix of normal equations not positive "
923 "definite"));
924
925 // Store the inverse in the diagonal since that is the quantity
926 // needed later in the factorization as well as the forward and
927 // backward substitution, minimizing the number of divisions.
928 matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
929 }
930
931 // Solve lower triangular part, L * y = rhs.
932 for (unsigned int i = 0; i < n_functions; ++i)
933 {
934 Point<dim> sum = coefficients[i];
935 for (unsigned int j = 0; j < i; ++j)
936 sum -= matrix[i][j] * coefficients[j];
937 coefficients[i] = sum * matrix[i][i];
938 }
939
940 // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
941 for (unsigned int i = n_functions; i > 0;)
942 {
943 --i;
944 Point<dim> sum = coefficients[i];
945 for (unsigned int j = i + 1; j < n_functions; ++j)
946 sum -= matrix[j][i] * coefficients[j];
947 coefficients[i] = sum * matrix[i][i];
948 }
949
950 // Check whether the approximation is indeed affine, allowing to
951 // skip the quadratic terms.
952 is_affine = true;
953 for (unsigned int i = dim + 1; i < n_functions; ++i)
954 if (coefficients[i].norm_square() > 1e-20)
955 {
956 is_affine = false;
957 break;
958 }
959 }
960
965 default;
966
970 template <typename Number>
973 {
974 Point<dim, Number> result;
975 for (unsigned int d = 0; d < dim; ++d)
976 result[d] = coefficients[0][d];
977
978 // Apply the normalization to ensure a good conditioning. Since Number
979 // might be a vectorized array whereas the normalization is a point of
980 // doubles, we cannot use the overload of operator- and must instead
981 // loop over the components of the point.
983 for (unsigned int d = 0; d < spacedim; ++d)
984 p_scaled[d] = (p[d] - normalization_shift[d]) * normalization_length;
985
986 for (unsigned int d = 0; d < spacedim; ++d)
987 result += coefficients[1 + d] * p_scaled[d];
988
989 if (!is_affine)
990 {
991 Point<dim, Number> result_affine = result;
992 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
993 for (unsigned int e = 0; e <= d; ++e, ++c)
994 result +=
995 coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
996
997 // Check if the quadratic approximation ends up considerably
998 // farther outside the unit cell on some or all SIMD lanes than
999 // the affine approximation - in that case, we switch those
1000 // components back to the affine approximation. Note that the
1001 // quadratic approximation will grow more quickly away from the
1002 // unit cell. We make the selection for each SIMD lane with a
1003 // ternary operation.
1004 const Number distance_to_unit_cell = result.distance_square(
1006 const Number affine_distance_to_unit_cell =
1007 result_affine.distance_square(
1009 for (unsigned int d = 0; d < dim; ++d)
1011 distance_to_unit_cell,
1012 affine_distance_to_unit_cell + 0.5,
1013 result_affine[d],
1014 result[d]);
1015 }
1016 return result;
1017 }
1018
1019 private:
1028
1033
1037 std::array<Point<dim>, n_functions> coefficients;
1038
1045 };
1046
1047
1048
1054 template <int dim, int spacedim>
1055 inline void
1057 const CellSimilarity::Similarity cell_similarity,
1058 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1059 std::vector<Point<spacedim>> &quadrature_points,
1060 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1061 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians,
1062 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1063 {
1064 const UpdateFlags update_flags = data.update_each;
1065
1066 using VectorizedArrayType =
1067 typename ::MappingQ<dim,
1068 spacedim>::InternalData::VectorizedArrayType;
1069 const unsigned int n_shape_values = data.n_shape_functions;
1070 const unsigned int n_q_points = data.shape_info.n_q_points;
1071 constexpr unsigned int n_lanes = VectorizedArrayType::size();
1072 constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1073 constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1074
1075 // Since MappingQ::InternalData does not have separate arrays for the
1076 // covariant and contravariant transformations, but uses the arrays in
1077 // the `MappingRelatedData`, it can happen that vectors do not have the
1078 // right size
1079 if (update_flags & update_contravariant_transformation)
1080 jacobians.resize(n_q_points);
1081 if (update_flags & update_covariant_transformation)
1082 inverse_jacobians.resize(n_q_points);
1083
1084 EvaluationFlags::EvaluationFlags evaluation_flag =
1087 ((cell_similarity != CellSimilarity::translation) &&
1088 (update_flags & update_contravariant_transformation) ?
1091 ((cell_similarity != CellSimilarity::translation) &&
1092 (update_flags & update_jacobian_grads) ?
1095
1096 Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1098 Assert(!(evaluation_flag & EvaluationFlags::values) ||
1099 n_q_points == quadrature_points.size(),
1100 ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1101 Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1102 data.n_shape_functions > 0,
1104 Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
1105 n_q_points == jacobian_grads.size(),
1106 ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1107
1108 // shortcut in case we have an identity interpolation and only request
1109 // the quadrature points
1110 if (evaluation_flag == EvaluationFlags::values &&
1111 data.shape_info.element_type ==
1113 {
1114 for (unsigned int q = 0; q < n_q_points; ++q)
1115 quadrature_points[q] =
1116 data.mapping_support_points[data.shape_info
1117 .lexicographic_numbering[q]];
1118 return;
1119 }
1120
1122
1123 // prepare arrays
1124 if (evaluation_flag != EvaluationFlags::nothing)
1125 {
1126 eval.set_data_pointers(&data.scratch, n_comp);
1127
1128 // make sure to initialize on all lanes also when some are unused in
1129 // the code below
1130 for (unsigned int i = 0; i < n_shape_values * n_comp; ++i)
1131 eval.begin_dof_values()[i] = VectorizedArrayType();
1132
1133 const std::vector<unsigned int> &renumber_to_lexicographic =
1134 data.shape_info.lexicographic_numbering;
1135 for (unsigned int i = 0; i < n_shape_values; ++i)
1136 for (unsigned int d = 0; d < spacedim; ++d)
1137 {
1138 const unsigned int in_comp = d % n_lanes;
1139 const unsigned int out_comp = d / n_lanes;
1140 eval
1141 .begin_dof_values()[out_comp * n_shape_values + i][in_comp] =
1142 data.mapping_support_points[renumber_to_lexicographic[i]][d];
1143 }
1144
1145 // do the actual tensorized evaluation
1147 n_comp, evaluation_flag, eval.begin_dof_values(), eval);
1148 }
1149
1150 // do the postprocessing
1151 if (evaluation_flag & EvaluationFlags::values)
1152 {
1153 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1154 for (unsigned int i = 0; i < n_q_points; ++i)
1155 for (unsigned int in_comp = 0;
1156 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1157 ++in_comp)
1158 quadrature_points[i][out_comp * n_lanes + in_comp] =
1159 eval.begin_values()[out_comp * n_q_points + i][in_comp];
1160 }
1161
1162 if (evaluation_flag & EvaluationFlags::gradients)
1163 {
1164 // We need to reinterpret the data after evaluate has been applied.
1165 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1166 for (unsigned int point = 0; point < n_q_points; ++point)
1167 for (unsigned int in_comp = 0;
1168 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1169 ++in_comp)
1170 for (unsigned int j = 0; j < dim; ++j)
1171 {
1172 jacobians[point][out_comp * n_lanes + in_comp][j] =
1173 eval.begin_gradients()[(out_comp * n_q_points + point) *
1174 dim +
1175 j][in_comp];
1176 }
1177 }
1178
1179 if (update_flags & update_volume_elements)
1180 if (cell_similarity != CellSimilarity::translation)
1181 for (unsigned int point = 0; point < n_q_points; ++point)
1182 data.volume_elements[point] = jacobians[point].determinant();
1183
1184 // copy values from InternalData to vector given by reference
1185 if (update_flags & update_covariant_transformation)
1186 {
1187 AssertDimension(jacobians.size(), n_q_points);
1188 AssertDimension(inverse_jacobians.size(), n_q_points);
1189 if (cell_similarity != CellSimilarity::translation)
1190 for (unsigned int point = 0; point < n_q_points; ++point)
1191 inverse_jacobians[point] =
1192 jacobians[point].covariant_form().transpose();
1193 }
1194
1195 if (evaluation_flag & EvaluationFlags::hessians)
1196 {
1197 constexpr int desymmetrize_3d[6][2] = {
1198 {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1199 constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1200
1201 // We need to reinterpret the data after evaluate has been applied.
1202 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1203 for (unsigned int point = 0; point < n_q_points; ++point)
1204 for (unsigned int j = 0; j < n_hessians; ++j)
1205 for (unsigned int in_comp = 0;
1206 in_comp < n_lanes &&
1207 in_comp < spacedim - out_comp * n_lanes;
1208 ++in_comp)
1209 {
1210 const unsigned int total_number = point * n_hessians + j;
1211 const unsigned int new_point = total_number % n_q_points;
1212 const unsigned int new_hessian_comp =
1213 total_number / n_q_points;
1214 const unsigned int new_hessian_comp_i =
1215 dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1216 desymmetrize_3d[new_hessian_comp][0];
1217 const unsigned int new_hessian_comp_j =
1218 dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1219 desymmetrize_3d[new_hessian_comp][1];
1220 const double value =
1221 eval.begin_hessians()[(out_comp * n_q_points + point) *
1222 n_hessians +
1223 j][in_comp];
1224 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1225 [new_hessian_comp_i][new_hessian_comp_j] =
1226 value;
1227 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1228 [new_hessian_comp_j][new_hessian_comp_i] =
1229 value;
1230 }
1231 }
1232 }
1233
1234
1235
1236 template <int dim, int spacedim>
1237 inline void
1239 const CellSimilarity::Similarity cell_similarity,
1240 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1241 const ArrayView<const Point<dim>> &unit_points,
1242 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1243 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1244 std::vector<Point<spacedim>> &quadrature_points,
1245 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1246 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
1247 {
1248 const UpdateFlags update_flags = data.update_each;
1249 const ArrayView<const Point<spacedim>> support_points(
1250 data.mapping_support_points);
1251
1252 const unsigned int n_points = unit_points.size();
1253 const unsigned int n_lanes = VectorizedArray<double>::size();
1254
1255 // Since MappingQ::InternalData does not have separate arrays for the
1256 // covariant and contravariant transformations, but uses the arrays in
1257 // the `MappingRelatedData`, it can happen that vectors do not have the
1258 // right size
1259 if (update_flags & update_contravariant_transformation)
1260 jacobians.resize(n_points);
1261 if (update_flags & update_covariant_transformation)
1262 inverse_jacobians.resize(n_points);
1263
1264 const bool is_translation =
1265 cell_similarity == CellSimilarity::translation;
1266
1267 const bool needs_gradient =
1268 !is_translation &&
1269 (update_flags &
1272
1273 // Use the more heavy VectorizedArray code path if there is more than
1274 // one point left to compute
1275 for (unsigned int i = 0; i < n_points; i += n_lanes)
1276 if (n_points - i > 1)
1277 {
1279 for (unsigned int j = 0; j < n_lanes; ++j)
1280 if (i + j < n_points)
1281 for (unsigned int d = 0; d < dim; ++d)
1282 p_vec[d][j] = unit_points[i + j][d];
1283 else
1284 for (unsigned int d = 0; d < dim; ++d)
1285 p_vec[d][j] = unit_points[i][d];
1286
1289 derivative;
1290 if (needs_gradient)
1291 {
1292 const auto result =
1294 polynomials_1d,
1295 support_points,
1296 p_vec,
1297 polynomials_1d.size() == 2,
1298 renumber_lexicographic_to_hierarchic);
1299
1300 value = result.first;
1301
1302 for (unsigned int d = 0; d < spacedim; ++d)
1303 for (unsigned int e = 0; e < dim; ++e)
1304 derivative[d][e] = result.second[e][d];
1305 }
1306 else
1308 polynomials_1d,
1309 support_points,
1310 p_vec,
1311 polynomials_1d.size() == 2,
1312 renumber_lexicographic_to_hierarchic);
1313
1314 if (update_flags & update_quadrature_points)
1315 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1316 for (unsigned int d = 0; d < spacedim; ++d)
1317 quadrature_points[i + j][d] = value[d][j];
1318
1319 if (is_translation)
1320 continue;
1321
1322 if (update_flags & update_contravariant_transformation)
1323 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1324 for (unsigned int d = 0; d < spacedim; ++d)
1325 for (unsigned int e = 0; e < dim; ++e)
1326 jacobians[i + j][d][e] = derivative[d][e][j];
1327
1328 if (update_flags & update_volume_elements)
1329 {
1330 const auto determinant = derivative.determinant();
1331 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1332 data.volume_elements[i + j] = determinant[j];
1333 }
1334
1335 if (update_flags & update_covariant_transformation)
1336 {
1337 const auto covariant = derivative.covariant_form();
1338 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1339 for (unsigned int d = 0; d < dim; ++d)
1340 for (unsigned int e = 0; e < spacedim; ++e)
1341 inverse_jacobians[i + j][d][e] = covariant[e][d][j];
1342 }
1343 }
1344 else
1345 {
1348 if (needs_gradient)
1349 {
1350 const auto result =
1352 polynomials_1d,
1353 support_points,
1354 unit_points[i],
1355 polynomials_1d.size() == 2,
1356 renumber_lexicographic_to_hierarchic);
1357
1358 value = result.first;
1359
1360 for (unsigned int d = 0; d < spacedim; ++d)
1361 for (unsigned int e = 0; e < dim; ++e)
1362 derivative[d][e] = result.second[e][d];
1363 }
1364 else
1366 polynomials_1d,
1367 support_points,
1368 unit_points[i],
1369 polynomials_1d.size() == 2,
1370 renumber_lexicographic_to_hierarchic);
1371
1372 if (update_flags & update_quadrature_points)
1373 quadrature_points[i] = value;
1374
1375 if (is_translation)
1376 continue;
1377
1378 if (dim == spacedim && update_flags & update_volume_elements)
1379 data.volume_elements[i] = derivative.determinant();
1380
1381 if (update_flags & update_contravariant_transformation)
1382 jacobians[i] = derivative;
1383
1384 if (update_flags & update_covariant_transformation)
1385 inverse_jacobians[i] = jacobians[i].covariant_form().transpose();
1386 }
1387 }
1388
1389
1390
1397 template <int dim, int spacedim>
1398 inline void
1400 const CellSimilarity::Similarity cell_similarity,
1401 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1402 const ArrayView<const Point<dim>> &unit_points,
1403 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1404 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1405 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1406 {
1407 if (data.update_each & update_jacobian_grads)
1408 {
1409 const ArrayView<const Point<spacedim>> support_points(
1410 data.mapping_support_points);
1411 const unsigned int n_q_points = jacobian_grads.size();
1412
1413 if (cell_similarity != CellSimilarity::translation)
1414 for (unsigned int point = 0; point < n_q_points; ++point)
1415 {
1418 polynomials_1d,
1419 support_points,
1420 unit_points[point],
1421 renumber_lexicographic_to_hierarchic);
1422
1423 for (unsigned int i = 0; i < spacedim; ++i)
1424 for (unsigned int j = 0; j < dim; ++j)
1425 for (unsigned int l = 0; l < dim; ++l)
1426 jacobian_grads[point][i][j][l] = second[j][l][i];
1427 }
1428 }
1429 }
1430
1431
1432
1439 template <int dim, int spacedim>
1440 inline void
1442 const CellSimilarity::Similarity cell_similarity,
1443 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1444 const ArrayView<const Point<dim>> &unit_points,
1445 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1446 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1447 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1448 {
1449 if (data.update_each & update_jacobian_pushed_forward_grads)
1450 {
1451 const ArrayView<const Point<spacedim>> support_points(
1452 data.mapping_support_points);
1453 const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
1454
1455 if (cell_similarity != CellSimilarity::translation)
1456 {
1457 double tmp[spacedim][spacedim][spacedim];
1458 for (unsigned int point = 0; point < n_q_points; ++point)
1459 {
1462 polynomials_1d,
1463 support_points,
1464 unit_points[point],
1465 renumber_lexicographic_to_hierarchic);
1466 const DerivativeForm<1, dim, spacedim> covariant =
1467 data.output_data->inverse_jacobians[point].transpose();
1468
1469 // first push forward the j-components
1470 for (unsigned int i = 0; i < spacedim; ++i)
1471 for (unsigned int j = 0; j < spacedim; ++j)
1472 for (unsigned int l = 0; l < dim; ++l)
1473 {
1474 tmp[i][j][l] = second[0][l][i] * covariant[j][0];
1475 for (unsigned int jr = 1; jr < dim; ++jr)
1476 {
1477 tmp[i][j][l] +=
1478 second[jr][l][i] * covariant[j][jr];
1479 }
1480 }
1481
1482 // now, pushing forward the l-components
1483 for (unsigned int i = 0; i < spacedim; ++i)
1484 for (unsigned int j = 0; j < spacedim; ++j)
1485 for (unsigned int l = 0; l < spacedim; ++l)
1486 {
1487 jacobian_pushed_forward_grads[point][i][j][l] =
1488 tmp[i][j][0] * covariant[l][0];
1489 for (unsigned int lr = 1; lr < dim; ++lr)
1490 {
1491 jacobian_pushed_forward_grads[point][i][j][l] +=
1492 tmp[i][j][lr] * covariant[l][lr];
1493 }
1494 }
1495 }
1496 }
1497 }
1498 }
1499
1500
1501
1502 template <int dim, int spacedim, int length_tensor>
1505 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1506 {
1507 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1509 for (unsigned int i = 0; i < spacedim; ++i)
1510 {
1511 if (dim == 1)
1512 result[i][0][0][0] = compressed[0][i];
1513 else if (dim >= 2)
1514 {
1515 for (unsigned int d = 0; d < 2; ++d)
1516 for (unsigned int e = 0; e < 2; ++e)
1517 for (unsigned int f = 0; f < 2; ++f)
1518 result[i][d][e][f] = compressed[d + e + f][i];
1519 if (dim == 3)
1520 {
1521 AssertDimension(length_tensor, 10);
1522
1523 // the derivatives are ordered in rows by increasing z
1524 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1525 // j walks through the columns
1526 for (unsigned int d = 0; d < 2; ++d)
1527 for (unsigned int e = 0; e < 2; ++e)
1528 {
1529 result[i][d][e][2] = compressed[4 + d + e][i];
1530 result[i][d][2][e] = compressed[4 + d + e][i];
1531 result[i][2][d][e] = compressed[4 + d + e][i];
1532 }
1533 for (unsigned int d = 0; d < 2; ++d)
1534 {
1535 result[i][d][2][2] = compressed[7 + d][i];
1536 result[i][2][d][2] = compressed[7 + d][i];
1537 result[i][2][2][d] = compressed[7 + d][i];
1538 }
1539 result[i][2][2][2] = compressed[9][i];
1540 }
1541 }
1542 }
1543 return result;
1544 }
1545
1546
1547
1554 template <int dim, int spacedim>
1555 inline void
1557 const CellSimilarity::Similarity cell_similarity,
1558 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1559 const ArrayView<const Point<dim>> &unit_points,
1560 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1561 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1562 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1563 {
1564 if (data.update_each & update_jacobian_2nd_derivatives)
1565 {
1566 const ArrayView<const Point<spacedim>> support_points(
1567 data.mapping_support_points);
1568 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1569
1570 if (cell_similarity != CellSimilarity::translation)
1571 {
1572 for (unsigned int point = 0; point < n_q_points; ++point)
1573 {
1574 jacobian_2nd_derivatives[point] = expand_3rd_derivatives<dim>(
1576 polynomials_1d,
1577 support_points,
1578 unit_points[point],
1579 renumber_lexicographic_to_hierarchic));
1580 }
1581 }
1582 }
1583 }
1584
1585
1586
1594 template <int dim, int spacedim>
1595 inline void
1597 const CellSimilarity::Similarity cell_similarity,
1598 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1599 const ArrayView<const Point<dim>> &unit_points,
1600 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1601 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1602 std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
1603 {
1605 {
1606 const ArrayView<const Point<spacedim>> support_points(
1607 data.mapping_support_points);
1608 const unsigned int n_q_points =
1609 jacobian_pushed_forward_2nd_derivatives.size();
1610
1611 if (cell_similarity != CellSimilarity::translation)
1612 {
1615 for (unsigned int point = 0; point < n_q_points; ++point)
1616 {
1620 polynomials_1d,
1621 support_points,
1622 unit_points[point],
1623 renumber_lexicographic_to_hierarchic));
1624 const DerivativeForm<1, dim, spacedim> covariant =
1625 data.output_data->inverse_jacobians[point].transpose();
1626
1627 // push forward the j-coordinate
1628 for (unsigned int i = 0; i < spacedim; ++i)
1629 for (unsigned int j = 0; j < spacedim; ++j)
1630 for (unsigned int l = 0; l < dim; ++l)
1631 for (unsigned int m = 0; m < dim; ++m)
1632 {
1633 tmp[i][j][l][m] =
1634 third[i][0][l][m] * covariant[j][0];
1635 for (unsigned int jr = 1; jr < dim; ++jr)
1636 tmp[i][j][l][m] +=
1637 third[i][jr][l][m] * covariant[j][jr];
1638 }
1639
1640 // push forward the l-coordinate
1641 for (unsigned int i = 0; i < spacedim; ++i)
1642 for (unsigned int j = 0; j < spacedim; ++j)
1643 for (unsigned int l = 0; l < spacedim; ++l)
1644 for (unsigned int m = 0; m < dim; ++m)
1645 {
1646 tmp2[i][j][l][m] =
1647 tmp[i][j][0][m] * covariant[l][0];
1648 for (unsigned int lr = 1; lr < dim; ++lr)
1649 tmp2[i][j][l][m] +=
1650 tmp[i][j][lr][m] * covariant[l][lr];
1651 }
1652
1653 // push forward the m-coordinate
1654 for (unsigned int i = 0; i < spacedim; ++i)
1655 for (unsigned int j = 0; j < spacedim; ++j)
1656 for (unsigned int l = 0; l < spacedim; ++l)
1657 for (unsigned int m = 0; m < spacedim; ++m)
1658 {
1659 jacobian_pushed_forward_2nd_derivatives
1660 [point][i][j][l][m] =
1661 tmp2[i][j][l][0] * covariant[m][0];
1662 for (unsigned int mr = 1; mr < dim; ++mr)
1663 jacobian_pushed_forward_2nd_derivatives[point][i]
1664 [j][l]
1665 [m] +=
1666 tmp2[i][j][l][mr] * covariant[m][mr];
1667 }
1668 }
1669 }
1670 }
1671 }
1672
1673
1674
1675 template <int dim, int spacedim, int length_tensor>
1678 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1679 {
1680 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1682 for (unsigned int i = 0; i < spacedim; ++i)
1683 {
1684 if (dim == 1)
1685 result[i][0][0][0][0] = compressed[0][i];
1686 else if (dim >= 2)
1687 {
1688 for (unsigned int d = 0; d < 2; ++d)
1689 for (unsigned int e = 0; e < 2; ++e)
1690 for (unsigned int f = 0; f < 2; ++f)
1691 for (unsigned int g = 0; g < 2; ++g)
1692 result[i][d][e][f][g] = compressed[d + e + f + g][i];
1693 if (dim == 3)
1694 {
1695 AssertDimension(length_tensor, 15);
1696
1697 // the derivatives are ordered in rows by increasing z
1698 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1699 // j walks through the columns
1700 for (unsigned int d = 0; d < 2; ++d)
1701 for (unsigned int e = 0; e < 2; ++e)
1702 for (unsigned int f = 0; f < 2; ++f)
1703 {
1704 result[i][d][e][f][2] = compressed[5 + d + e + f][i];
1705 result[i][d][e][2][f] = compressed[5 + d + e + f][i];
1706 result[i][d][2][e][f] = compressed[5 + d + e + f][i];
1707 result[i][2][d][e][f] = compressed[5 + d + e + f][i];
1708 }
1709 for (unsigned int d = 0; d < 2; ++d)
1710 for (unsigned int e = 0; e < 2; ++e)
1711 {
1712 result[i][d][e][2][2] = compressed[9 + d + e][i];
1713 result[i][d][2][e][2] = compressed[9 + d + e][i];
1714 result[i][d][2][2][e] = compressed[9 + d + e][i];
1715 result[i][2][d][e][2] = compressed[9 + d + e][i];
1716 result[i][2][d][2][e] = compressed[9 + d + e][i];
1717 result[i][2][2][d][e] = compressed[9 + d + e][i];
1718 }
1719 for (unsigned int d = 0; d < 2; ++d)
1720 {
1721 result[i][d][2][2][2] = compressed[12 + d][i];
1722 result[i][2][d][2][2] = compressed[12 + d][i];
1723 result[i][2][2][d][2] = compressed[12 + d][i];
1724 result[i][2][2][2][d] = compressed[12 + d][i];
1725 }
1726 result[i][2][2][2][2] = compressed[14][i];
1727 }
1728 }
1729 }
1730 return result;
1731 }
1732
1733
1734
1741 template <int dim, int spacedim>
1742 inline void
1744 const CellSimilarity::Similarity cell_similarity,
1745 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1746 const ArrayView<const Point<dim>> &unit_points,
1747 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1748 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1749 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1750 {
1751 if (data.update_each & update_jacobian_3rd_derivatives)
1752 {
1753 const ArrayView<const Point<spacedim>> support_points(
1754 data.mapping_support_points);
1755 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1756
1757 if (cell_similarity != CellSimilarity::translation)
1758 {
1759 for (unsigned int point = 0; point < n_q_points; ++point)
1760 {
1761 jacobian_3rd_derivatives[point] = expand_4th_derivatives<dim>(
1763 polynomials_1d,
1764 support_points,
1765 unit_points[point],
1766 renumber_lexicographic_to_hierarchic));
1767 }
1768 }
1769 }
1770 }
1771
1772
1773
1781 template <int dim, int spacedim>
1782 inline void
1784 const CellSimilarity::Similarity cell_similarity,
1785 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1786 const ArrayView<const Point<dim>> &unit_points,
1787 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1788 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1789 std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
1790 {
1792 {
1793 const ArrayView<const Point<spacedim>> support_points(
1794 data.mapping_support_points);
1795 const unsigned int n_q_points =
1796 jacobian_pushed_forward_3rd_derivatives.size();
1797
1798 if (cell_similarity != CellSimilarity::translation)
1799 {
1800 ::
1801 ndarray<double, spacedim, spacedim, spacedim, spacedim, dim>
1802 tmp;
1804 tmp2;
1805 for (unsigned int point = 0; point < n_q_points; ++point)
1806 {
1810 polynomials_1d,
1811 support_points,
1812 unit_points[point],
1813 renumber_lexicographic_to_hierarchic));
1814
1815 const DerivativeForm<1, dim, spacedim> covariant =
1816 data.output_data->inverse_jacobians[point].transpose();
1817
1818 // push-forward the j-coordinate
1819 for (unsigned int i = 0; i < spacedim; ++i)
1820 for (unsigned int j = 0; j < spacedim; ++j)
1821 for (unsigned int l = 0; l < dim; ++l)
1822 for (unsigned int m = 0; m < dim; ++m)
1823 for (unsigned int n = 0; n < dim; ++n)
1824 {
1825 tmp[i][j][l][m][n] =
1826 fourth[i][0][l][m][n] * covariant[j][0];
1827 for (unsigned int jr = 1; jr < dim; ++jr)
1828 tmp[i][j][l][m][n] +=
1829 fourth[i][jr][l][m][n] * covariant[j][jr];
1830 }
1831
1832 // push-forward the l-coordinate
1833 for (unsigned int i = 0; i < spacedim; ++i)
1834 for (unsigned int j = 0; j < spacedim; ++j)
1835 for (unsigned int l = 0; l < spacedim; ++l)
1836 for (unsigned int m = 0; m < dim; ++m)
1837 for (unsigned int n = 0; n < dim; ++n)
1838 {
1839 tmp2[i][j][l][m][n] =
1840 tmp[i][j][0][m][n] * covariant[l][0];
1841 for (unsigned int lr = 1; lr < dim; ++lr)
1842 tmp2[i][j][l][m][n] +=
1843 tmp[i][j][lr][m][n] * covariant[l][lr];
1844 }
1845
1846 // push-forward the m-coordinate
1847 for (unsigned int i = 0; i < spacedim; ++i)
1848 for (unsigned int j = 0; j < spacedim; ++j)
1849 for (unsigned int l = 0; l < spacedim; ++l)
1850 for (unsigned int m = 0; m < spacedim; ++m)
1851 for (unsigned int n = 0; n < dim; ++n)
1852 {
1853 tmp[i][j][l][m][n] =
1854 tmp2[i][j][l][0][n] * covariant[m][0];
1855 for (unsigned int mr = 1; mr < dim; ++mr)
1856 tmp[i][j][l][m][n] +=
1857 tmp2[i][j][l][mr][n] * covariant[m][mr];
1858 }
1859
1860 // push-forward the n-coordinate
1861 for (unsigned int i = 0; i < spacedim; ++i)
1862 for (unsigned int j = 0; j < spacedim; ++j)
1863 for (unsigned int l = 0; l < spacedim; ++l)
1864 for (unsigned int m = 0; m < spacedim; ++m)
1865 for (unsigned int n = 0; n < spacedim; ++n)
1866 {
1867 jacobian_pushed_forward_3rd_derivatives
1868 [point][i][j][l][m][n] =
1869 tmp[i][j][l][m][0] * covariant[n][0];
1870 for (unsigned int nr = 1; nr < dim; ++nr)
1871 jacobian_pushed_forward_3rd_derivatives[point]
1872 [i][j][l]
1873 [m][n] +=
1874 tmp[i][j][l][m][nr] * covariant[n][nr];
1875 }
1876 }
1877 }
1878 }
1879 }
1880
1881
1882
1892 template <int dim, int spacedim>
1893 inline void
1895 const ::MappingQ<dim, spacedim> &mapping,
1896 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
1897 const unsigned int face_no,
1898 const unsigned int subface_no,
1899 const unsigned int n_q_points,
1900 const std::vector<double> &weights,
1901 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1903 &output_data)
1904 {
1905 const UpdateFlags update_flags = data.update_each;
1906
1907 if (update_flags &
1909 {
1910 if (update_flags & update_boundary_forms)
1911 AssertDimension(output_data.boundary_forms.size(), n_q_points);
1912 if (update_flags & update_normal_vectors)
1913 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1914 if (update_flags & update_JxW_values)
1915 AssertDimension(output_data.JxW_values.size(), n_q_points);
1916
1917 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1918
1919 // first compute some common data that is used for evaluating
1920 // all of the flags below
1921
1922 // map the unit tangentials to the real cell. checking for d!=dim-1
1923 // eliminates compiler warnings regarding unsigned int expressions <
1924 // 0.
1925 for (unsigned int d = 0; d != dim - 1; ++d)
1926 {
1927 const unsigned int vector_index =
1929 Assert(vector_index < data.unit_tangentials.size(),
1931 Assert(data.aux[d].size() <=
1932 data.unit_tangentials[vector_index].size(),
1934 mapping.transform(make_array_view(
1935 data.unit_tangentials[vector_index]),
1937 data,
1938 make_array_view(data.aux[d]));
1939 }
1940
1941 if (update_flags & update_boundary_forms)
1942 {
1943 // if dim==spacedim, we can use the unit tangentials to compute
1944 // the boundary form by simply taking the cross product
1945 if (dim == spacedim)
1946 {
1947 for (unsigned int i = 0; i < n_q_points; ++i)
1948 switch (dim)
1949 {
1950 case 1:
1951 // in 1d, we don't have access to any of the
1952 // data.aux fields (because it has only dim-1
1953 // components), but we can still compute the
1954 // boundary form by simply looking at the number of
1955 // the face
1956 output_data.boundary_forms[i][0] =
1957 (face_no == 0 ? -1 : +1);
1958 break;
1959 case 2:
1960 output_data.boundary_forms[i] =
1961 cross_product_2d(data.aux[0][i]);
1962 break;
1963 case 3:
1964 output_data.boundary_forms[i] =
1965 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1966 break;
1967 default:
1969 }
1970 }
1971 else //(dim < spacedim)
1972 {
1973 // in the codim-one case, the boundary form results from the
1974 // cross product of all the face tangential vectors and the
1975 // cell normal vector
1976 //
1977 // to compute the cell normal, use the same method used in
1978 // fill_fe_values for cells above
1979 AssertDimension(data.output_data->jacobians.size(),
1980 n_q_points);
1981
1982 for (unsigned int point = 0; point < n_q_points; ++point)
1983 {
1984 const DerivativeForm<1, dim, spacedim> contravariant =
1985 data.output_data->jacobians[point];
1986 if (dim == 1)
1987 {
1988 // J is a tangent vector
1989 output_data.boundary_forms[point] =
1990 contravariant.transpose()[0];
1991 output_data.boundary_forms[point] /=
1992 (face_no == 0 ? -1. : +1.) *
1993 output_data.boundary_forms[point].norm();
1994 }
1995
1996 if (dim == 2)
1997 {
1999 contravariant.transpose();
2000
2001 Tensor<1, spacedim> cell_normal =
2002 cross_product_3d(DX_t[0], DX_t[1]);
2003 cell_normal /= cell_normal.norm();
2004
2005 // then compute the face normal from the face
2006 // tangent and the cell normal:
2007 output_data.boundary_forms[point] =
2008 cross_product_3d(data.aux[0][point], cell_normal);
2009 }
2010 }
2011 }
2012 }
2013
2014 if (update_flags & update_JxW_values)
2015 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
2016 {
2017 output_data.JxW_values[i] =
2018 output_data.boundary_forms[i].norm() * weights[i];
2019
2020 if (subface_no != numbers::invalid_unsigned_int)
2021 {
2022 const double area_ratio = GeometryInfo<dim>::subface_ratio(
2023 cell->subface_case(face_no), subface_no);
2024 output_data.JxW_values[i] *= area_ratio;
2025 }
2026 }
2027
2028 if (update_flags & update_normal_vectors)
2029 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
2030 output_data.normal_vectors[i] =
2031 Point<spacedim>(output_data.boundary_forms[i] /
2032 output_data.boundary_forms[i].norm());
2033 }
2034 }
2035
2036
2037
2044 template <int dim, int spacedim>
2045 inline void
2047 const ::MappingQ<dim, spacedim> &mapping,
2048 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
2049 const unsigned int face_no,
2050 const unsigned int subface_no,
2051 const typename QProjector<dim>::DataSetDescriptor data_set,
2052 const Quadrature<dim - 1> &quadrature,
2053 const typename ::MappingQ<dim, spacedim>::InternalData &data,
2054 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
2055 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
2057 &output_data)
2058 {
2059 const ArrayView<const Point<dim>> quadrature_points(
2060 &data.quadrature_points[data_set], quadrature.size());
2061 if (dim > 1 && data.tensor_product_quadrature)
2062 {
2065 data,
2066 output_data.quadrature_points,
2067 output_data.jacobians,
2068 output_data.inverse_jacobians,
2069 output_data.jacobian_grads);
2070 }
2071 else
2072 {
2076 data,
2077 quadrature_points,
2078 polynomials_1d,
2079 renumber_lexicographic_to_hierarchic,
2080 output_data.quadrature_points,
2081 output_data.jacobians,
2082 output_data.inverse_jacobians);
2085 data,
2086 quadrature_points,
2087 polynomials_1d,
2088 renumber_lexicographic_to_hierarchic,
2089 output_data.jacobian_grads);
2090 }
2093 data,
2094 quadrature_points,
2095 polynomials_1d,
2096 renumber_lexicographic_to_hierarchic,
2097 output_data.jacobian_pushed_forward_grads);
2100 data,
2101 quadrature_points,
2102 polynomials_1d,
2103 renumber_lexicographic_to_hierarchic,
2104 output_data.jacobian_2nd_derivatives);
2107 data,
2108 quadrature_points,
2109 polynomials_1d,
2110 renumber_lexicographic_to_hierarchic,
2114 data,
2115 quadrature_points,
2116 polynomials_1d,
2117 renumber_lexicographic_to_hierarchic,
2118 output_data.jacobian_3rd_derivatives);
2121 data,
2122 quadrature_points,
2123 polynomials_1d,
2124 renumber_lexicographic_to_hierarchic,
2126
2128 cell,
2129 face_no,
2130 subface_no,
2131 quadrature.size(),
2132 quadrature.get_weights(),
2133 data,
2134 output_data);
2135 }
2136
2137
2138
2142 template <int dim, int spacedim, int rank>
2143 inline void
2145 const ArrayView<const Tensor<rank, dim>> &input,
2146 const MappingKind mapping_kind,
2147 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2148 const ArrayView<Tensor<rank, spacedim>> &output)
2149 {
2150 AssertDimension(input.size(), output.size());
2151 Assert((dynamic_cast<
2152 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2153 &mapping_data) != nullptr),
2155 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2156 static_cast<
2157 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2158 mapping_data);
2159
2160 switch (mapping_kind)
2161 {
2163 {
2166 "update_contravariant_transformation"));
2167
2168 for (unsigned int i = 0; i < output.size(); ++i)
2169 output[i] = apply_transformation(data.output_data->jacobians[i],
2170 input[i]);
2171
2172 return;
2173 }
2174
2175 case mapping_piola:
2176 {
2179 "update_contravariant_transformation"));
2180 Assert(data.update_each & update_volume_elements,
2182 "update_volume_elements"));
2183 Assert(rank == 1, ExcMessage("Only for rank 1"));
2184 if (rank != 1)
2185 return;
2186
2187 for (unsigned int i = 0; i < output.size(); ++i)
2188 {
2189 output[i] =
2190 apply_transformation(data.output_data->jacobians[i],
2191 input[i]);
2192 output[i] /= data.volume_elements[i];
2193 }
2194 return;
2195 }
2196 // We still allow this operation as in the
2197 // reference cell Derivatives are Tensor
2198 // rather than DerivativeForm
2199 case mapping_covariant:
2200 {
2203 "update_covariant_transformation"));
2204
2205 for (unsigned int i = 0; i < output.size(); ++i)
2206 output[i] = apply_transformation(
2207 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2208
2209 return;
2210 }
2211
2212 default:
2214 }
2215 }
2216
2217
2218
2222 template <int dim, int spacedim, int rank>
2223 inline void
2225 const ArrayView<const Tensor<rank, dim>> &input,
2226 const MappingKind mapping_kind,
2227 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2228 const ArrayView<Tensor<rank, spacedim>> &output)
2229 {
2230 AssertDimension(input.size(), output.size());
2231 Assert((dynamic_cast<
2232 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2233 &mapping_data) != nullptr),
2235 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2236 static_cast<
2237 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2238 mapping_data);
2239
2240 switch (mapping_kind)
2241 {
2243 {
2244 Assert(data.update_each & update_covariant_transformation,
2246 "update_covariant_transformation"));
2249 "update_contravariant_transformation"));
2250 Assert(rank == 2, ExcMessage("Only for rank 2"));
2251
2252 for (unsigned int i = 0; i < output.size(); ++i)
2253 {
2255 apply_transformation(data.output_data->jacobians[i],
2256 transpose(input[i]));
2257 output[i] = apply_transformation(
2258 data.output_data->inverse_jacobians[i].transpose(),
2259 A.transpose());
2260 }
2261
2262 return;
2263 }
2264
2266 {
2267 Assert(data.update_each & update_covariant_transformation,
2269 "update_covariant_transformation"));
2270 Assert(rank == 2, ExcMessage("Only for rank 2"));
2271
2272 for (unsigned int i = 0; i < output.size(); ++i)
2273 {
2274 const DerivativeForm<1, dim, spacedim> covariant =
2275 data.output_data->inverse_jacobians[i].transpose();
2277 apply_transformation(covariant, transpose(input[i]));
2278 output[i] = apply_transformation(covariant, A.transpose());
2279 }
2280
2281 return;
2282 }
2283
2285 {
2286 Assert(data.update_each & update_covariant_transformation,
2288 "update_covariant_transformation"));
2291 "update_contravariant_transformation"));
2292 Assert(data.update_each & update_volume_elements,
2294 "update_volume_elements"));
2295 Assert(rank == 2, ExcMessage("Only for rank 2"));
2296
2297 for (unsigned int i = 0; i < output.size(); ++i)
2298 {
2299 const DerivativeForm<1, dim, spacedim> covariant =
2300 data.output_data->inverse_jacobians[i].transpose();
2302 apply_transformation(covariant, input[i]);
2303 const Tensor<2, spacedim> T =
2304 apply_transformation(data.output_data->jacobians[i],
2305 A.transpose());
2306
2307 output[i] = transpose(T);
2308 output[i] /= data.volume_elements[i];
2309 }
2310
2311 return;
2312 }
2313
2314 default:
2316 }
2317 }
2318
2319
2320
2324 template <int dim, int spacedim>
2325 inline void
2327 const ArrayView<const Tensor<3, dim>> &input,
2328 const MappingKind mapping_kind,
2329 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2330 const ArrayView<Tensor<3, spacedim>> &output)
2331 {
2332 AssertDimension(input.size(), output.size());
2333 Assert((dynamic_cast<
2334 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2335 &mapping_data) != nullptr),
2337 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2338 static_cast<
2339 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2340 mapping_data);
2341
2342 switch (mapping_kind)
2343 {
2345 {
2346 Assert(data.update_each & update_covariant_transformation,
2348 "update_covariant_transformation"));
2351 "update_contravariant_transformation"));
2352
2353 for (unsigned int q = 0; q < output.size(); ++q)
2354 {
2355 const DerivativeForm<1, dim, spacedim> covariant =
2356 data.output_data->inverse_jacobians[q].transpose();
2357 const DerivativeForm<1, dim, spacedim> contravariant =
2358 data.output_data->jacobians[q];
2359
2360 for (unsigned int i = 0; i < spacedim; ++i)
2361 {
2362 double tmp1[dim][dim];
2363 for (unsigned int J = 0; J < dim; ++J)
2364 for (unsigned int K = 0; K < dim; ++K)
2365 {
2366 tmp1[J][K] =
2367 contravariant[i][0] * input[q][0][J][K];
2368 for (unsigned int I = 1; I < dim; ++I)
2369 tmp1[J][K] +=
2370 contravariant[i][I] * input[q][I][J][K];
2371 }
2372 for (unsigned int j = 0; j < spacedim; ++j)
2373 {
2374 double tmp2[dim];
2375 for (unsigned int K = 0; K < dim; ++K)
2376 {
2377 tmp2[K] = covariant[j][0] * tmp1[0][K];
2378 for (unsigned int J = 1; J < dim; ++J)
2379 tmp2[K] += covariant[j][J] * tmp1[J][K];
2380 }
2381 for (unsigned int k = 0; k < spacedim; ++k)
2382 {
2383 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2384 for (unsigned int K = 1; K < dim; ++K)
2385 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2386 }
2387 }
2388 }
2389 }
2390 return;
2391 }
2392
2394 {
2395 Assert(data.update_each & update_covariant_transformation,
2397 "update_covariant_transformation"));
2398
2399 for (unsigned int q = 0; q < output.size(); ++q)
2400 {
2401 const DerivativeForm<1, dim, spacedim> covariant =
2402 data.output_data->inverse_jacobians[q].transpose();
2403
2404 for (unsigned int i = 0; i < spacedim; ++i)
2405 {
2406 double tmp1[dim][dim];
2407 for (unsigned int J = 0; J < dim; ++J)
2408 for (unsigned int K = 0; K < dim; ++K)
2409 {
2410 tmp1[J][K] = covariant[i][0] * input[q][0][J][K];
2411 for (unsigned int I = 1; I < dim; ++I)
2412 tmp1[J][K] += covariant[i][I] * input[q][I][J][K];
2413 }
2414 for (unsigned int j = 0; j < spacedim; ++j)
2415 {
2416 double tmp2[dim];
2417 for (unsigned int K = 0; K < dim; ++K)
2418 {
2419 tmp2[K] = covariant[j][0] * tmp1[0][K];
2420 for (unsigned int J = 1; J < dim; ++J)
2421 tmp2[K] += covariant[j][J] * tmp1[J][K];
2422 }
2423 for (unsigned int k = 0; k < spacedim; ++k)
2424 {
2425 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2426 for (unsigned int K = 1; K < dim; ++K)
2427 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2428 }
2429 }
2430 }
2431 }
2432
2433 return;
2434 }
2435
2437 {
2438 Assert(data.update_each & update_covariant_transformation,
2440 "update_covariant_transformation"));
2443 "update_contravariant_transformation"));
2444 Assert(data.update_each & update_volume_elements,
2446 "update_volume_elements"));
2447
2448 for (unsigned int q = 0; q < output.size(); ++q)
2449 {
2450 const DerivativeForm<1, dim, spacedim> covariant =
2451 data.output_data->inverse_jacobians[q].transpose();
2452 const DerivativeForm<1, dim, spacedim> contravariant =
2453 data.output_data->jacobians[q];
2454 for (unsigned int i = 0; i < spacedim; ++i)
2455 {
2456 double factor[dim];
2457 for (unsigned int I = 0; I < dim; ++I)
2458 factor[I] =
2459 contravariant[i][I] * (1. / data.volume_elements[q]);
2460 double tmp1[dim][dim];
2461 for (unsigned int J = 0; J < dim; ++J)
2462 for (unsigned int K = 0; K < dim; ++K)
2463 {
2464 tmp1[J][K] = factor[0] * input[q][0][J][K];
2465 for (unsigned int I = 1; I < dim; ++I)
2466 tmp1[J][K] += factor[I] * input[q][I][J][K];
2467 }
2468 for (unsigned int j = 0; j < spacedim; ++j)
2469 {
2470 double tmp2[dim];
2471 for (unsigned int K = 0; K < dim; ++K)
2472 {
2473 tmp2[K] = covariant[j][0] * tmp1[0][K];
2474 for (unsigned int J = 1; J < dim; ++J)
2475 tmp2[K] += covariant[j][J] * tmp1[J][K];
2476 }
2477 for (unsigned int k = 0; k < spacedim; ++k)
2478 {
2479 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2480 for (unsigned int K = 1; K < dim; ++K)
2481 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2482 }
2483 }
2484 }
2485 }
2486
2487 return;
2488 }
2489
2490 default:
2492 }
2493 }
2494
2495
2496
2501 template <int dim, int spacedim, int rank>
2502 inline void
2505 const MappingKind mapping_kind,
2506 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2508 {
2509 AssertDimension(input.size(), output.size());
2510 Assert((dynamic_cast<
2511 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2512 &mapping_data) != nullptr),
2514 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2515 static_cast<
2516 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2517 mapping_data);
2518
2519 switch (mapping_kind)
2520 {
2521 case mapping_covariant:
2522 {
2525 "update_covariant_transformation"));
2526
2527 for (unsigned int i = 0; i < output.size(); ++i)
2528 output[i] = apply_transformation(
2529 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2530
2531 return;
2532 }
2533 default:
2535 }
2536 }
2537 } // namespace MappingQImplementation
2538} // namespace internal
2539
2541
2542#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
std::size_t size() const
Definition array_view.h:684
Number determinant() const
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
DerivativeForm< 1, spacedim, dim, Number > transpose() const
void set_data_pointers(AlignedVector< Number > *scratch_data, const unsigned int n_components)
const Number * begin_gradients() const
const Number * begin_values() const
const Number * begin_dof_values() const
const Number * begin_hessians() const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
const std::vector< double > & get_weights() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
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
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)
InverseQuadraticApproximation(const InverseQuadraticApproximation &)=default
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Point< 2 > second
Definition grid_out.cc:4624
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
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:79
@ mapping_piola
Definition mapping.h:114
@ mapping_covariant_gradient
Definition mapping.h:100
@ mapping_covariant
Definition mapping.h:89
@ mapping_contravariant
Definition mapping.h:94
@ mapping_contravariant_hessian
Definition mapping.h:156
@ mapping_covariant_hessian
Definition mapping.h:150
@ mapping_contravariant_gradient
Definition mapping.h:106
@ mapping_piola_gradient
Definition mapping.h:120
@ mapping_piola_hessian
Definition mapping.h:162
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
LogStream deallog
Definition logstream.cc:36
EvaluationFlags
The EvaluationFlags enum.
std::vector< unsigned int > hierarchic_to_lexicographic_numbering(unsigned int degree)
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim > > &vertices)
constexpr T fixed_power(const T t)
Definition utilities.h:942
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
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={})
Tensor< 1, 1, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_higher_derivatives(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< 1, Number2 > &p, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition types.h:220
::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 > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
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 > &)
Number compare_and_apply_mask(const Number &left, const Number &right, const Number &true_value, const Number &false_value)