Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mapping_q_internal.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_mapping_q_internal_h
17#define dealii_mapping_q_internal_h
18
19#include <deal.II/base/config.h>
20
24#include <deal.II/base/point.h>
26#include <deal.II/base/table.h>
27#include <deal.II/base/tensor.h>
28
29#include <deal.II/fe/fe_tools.h>
33
35
42
43#include <array>
44#include <limits>
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)),
151 std::max(std::abs(x2), std::abs(x3)));
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 =
162 std::max(std::abs(y2), std::abs(y3)));
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.
181 Assert(false, ExcInternalError());
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
196 Assert(false, ExcInternalError());
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 {
250 ::Table<2, double> loqvs;
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 {
302 ::Table<2, double> lohvs;
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) =
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.
463 Point<spacedim> p_real;
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>
479 const Point<spacedim, Number> & p,
480 const Point<dim, Number> & initial_p_unit,
481 const std::vector<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 {
486 if (print_iterations_to_deallog)
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
503 Point<dim, Number> p_unit = initial_p_unit;
505 polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber);
506
507 Tensor<1, spacedim, Number> f = p_real.first - p;
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
555 Point<dim, Number> invalid_point;
556 invalid_point[0] = std::numeric_limits<double>::infinity();
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 {
565 if (print_iterations_to_deallog)
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 determinand(df) > 0 on all SIMD lanes
576 if (!(std::min(determinant(df),
577 Number(std::numeric_limits<double>::min())) ==
578 Number(std::numeric_limits<double>::min())))
579 {
580 // We allow to enter this function with an initial guess
581 // outside the unit cell. We might have run into invalid
582 // Jacobians due to that, so we should at least try once to go
583 // back to the unit cell and go on with the Newton iteration
584 // from there. Since the outside case is unlikely, we can
585 // afford spending the extra effort at this place.
586 if (tried_project_to_unit_cell == false)
587 {
590 polynomials_1d,
591 points,
592 p_unit,
593 polynomials_1d.size() == 2,
594 renumber);
595 f = p_real.first - p;
596 f_weighted_norm_square = 1.;
597 last_f_weighted_norm_square = 1;
598 tried_project_to_unit_cell = true;
599 continue;
600 }
601 else
602 return invalid_point;
603 }
604
605 // Solve [f'(x)]d=f(x)
606 const Tensor<2, spacedim, Number> df_inverse = invert(df);
607 const Tensor<1, spacedim, Number> delta = df_inverse * f;
608 last_f_weighted_norm_square = delta.norm_square();
609
610 if (print_iterations_to_deallog)
611 deallog << " delta=" << delta << std::endl;
612
613 // do a line search
614 double step_length = 1.0;
615 do
616 {
617 // update of p_unit. The spacedim-th component of transformed
618 // point is simply ignored in codimension one case. When this
619 // component is not zero, then we are projecting the point to
620 // the surface or curve identified by the cell.
621 Point<dim, Number> p_unit_trial = p_unit;
622 for (unsigned int i = 0; i < dim; ++i)
623 p_unit_trial[i] -= step_length * delta[i];
624
625 // shape values and derivatives at new p_unit point
626 const auto p_real_trial =
628 polynomials_1d,
629 points,
630 p_unit_trial,
631 polynomials_1d.size() == 2,
632 renumber);
633 const Tensor<1, spacedim, Number> f_trial =
634 p_real_trial.first - p;
635 f_weighted_norm_square = (df_inverse * f_trial).norm_square();
636
637 if (print_iterations_to_deallog)
638 {
639 deallog << " step_length=" << step_length << std::endl;
640 if (step_length == 1.0)
641 deallog << " ||f || =" << f.norm() << std::endl;
642 deallog << " ||f*|| =" << f_trial.norm() << std::endl
643 << " ||f*||_A ="
644 << std::sqrt(f_weighted_norm_square) << std::endl;
645 }
646
647 // See if we are making progress with the current step length
648 // and if not, reduce it by a factor of two and try again.
649 //
650 // Strictly speaking, we should probably use the same norm as we
651 // use for the outer algorithm. In practice, line search is just
652 // a crutch to find a "reasonable" step length, and so using the
653 // l2 norm is probably just fine.
654 //
655 // check f_trial.norm() < f.norm() in SIMD form. This is a bit
656 // more complicated because some SIMD lanes might not be doing
657 // any progress any more as they have already reached roundoff
658 // accuracy. We define that as the case of updates less than
659 // 1e-6. The tolerance might seem coarse but since we are
660 // dealing with a Newton iteration of a polynomial function we
661 // either converge quadratically or not any more. Thus, our
662 // selection is to terminate if either the norm of f is
663 // decreasing or that threshold of 1e-6 is reached.
664 if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) *
665 std::max(f_trial.norm_square() - f.norm_square(),
666 Number(0.)) ==
667 Number(0.))
668 {
669 p_real = p_real_trial;
670 p_unit = p_unit_trial;
671 f = f_trial;
672 break;
673 }
674 else if (step_length > 0.05)
675 step_length *= 0.5;
676 else
677 break;
678 }
679 while (true);
680
681 // In case we terminated the line search due to the step becoming
682 // too small, we give the iteration another try with the
683 // projection of the initial guess to the unit cell before we give
684 // up, just like for the negative determinant case.
685 if (step_length <= 0.05 && tried_project_to_unit_cell == false)
686 {
689 polynomials_1d,
690 points,
691 p_unit,
692 polynomials_1d.size() == 2,
693 renumber);
694 f = p_real.first - p;
695 f_weighted_norm_square = 1.;
696 last_f_weighted_norm_square = 1;
697 tried_project_to_unit_cell = true;
698 continue;
699 }
700 else if (step_length <= 0.05)
701 return invalid_point;
702
703 ++newton_iteration;
704 if (newton_iteration > newton_iteration_limit)
705 return invalid_point;
706 }
707 // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the
708 // weighted norm is less than 1e-6 and the improvement against the
709 // previous step was less than a factor of 10 (in that regime, we
710 // either have quadratic convergence or roundoff errors due to a bad
711 // mapping)
712 while (
713 !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) *
714 std::max(last_f_weighted_norm_square -
715 std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) *
716 100.,
717 Number(0.)) ==
718 Number(0.)));
719
720 if (print_iterations_to_deallog)
721 deallog << "Iteration converged for p_unit = [ " << p_unit
722 << " ] and iteration error "
723 << std::sqrt(f_weighted_norm_square) << std::endl;
724
725 return p_unit;
726 }
727
728
729
733 template <int dim>
734 inline Point<dim>
736 const Point<dim + 1> & p,
737 const Point<dim> & initial_p_unit,
738 const std::vector<Point<dim + 1>> & points,
739 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
740 const std::vector<unsigned int> & renumber)
741 {
742 const unsigned int spacedim = dim + 1;
743
744 AssertDimension(points.size(),
745 Utilities::pow(polynomials_1d.size(), dim));
746
747 Point<dim> p_unit = initial_p_unit;
748
749 const double eps = 1.e-12;
750 const unsigned int loop_limit = 10;
751
752 unsigned int loop = 0;
753 double f_weighted_norm_square = 1.;
754
755 while (f_weighted_norm_square > eps * eps && loop++ < loop_limit)
756 {
757 const auto p_real =
759 polynomials_1d,
760 points,
761 p_unit,
762 polynomials_1d.size() == 2,
763 renumber);
764 Tensor<1, spacedim> p_minus_F = p - p_real.first;
765 const DerivativeForm<1, spacedim, dim> DF = p_real.second;
766
768 polynomials_1d, points, p_unit, renumber);
769 Point<dim> f;
771 for (unsigned int j = 0; j < dim; ++j)
772 {
773 f[j] = DF[j] * p_minus_F;
774 for (unsigned int l = 0; l < dim; ++l)
775 df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F;
776 }
777
778 // Solve [df(x)]d=f(x)
779 const Tensor<1, dim> d =
780 invert(df) * static_cast<const Tensor<1, dim> &>(f);
781 f_weighted_norm_square = d.norm_square();
782 p_unit -= d;
783 }
784
785
786 // Here we check that in the last execution of while the first
787 // condition was already wrong, meaning the residual was below
788 // eps. Only if the first condition failed, loop will have been
789 // increased and tested, and thus have reached the limit.
790 AssertThrow(loop < loop_limit,
792
793 return p_unit;
794 }
795
796
797
815 template <int dim, int spacedim>
817 {
818 public:
822 static constexpr unsigned int n_functions =
823 (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
824
837 const std::vector<Point<spacedim>> &real_support_points,
838 const std::vector<Point<dim>> & unit_support_points)
839 : normalization_shift(real_support_points[0])
841 1. / real_support_points[0].distance(real_support_points[1]))
842 , is_affine(true)
843 {
844 AssertDimension(real_support_points.size(), unit_support_points.size());
845
846 // For the bi-/trilinear approximation, we cannot build a quadratic
847 // polynomial due to a lack of points (interpolation matrix would
848 // get singular). Similarly, it is not entirely clear how to gather
849 // enough information for the case dim < spacedim.
850 //
851 // In both cases we require the vector real_support_points to
852 // contain the vertex positions and fall back to an affine
853 // approximation:
854 Assert(dim == spacedim || real_support_points.size() ==
857 if (real_support_points.size() == GeometryInfo<dim>::vertices_per_cell)
858 {
859 const auto affine = GridTools::affine_cell_approximation<dim>(
860 make_array_view(real_support_points));
862 affine.first.covariant_form().transpose();
863
864 // The code for evaluation assumes an additional transformation of
865 // the form (x - normalization_shift) * normalization_length --
866 // account for this in the definition of the coefficients.
867 coefficients[0] =
868 apply_transformation(A_inv, normalization_shift - affine.second);
869 for (unsigned int d = 0; d < spacedim; ++d)
870 for (unsigned int e = 0; e < dim; ++e)
871 coefficients[1 + d][e] =
872 A_inv[e][d] * (1.0 / normalization_length);
873 is_affine = true;
874 return;
875 }
876
878 std::array<double, n_functions> shape_values;
879 for (unsigned int q = 0; q < unit_support_points.size(); ++q)
880 {
881 // Evaluate quadratic shape functions in point, with the
882 // normalization applied in order to avoid roundoff issues with
883 // scaling far away from 1.
884 shape_values[0] = 1.;
885 const Tensor<1, spacedim> p_scaled =
887 (real_support_points[q] - normalization_shift);
888 for (unsigned int d = 0; d < spacedim; ++d)
889 shape_values[1 + d] = p_scaled[d];
890 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
891 for (unsigned int e = 0; e <= d; ++e, ++c)
892 shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
893
894 // Build lower diagonal of least squares matrix and rhs, the
895 // essential part being that we construct the matrix with the
896 // real points and the right hand side by comparing to the
897 // reference point positions which sets up an inverse
898 // interpolation.
899 for (unsigned int i = 0; i < n_functions; ++i)
900 for (unsigned int j = 0; j < n_functions; ++j)
901 matrix[i][j] += shape_values[i] * shape_values[j];
902 for (unsigned int i = 0; i < n_functions; ++i)
903 coefficients[i] += shape_values[i] * unit_support_points[q];
904 }
905
906 // Factorize the matrix A = L * L^T in-place with the
907 // Cholesky-Banachiewicz algorithm. The implementation is similar to
908 // FullMatrix::cholesky() but re-implemented to avoid memory
909 // allocations and some unnecessary divisions which we can do here as
910 // we only need to solve with dim right hand sides.
911 for (unsigned int i = 0; i < n_functions; ++i)
912 {
913 double Lij_sum = 0;
914 for (unsigned int j = 0; j < i; ++j)
915 {
916 double Lik_Ljk_sum = 0;
917 for (unsigned int k = 0; k < j; ++k)
918 Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
919 matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
920 Lij_sum += matrix[i][j] * matrix[i][j];
921 }
922 AssertThrow(matrix[i][i] - Lij_sum >= 0,
923 ExcMessage("Matrix of normal equations not positive "
924 "definite"));
925
926 // Store the inverse in the diagonal since that is the quantity
927 // needed later in the factorization as well as the forward and
928 // backward substitution, minimizing the number of divisions.
929 matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
930 }
931
932 // Solve lower triangular part, L * y = rhs.
933 for (unsigned int i = 0; i < n_functions; ++i)
934 {
935 Point<dim> sum = coefficients[i];
936 for (unsigned int j = 0; j < i; ++j)
937 sum -= matrix[i][j] * coefficients[j];
938 coefficients[i] = sum * matrix[i][i];
939 }
940
941 // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
942 for (unsigned int i = n_functions; i > 0;)
943 {
944 --i;
945 Point<dim> sum = coefficients[i];
946 for (unsigned int j = i + 1; j < n_functions; ++j)
947 sum -= matrix[j][i] * coefficients[j];
948 coefficients[i] = sum * matrix[i][i];
949 }
950
951 // Check whether the approximation is indeed affine, allowing to
952 // skip the quadratic terms.
953 is_affine = true;
954 for (unsigned int i = dim + 1; i < n_functions; ++i)
955 if (coefficients[i].norm_square() > 1e-20)
956 {
957 is_affine = false;
958 break;
959 }
960 }
961
966 default;
967
971 template <typename Number>
974 {
975 Point<dim, Number> result;
976 for (unsigned int d = 0; d < dim; ++d)
977 result[d] = coefficients[0][d];
978
979 // Apply the normalization to ensure a good conditioning. Since Number
980 // might be a vectorized array whereas the normalization is a point of
981 // doubles, we cannot use the overload of operator- and must instead
982 // loop over the components of the point.
984 for (unsigned int d = 0; d < spacedim; ++d)
985 p_scaled[d] = (p[d] - normalization_shift[d]) * normalization_length;
986
987 for (unsigned int d = 0; d < spacedim; ++d)
988 result += coefficients[1 + d] * p_scaled[d];
989
990 if (!is_affine)
991 {
992 Point<dim, Number> result_affine = result;
993 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
994 for (unsigned int e = 0; e <= d; ++e, ++c)
995 result +=
996 coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
997
998 // Check if the quadratic approximation ends up considerably
999 // farther outside the unit cell on some or all SIMD lanes than
1000 // the affine approximation - in that case, we switch those
1001 // components back to the affine approximation. Note that the
1002 // quadratic approximation will grow more quickly away from the
1003 // unit cell. We make the selection for each SIMD lane with a
1004 // ternary operation.
1005 const Number distance_to_unit_cell = result.distance_square(
1007 const Number affine_distance_to_unit_cell =
1008 result_affine.distance_square(
1010 for (unsigned int d = 0; d < dim; ++d)
1011 result[d] = compare_and_apply_mask<SIMDComparison::greater_than>(
1012 distance_to_unit_cell,
1013 affine_distance_to_unit_cell + 0.5,
1014 result_affine[d],
1015 result[d]);
1016 }
1017 return result;
1018 }
1019
1020 private:
1029
1034
1038 std::array<Point<dim>, n_functions> coefficients;
1039
1046 };
1047
1048
1049
1055 template <int dim, int spacedim>
1056 inline void
1058 const CellSimilarity::Similarity cell_similarity,
1059 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1060 std::vector<Point<spacedim>> & quadrature_points,
1061 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1062 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians,
1063 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1064 {
1065 const UpdateFlags update_flags = data.update_each;
1066
1067 using VectorizedArrayType =
1068 typename ::MappingQ<dim,
1069 spacedim>::InternalData::VectorizedArrayType;
1070 const unsigned int n_shape_values = data.n_shape_functions;
1071 const unsigned int n_q_points = data.shape_info.n_q_points;
1072 constexpr unsigned int n_lanes = VectorizedArrayType::size();
1073 constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1074 constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1075
1076 // Since MappingQ::InternalData does not have separate arrays for the
1077 // covariant and contravariant transformations, but uses the arrays in
1078 // the `MappingRelatedData`, it can happen that vectors do not have the
1079 // right size
1080 if (update_flags & update_contravariant_transformation)
1081 jacobians.resize(n_q_points);
1082 if (update_flags & update_covariant_transformation)
1083 inverse_jacobians.resize(n_q_points);
1084
1085 EvaluationFlags::EvaluationFlags evaluation_flag =
1088 ((cell_similarity != CellSimilarity::translation) &&
1089 (update_flags & update_contravariant_transformation) ?
1092 ((cell_similarity != CellSimilarity::translation) &&
1093 (update_flags & update_jacobian_grads) ?
1096
1097 Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1099 Assert(!(evaluation_flag & EvaluationFlags::values) ||
1100 n_q_points == quadrature_points.size(),
1101 ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1102 Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1103 data.n_shape_functions > 0,
1105 Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
1106 n_q_points == jacobian_grads.size(),
1107 ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1108
1109 // shortcut in case we have an identity interpolation and only request
1110 // the quadrature points
1111 if (evaluation_flag == EvaluationFlags::values &&
1112 data.shape_info.element_type ==
1114 {
1115 for (unsigned int q = 0; q < n_q_points; ++q)
1116 quadrature_points[q] =
1117 data.mapping_support_points[data.shape_info
1118 .lexicographic_numbering[q]];
1119 return;
1120 }
1121
1123
1124 // prepare arrays
1125 if (evaluation_flag != EvaluationFlags::nothing)
1126 {
1127 eval.set_data_pointers(&data.scratch, n_comp);
1128
1129 // make sure to initialize on all lanes also when some are unused in
1130 // the code below
1131 for (unsigned int i = 0; i < n_shape_values * n_comp; ++i)
1132 eval.begin_dof_values()[i] = VectorizedArrayType();
1133
1134 const std::vector<unsigned int> &renumber_to_lexicographic =
1135 data.shape_info.lexicographic_numbering;
1136 for (unsigned int i = 0; i < n_shape_values; ++i)
1137 for (unsigned int d = 0; d < spacedim; ++d)
1138 {
1139 const unsigned int in_comp = d % n_lanes;
1140 const unsigned int out_comp = d / n_lanes;
1141 eval
1142 .begin_dof_values()[out_comp * n_shape_values + i][in_comp] =
1143 data.mapping_support_points[renumber_to_lexicographic[i]][d];
1144 }
1145
1146 // do the actual tensorized evaluation
1148 n_comp, evaluation_flag, eval.begin_dof_values(), eval);
1149 }
1150
1151 // do the postprocessing
1152 if (evaluation_flag & EvaluationFlags::values)
1153 {
1154 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1155 for (unsigned int i = 0; i < n_q_points; ++i)
1156 for (unsigned int in_comp = 0;
1157 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1158 ++in_comp)
1159 quadrature_points[i][out_comp * n_lanes + in_comp] =
1160 eval.begin_values()[out_comp * n_q_points + i][in_comp];
1161 }
1162
1163 if (evaluation_flag & EvaluationFlags::gradients)
1164 {
1165 // We need to reinterpret the data after evaluate has been applied.
1166 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1167 for (unsigned int point = 0; point < n_q_points; ++point)
1168 for (unsigned int j = 0; j < dim; ++j)
1169 for (unsigned int in_comp = 0;
1170 in_comp < n_lanes &&
1171 in_comp < spacedim - out_comp * n_lanes;
1172 ++in_comp)
1173 {
1174 const unsigned int total_number = point * dim + j;
1175 const unsigned int new_comp = total_number / n_q_points;
1176 const unsigned int new_point = total_number % n_q_points;
1177 jacobians[new_point][out_comp * n_lanes +
1178 in_comp][new_comp] =
1179 eval.begin_gradients()[(out_comp * n_q_points + point) *
1180 dim +
1181 j][in_comp];
1182 }
1183 }
1184
1185 if (update_flags & update_volume_elements)
1186 if (cell_similarity != CellSimilarity::translation)
1187 for (unsigned int point = 0; point < n_q_points; ++point)
1188 data.volume_elements[point] = jacobians[point].determinant();
1189
1190 // copy values from InternalData to vector given by reference
1191 if (update_flags & update_covariant_transformation)
1192 {
1193 AssertDimension(jacobians.size(), n_q_points);
1194 AssertDimension(inverse_jacobians.size(), n_q_points);
1195 if (cell_similarity != CellSimilarity::translation)
1196 for (unsigned int point = 0; point < n_q_points; ++point)
1197 inverse_jacobians[point] =
1198 jacobians[point].covariant_form().transpose();
1199 }
1200
1201 if (evaluation_flag & EvaluationFlags::hessians)
1202 {
1203 constexpr int desymmetrize_3d[6][2] = {
1204 {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1205 constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1206
1207 // We need to reinterpret the data after evaluate has been applied.
1208 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1209 for (unsigned int point = 0; point < n_q_points; ++point)
1210 for (unsigned int j = 0; j < n_hessians; ++j)
1211 for (unsigned int in_comp = 0;
1212 in_comp < n_lanes &&
1213 in_comp < spacedim - out_comp * n_lanes;
1214 ++in_comp)
1215 {
1216 const unsigned int total_number = point * n_hessians + j;
1217 const unsigned int new_point = total_number % n_q_points;
1218 const unsigned int new_hessian_comp =
1219 total_number / n_q_points;
1220 const unsigned int new_hessian_comp_i =
1221 dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1222 desymmetrize_3d[new_hessian_comp][0];
1223 const unsigned int new_hessian_comp_j =
1224 dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1225 desymmetrize_3d[new_hessian_comp][1];
1226 const double value =
1227 eval.begin_hessians()[(out_comp * n_q_points + point) *
1228 n_hessians +
1229 j][in_comp];
1230 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1231 [new_hessian_comp_i][new_hessian_comp_j] =
1232 value;
1233 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1234 [new_hessian_comp_j][new_hessian_comp_i] =
1235 value;
1236 }
1237 }
1238 }
1239
1240
1241
1242 template <int dim, int spacedim>
1243 inline DEAL_II_ALWAYS_INLINE void
1245 const unsigned int n_points,
1246 const unsigned int cur_index,
1247 const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
1248 & derivative,
1249 std::vector<DerivativeForm<1, dim, spacedim>> &result_array)
1250 {
1251 AssertDimension(result_array.size(), n_points);
1252 constexpr unsigned int n_lanes = VectorizedArray<double>::size();
1253 if (cur_index + n_lanes <= n_points)
1254 {
1255 std::array<unsigned int, n_lanes> indices;
1256 for (unsigned int j = 0; j < n_lanes; ++j)
1257 indices[j] = j * dim * spacedim;
1258 const unsigned int even = (dim * spacedim) / 4 * 4;
1259 double * result_ptr = &result_array[cur_index][0][0];
1260 const VectorizedArray<double> *derivative_ptr = &derivative[0][0];
1262 false, even, derivative_ptr, indices.data(), result_ptr);
1263 for (unsigned int d = even; d < dim * spacedim; ++d)
1264 for (unsigned int j = 0; j < n_lanes; ++j)
1265 result_ptr[j * dim * spacedim + d] = derivative_ptr[d][j];
1266 }
1267 else
1268 for (unsigned int j = 0; j < n_lanes && cur_index + j < n_points; ++j)
1269 for (unsigned int d = 0; d < spacedim; ++d)
1270 for (unsigned int e = 0; e < dim; ++e)
1271 result_array[cur_index + j][d][e] = derivative[d][e][j];
1272 }
1273
1274
1275
1276 template <int dim, int spacedim>
1277 inline void
1279 const CellSimilarity::Similarity cell_similarity,
1280 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1281 const ArrayView<const Point<dim>> & unit_points,
1282 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1283 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1284 std::vector<Point<spacedim>> & quadrature_points,
1285 std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
1286 std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
1287 {
1288 const UpdateFlags update_flags = data.update_each;
1289 const std::vector<Point<spacedim>> &support_points =
1290 data.mapping_support_points;
1291
1292 const unsigned int n_points = unit_points.size();
1293 const unsigned int n_lanes = VectorizedArray<double>::size();
1294
1295 // Since MappingQ::InternalData does not have separate arrays for the
1296 // covariant and contravariant transformations, but uses the arrays in
1297 // the `MappingRelatedData`, it can happen that vectors do not have the
1298 // right size
1299 if (update_flags & update_contravariant_transformation)
1300 jacobians.resize(n_points);
1301 if (update_flags & update_covariant_transformation)
1302 inverse_jacobians.resize(n_points);
1303
1304 const bool is_translation =
1305 cell_similarity == CellSimilarity::translation;
1306
1307 const bool needs_gradient =
1308 !is_translation &&
1309 (update_flags &
1312
1313 // Use the more heavy VectorizedArray code path if there is more than
1314 // one point left to compute
1315 for (unsigned int i = 0; i < n_points; i += n_lanes)
1316 if (n_points - i > 1)
1317 {
1319 for (unsigned int j = 0; j < n_lanes; ++j)
1320 if (i + j < n_points)
1321 for (unsigned int d = 0; d < dim; ++d)
1322 p_vec[d][j] = unit_points[i + j][d];
1323 else
1324 for (unsigned int d = 0; d < dim; ++d)
1325 p_vec[d][j] = unit_points[i][d];
1326
1329 derivative;
1330 if (needs_gradient)
1331 {
1332 const auto result =
1334 polynomials_1d,
1335 support_points,
1336 p_vec,
1337 polynomials_1d.size() == 2,
1338 renumber_lexicographic_to_hierarchic);
1339
1340 value = result.first;
1341
1342 for (unsigned int d = 0; d < spacedim; ++d)
1343 for (unsigned int e = 0; e < dim; ++e)
1344 derivative[d][e] = result.second[e][d];
1345 }
1346 else
1348 polynomials_1d,
1349 support_points,
1350 p_vec,
1351 polynomials_1d.size() == 2,
1352 renumber_lexicographic_to_hierarchic);
1353
1354 if (update_flags & update_quadrature_points)
1355 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1356 for (unsigned int d = 0; d < spacedim; ++d)
1357 quadrature_points[i + j][d] = value[d][j];
1358
1359 if (is_translation)
1360 continue;
1361
1362 if (update_flags & update_contravariant_transformation)
1363 store_vectorized_tensor(n_points, i, derivative, jacobians);
1364
1365 if (update_flags & update_volume_elements)
1366 {
1367 const auto determinant = derivative.determinant();
1368 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1369 data.volume_elements[i + j] = determinant[j];
1370 }
1371
1372 if (update_flags & update_covariant_transformation)
1373 {
1374 const auto covariant = derivative.covariant_form();
1375 store_vectorized_tensor(n_points,
1376 i,
1377 covariant.transpose(),
1378 inverse_jacobians);
1379 }
1380 }
1381 else
1382 {
1385 if (needs_gradient)
1386 {
1387 const auto result =
1389 polynomials_1d,
1390 support_points,
1391 unit_points[i],
1392 polynomials_1d.size() == 2,
1393 renumber_lexicographic_to_hierarchic);
1394
1395 value = result.first;
1396
1397 for (unsigned int d = 0; d < spacedim; ++d)
1398 for (unsigned int e = 0; e < dim; ++e)
1399 derivative[d][e] = result.second[e][d];
1400 }
1401 else
1403 polynomials_1d,
1404 support_points,
1405 unit_points[i],
1406 polynomials_1d.size() == 2,
1407 renumber_lexicographic_to_hierarchic);
1408
1409 if (update_flags & update_quadrature_points)
1410 quadrature_points[i] = value;
1411
1412 if (is_translation)
1413 continue;
1414
1415 if (dim == spacedim && update_flags & update_volume_elements)
1416 data.volume_elements[i] = derivative.determinant();
1417
1418 if (update_flags & update_contravariant_transformation)
1419 jacobians[i] = derivative;
1420
1421 if (update_flags & update_covariant_transformation)
1422 inverse_jacobians[i] = jacobians[i].covariant_form().transpose();
1423 }
1424 }
1425
1426
1427
1434 template <int dim, int spacedim>
1435 inline void
1437 const CellSimilarity::Similarity cell_similarity,
1438 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1439 const ArrayView<const Point<dim>> & unit_points,
1440 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1441 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1442 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1443 {
1444 if (data.update_each & update_jacobian_grads)
1445 {
1446 const std::vector<Point<spacedim>> &support_points =
1447 data.mapping_support_points;
1448 const unsigned int n_q_points = jacobian_grads.size();
1449
1450 if (cell_similarity != CellSimilarity::translation)
1451 for (unsigned int point = 0; point < n_q_points; ++point)
1452 {
1455 polynomials_1d,
1456 support_points,
1457 unit_points[point],
1458 renumber_lexicographic_to_hierarchic);
1459
1460 for (unsigned int i = 0; i < spacedim; ++i)
1461 for (unsigned int j = 0; j < dim; ++j)
1462 for (unsigned int l = 0; l < dim; ++l)
1463 jacobian_grads[point][i][j][l] = second[j][l][i];
1464 }
1465 }
1466 }
1467
1468
1469
1476 template <int dim, int spacedim>
1477 inline void
1479 const CellSimilarity::Similarity cell_similarity,
1480 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1481 const ArrayView<const Point<dim>> & unit_points,
1482 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1483 const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
1484 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1485 {
1486 if (data.update_each & update_jacobian_pushed_forward_grads)
1487 {
1488 const std::vector<Point<spacedim>> &support_points =
1489 data.mapping_support_points;
1490 const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
1491
1492 if (cell_similarity != CellSimilarity::translation)
1493 {
1494 double tmp[spacedim][spacedim][spacedim];
1495 for (unsigned int point = 0; point < n_q_points; ++point)
1496 {
1499 polynomials_1d,
1500 support_points,
1501 unit_points[point],
1502 renumber_lexicographic_to_hierarchic);
1503 const DerivativeForm<1, dim, spacedim> covariant =
1504 data.output_data->inverse_jacobians[point].transpose();
1505
1506 // first push forward the j-components
1507 for (unsigned int i = 0; i < spacedim; ++i)
1508 for (unsigned int j = 0; j < spacedim; ++j)
1509 for (unsigned int l = 0; l < dim; ++l)
1510 {
1511 tmp[i][j][l] = second[0][l][i] * covariant[j][0];
1512 for (unsigned int jr = 1; jr < dim; ++jr)
1513 {
1514 tmp[i][j][l] +=
1515 second[jr][l][i] * covariant[j][jr];
1516 }
1517 }
1518
1519 // now, pushing forward the l-components
1520 for (unsigned int i = 0; i < spacedim; ++i)
1521 for (unsigned int j = 0; j < spacedim; ++j)
1522 for (unsigned int l = 0; l < spacedim; ++l)
1523 {
1524 jacobian_pushed_forward_grads[point][i][j][l] =
1525 tmp[i][j][0] * covariant[l][0];
1526 for (unsigned int lr = 1; lr < dim; ++lr)
1527 {
1528 jacobian_pushed_forward_grads[point][i][j][l] +=
1529 tmp[i][j][lr] * covariant[l][lr];
1530 }
1531 }
1532 }
1533 }
1534 }
1535 }
1536
1537
1538
1539 template <int dim, int spacedim, int length_tensor>
1542 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1543 {
1544 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1546 for (unsigned int i = 0; i < spacedim; ++i)
1547 {
1548 if (dim == 1)
1549 result[i][0][0][0] = compressed[0][i];
1550 else if (dim >= 2)
1551 {
1552 for (unsigned int d = 0; d < 2; ++d)
1553 for (unsigned int e = 0; e < 2; ++e)
1554 for (unsigned int f = 0; f < 2; ++f)
1555 result[i][d][e][f] = compressed[d + e + f][i];
1556 if (dim == 3)
1557 {
1558 AssertDimension(length_tensor, 10);
1559
1560 // the derivatives are ordered in rows by increasing z
1561 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1562 // j walks through the columns
1563 for (unsigned int d = 0; d < 2; ++d)
1564 for (unsigned int e = 0; e < 2; ++e)
1565 {
1566 result[i][d][e][2] = compressed[4 + d + e][i];
1567 result[i][d][2][e] = compressed[4 + d + e][i];
1568 result[i][2][d][e] = compressed[4 + d + e][i];
1569 }
1570 for (unsigned int d = 0; d < 2; ++d)
1571 {
1572 result[i][d][2][2] = compressed[7 + d][i];
1573 result[i][2][d][2] = compressed[7 + d][i];
1574 result[i][2][2][d] = compressed[7 + d][i];
1575 }
1576 result[i][2][2][2] = compressed[9][i];
1577 }
1578 }
1579 }
1580 return result;
1581 }
1582
1583
1584
1591 template <int dim, int spacedim>
1592 inline void
1594 const CellSimilarity::Similarity cell_similarity,
1595 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1596 const ArrayView<const Point<dim>> & unit_points,
1597 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1598 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1599 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1600 {
1601 if (data.update_each & update_jacobian_2nd_derivatives)
1602 {
1603 const std::vector<Point<spacedim>> &support_points =
1604 data.mapping_support_points;
1605 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1606
1607 if (cell_similarity != CellSimilarity::translation)
1608 {
1609 for (unsigned int point = 0; point < n_q_points; ++point)
1610 {
1611 jacobian_2nd_derivatives[point] = expand_3rd_derivatives<dim>(
1612 internal::evaluate_tensor_product_higher_derivatives<3>(
1613 polynomials_1d,
1614 support_points,
1615 unit_points[point],
1616 renumber_lexicographic_to_hierarchic));
1617 }
1618 }
1619 }
1620 }
1621
1622
1623
1631 template <int dim, int spacedim>
1632 inline void
1634 const CellSimilarity::Similarity cell_similarity,
1635 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1636 const ArrayView<const Point<dim>> & unit_points,
1637 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1638 const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
1639 std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
1640 {
1642 {
1643 const std::vector<Point<spacedim>> &support_points =
1644 data.mapping_support_points;
1645 const unsigned int n_q_points =
1646 jacobian_pushed_forward_2nd_derivatives.size();
1647
1648 if (cell_similarity != CellSimilarity::translation)
1649 {
1652 for (unsigned int point = 0; point < n_q_points; ++point)
1653 {
1655 expand_3rd_derivatives<dim>(
1656 internal::evaluate_tensor_product_higher_derivatives<3>(
1657 polynomials_1d,
1658 support_points,
1659 unit_points[point],
1660 renumber_lexicographic_to_hierarchic));
1661 const DerivativeForm<1, dim, spacedim> covariant =
1662 data.output_data->inverse_jacobians[point].transpose();
1663
1664 // push forward the j-coordinate
1665 for (unsigned int i = 0; i < spacedim; ++i)
1666 for (unsigned int j = 0; j < spacedim; ++j)
1667 for (unsigned int l = 0; l < dim; ++l)
1668 for (unsigned int m = 0; m < dim; ++m)
1669 {
1670 tmp[i][j][l][m] =
1671 third[i][0][l][m] * covariant[j][0];
1672 for (unsigned int jr = 1; jr < dim; ++jr)
1673 tmp[i][j][l][m] +=
1674 third[i][jr][l][m] * covariant[j][jr];
1675 }
1676
1677 // push forward the l-coordinate
1678 for (unsigned int i = 0; i < spacedim; ++i)
1679 for (unsigned int j = 0; j < spacedim; ++j)
1680 for (unsigned int l = 0; l < spacedim; ++l)
1681 for (unsigned int m = 0; m < dim; ++m)
1682 {
1683 tmp2[i][j][l][m] =
1684 tmp[i][j][0][m] * covariant[l][0];
1685 for (unsigned int lr = 1; lr < dim; ++lr)
1686 tmp2[i][j][l][m] +=
1687 tmp[i][j][lr][m] * covariant[l][lr];
1688 }
1689
1690 // push forward the m-coordinate
1691 for (unsigned int i = 0; i < spacedim; ++i)
1692 for (unsigned int j = 0; j < spacedim; ++j)
1693 for (unsigned int l = 0; l < spacedim; ++l)
1694 for (unsigned int m = 0; m < spacedim; ++m)
1695 {
1696 jacobian_pushed_forward_2nd_derivatives
1697 [point][i][j][l][m] =
1698 tmp2[i][j][l][0] * covariant[m][0];
1699 for (unsigned int mr = 1; mr < dim; ++mr)
1700 jacobian_pushed_forward_2nd_derivatives[point][i]
1701 [j][l]
1702 [m] +=
1703 tmp2[i][j][l][mr] * covariant[m][mr];
1704 }
1705 }
1706 }
1707 }
1708 }
1709
1710
1711
1712 template <int dim, int spacedim, int length_tensor>
1715 const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
1716 {
1717 Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
1719 for (unsigned int i = 0; i < spacedim; ++i)
1720 {
1721 if (dim == 1)
1722 result[i][0][0][0][0] = compressed[0][i];
1723 else if (dim >= 2)
1724 {
1725 for (unsigned int d = 0; d < 2; ++d)
1726 for (unsigned int e = 0; e < 2; ++e)
1727 for (unsigned int f = 0; f < 2; ++f)
1728 for (unsigned int g = 0; g < 2; ++g)
1729 result[i][d][e][f][g] = compressed[d + e + f + g][i];
1730 if (dim == 3)
1731 {
1732 AssertDimension(length_tensor, 15);
1733
1734 // the derivatives are ordered in rows by increasing z
1735 // derivative, and in each row we have x^{(n-j)} y^{(j)} as
1736 // j walks through the columns
1737 for (unsigned int d = 0; d < 2; ++d)
1738 for (unsigned int e = 0; e < 2; ++e)
1739 for (unsigned int f = 0; f < 2; ++f)
1740 {
1741 result[i][d][e][f][2] = compressed[5 + d + e + f][i];
1742 result[i][d][e][2][f] = compressed[5 + d + e + f][i];
1743 result[i][d][2][e][f] = compressed[5 + d + e + f][i];
1744 result[i][2][d][e][f] = compressed[5 + d + e + f][i];
1745 }
1746 for (unsigned int d = 0; d < 2; ++d)
1747 for (unsigned int e = 0; e < 2; ++e)
1748 {
1749 result[i][d][e][2][2] = compressed[9 + d + e][i];
1750 result[i][d][2][e][2] = compressed[9 + d + e][i];
1751 result[i][d][2][2][e] = compressed[9 + d + e][i];
1752 result[i][2][d][e][2] = compressed[9 + d + e][i];
1753 result[i][2][d][2][e] = compressed[9 + d + e][i];
1754 result[i][2][2][d][e] = compressed[9 + d + e][i];
1755 }
1756 for (unsigned int d = 0; d < 2; ++d)
1757 {
1758 result[i][d][2][2][2] = compressed[12 + d][i];
1759 result[i][2][d][2][2] = compressed[12 + d][i];
1760 result[i][2][2][d][2] = compressed[12 + d][i];
1761 result[i][2][2][2][d] = compressed[12 + d][i];
1762 }
1763 result[i][2][2][2][2] = compressed[14][i];
1764 }
1765 }
1766 }
1767 return result;
1768 }
1769
1770
1771
1778 template <int dim, int spacedim>
1779 inline void
1781 const CellSimilarity::Similarity cell_similarity,
1782 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1783 const ArrayView<const Point<dim>> & unit_points,
1784 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1785 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
1786 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1787 {
1788 if (data.update_each & update_jacobian_3rd_derivatives)
1789 {
1790 const std::vector<Point<spacedim>> &support_points =
1791 data.mapping_support_points;
1792 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1793
1794 if (cell_similarity != CellSimilarity::translation)
1795 {
1796 for (unsigned int point = 0; point < n_q_points; ++point)
1797 {
1798 jacobian_3rd_derivatives[point] = expand_4th_derivatives<dim>(
1799 internal::evaluate_tensor_product_higher_derivatives<4>(
1800 polynomials_1d,
1801 support_points,
1802 unit_points[point],
1803 renumber_lexicographic_to_hierarchic));
1804 }
1805 }
1806 }
1807 }
1808
1809
1810
1818 template <int dim, int spacedim>
1819 inline void
1821 const CellSimilarity::Similarity cell_similarity,
1822 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1823 const ArrayView<const Point<dim>> & unit_points,
1824 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
1825 const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
1826 std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
1827 {
1829 {
1830 const std::vector<Point<spacedim>> &support_points =
1831 data.mapping_support_points;
1832 const unsigned int n_q_points =
1833 jacobian_pushed_forward_3rd_derivatives.size();
1834
1835 if (cell_similarity != CellSimilarity::translation)
1836 {
1837 ::
1838 ndarray<double, spacedim, spacedim, spacedim, spacedim, dim>
1839 tmp;
1841 tmp2;
1842 for (unsigned int point = 0; point < n_q_points; ++point)
1843 {
1845 expand_4th_derivatives<dim>(
1846 internal::evaluate_tensor_product_higher_derivatives<4>(
1847 polynomials_1d,
1848 support_points,
1849 unit_points[point],
1850 renumber_lexicographic_to_hierarchic));
1851
1852 const DerivativeForm<1, dim, spacedim> covariant =
1853 data.output_data->inverse_jacobians[point].transpose();
1854
1855 // push-forward the j-coordinate
1856 for (unsigned int i = 0; i < spacedim; ++i)
1857 for (unsigned int j = 0; j < spacedim; ++j)
1858 for (unsigned int l = 0; l < dim; ++l)
1859 for (unsigned int m = 0; m < dim; ++m)
1860 for (unsigned int n = 0; n < dim; ++n)
1861 {
1862 tmp[i][j][l][m][n] =
1863 fourth[i][0][l][m][n] * covariant[j][0];
1864 for (unsigned int jr = 1; jr < dim; ++jr)
1865 tmp[i][j][l][m][n] +=
1866 fourth[i][jr][l][m][n] * covariant[j][jr];
1867 }
1868
1869 // push-forward the l-coordinate
1870 for (unsigned int i = 0; i < spacedim; ++i)
1871 for (unsigned int j = 0; j < spacedim; ++j)
1872 for (unsigned int l = 0; l < spacedim; ++l)
1873 for (unsigned int m = 0; m < dim; ++m)
1874 for (unsigned int n = 0; n < dim; ++n)
1875 {
1876 tmp2[i][j][l][m][n] =
1877 tmp[i][j][0][m][n] * covariant[l][0];
1878 for (unsigned int lr = 1; lr < dim; ++lr)
1879 tmp2[i][j][l][m][n] +=
1880 tmp[i][j][lr][m][n] * covariant[l][lr];
1881 }
1882
1883 // push-forward the m-coordinate
1884 for (unsigned int i = 0; i < spacedim; ++i)
1885 for (unsigned int j = 0; j < spacedim; ++j)
1886 for (unsigned int l = 0; l < spacedim; ++l)
1887 for (unsigned int m = 0; m < spacedim; ++m)
1888 for (unsigned int n = 0; n < dim; ++n)
1889 {
1890 tmp[i][j][l][m][n] =
1891 tmp2[i][j][l][0][n] * covariant[m][0];
1892 for (unsigned int mr = 1; mr < dim; ++mr)
1893 tmp[i][j][l][m][n] +=
1894 tmp2[i][j][l][mr][n] * covariant[m][mr];
1895 }
1896
1897 // push-forward the n-coordinate
1898 for (unsigned int i = 0; i < spacedim; ++i)
1899 for (unsigned int j = 0; j < spacedim; ++j)
1900 for (unsigned int l = 0; l < spacedim; ++l)
1901 for (unsigned int m = 0; m < spacedim; ++m)
1902 for (unsigned int n = 0; n < spacedim; ++n)
1903 {
1904 jacobian_pushed_forward_3rd_derivatives
1905 [point][i][j][l][m][n] =
1906 tmp[i][j][l][m][0] * covariant[n][0];
1907 for (unsigned int nr = 1; nr < dim; ++nr)
1908 jacobian_pushed_forward_3rd_derivatives[point]
1909 [i][j][l]
1910 [m][n] +=
1911 tmp[i][j][l][m][nr] * covariant[n][nr];
1912 }
1913 }
1914 }
1915 }
1916 }
1917
1918
1919
1929 template <int dim, int spacedim>
1930 inline void
1932 const ::MappingQ<dim, spacedim> &mapping,
1933 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
1934 const unsigned int face_no,
1935 const unsigned int subface_no,
1936 const unsigned int n_q_points,
1937 const std::vector<double> & weights,
1938 const typename ::MappingQ<dim, spacedim>::InternalData &data,
1940 &output_data)
1941 {
1942 const UpdateFlags update_flags = data.update_each;
1943
1944 if (update_flags &
1946 {
1947 if (update_flags & update_boundary_forms)
1948 AssertDimension(output_data.boundary_forms.size(), n_q_points);
1949 if (update_flags & update_normal_vectors)
1950 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1951 if (update_flags & update_JxW_values)
1952 AssertDimension(output_data.JxW_values.size(), n_q_points);
1953
1954 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1955
1956 // first compute some common data that is used for evaluating
1957 // all of the flags below
1958
1959 // map the unit tangentials to the real cell. checking for d!=dim-1
1960 // eliminates compiler warnings regarding unsigned int expressions <
1961 // 0.
1962 for (unsigned int d = 0; d != dim - 1; ++d)
1963 {
1964 const unsigned int vector_index =
1966 Assert(vector_index < data.unit_tangentials.size(),
1968 Assert(data.aux[d].size() <=
1969 data.unit_tangentials[vector_index].size(),
1971 mapping.transform(make_array_view(
1972 data.unit_tangentials[vector_index]),
1974 data,
1975 make_array_view(data.aux[d]));
1976 }
1977
1978 if (update_flags & update_boundary_forms)
1979 {
1980 // if dim==spacedim, we can use the unit tangentials to compute
1981 // the boundary form by simply taking the cross product
1982 if (dim == spacedim)
1983 {
1984 for (unsigned int i = 0; i < n_q_points; ++i)
1985 switch (dim)
1986 {
1987 case 1:
1988 // in 1d, we don't have access to any of the
1989 // data.aux fields (because it has only dim-1
1990 // components), but we can still compute the
1991 // boundary form by simply looking at the number of
1992 // the face
1993 output_data.boundary_forms[i][0] =
1994 (face_no == 0 ? -1 : +1);
1995 break;
1996 case 2:
1997 output_data.boundary_forms[i] =
1998 cross_product_2d(data.aux[0][i]);
1999 break;
2000 case 3:
2001 output_data.boundary_forms[i] =
2002 cross_product_3d(data.aux[0][i], data.aux[1][i]);
2003 break;
2004 default:
2005 Assert(false, ExcNotImplemented());
2006 }
2007 }
2008 else //(dim < spacedim)
2009 {
2010 // in the codim-one case, the boundary form results from the
2011 // cross product of all the face tangential vectors and the
2012 // cell normal vector
2013 //
2014 // to compute the cell normal, use the same method used in
2015 // fill_fe_values for cells above
2016 AssertDimension(data.output_data->jacobians.size(),
2017 n_q_points);
2018
2019 for (unsigned int point = 0; point < n_q_points; ++point)
2020 {
2021 const DerivativeForm<1, dim, spacedim> contravariant =
2022 data.output_data->jacobians[point];
2023 if (dim == 1)
2024 {
2025 // J is a tangent vector
2026 output_data.boundary_forms[point] =
2027 contravariant.transpose()[0];
2028 output_data.boundary_forms[point] /=
2029 (face_no == 0 ? -1. : +1.) *
2030 output_data.boundary_forms[point].norm();
2031 }
2032
2033 if (dim == 2)
2034 {
2036 contravariant.transpose();
2037
2038 Tensor<1, spacedim> cell_normal =
2039 cross_product_3d(DX_t[0], DX_t[1]);
2040 cell_normal /= cell_normal.norm();
2041
2042 // then compute the face normal from the face
2043 // tangent and the cell normal:
2044 output_data.boundary_forms[point] =
2045 cross_product_3d(data.aux[0][point], cell_normal);
2046 }
2047 }
2048 }
2049 }
2050
2051 if (update_flags & update_JxW_values)
2052 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
2053 {
2054 output_data.JxW_values[i] =
2055 output_data.boundary_forms[i].norm() * weights[i];
2056
2057 if (subface_no != numbers::invalid_unsigned_int)
2058 {
2059 const double area_ratio = GeometryInfo<dim>::subface_ratio(
2060 cell->subface_case(face_no), subface_no);
2061 output_data.JxW_values[i] *= area_ratio;
2062 }
2063 }
2064
2065 if (update_flags & update_normal_vectors)
2066 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
2067 output_data.normal_vectors[i] =
2068 Point<spacedim>(output_data.boundary_forms[i] /
2069 output_data.boundary_forms[i].norm());
2070 }
2071 }
2072
2073
2074
2081 template <int dim, int spacedim>
2082 inline void
2084 const ::MappingQ<dim, spacedim> &mapping,
2085 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
2086 const unsigned int face_no,
2087 const unsigned int subface_no,
2088 const typename QProjector<dim>::DataSetDescriptor data_set,
2089 const Quadrature<dim - 1> & quadrature,
2090 const typename ::MappingQ<dim, spacedim>::InternalData &data,
2091 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
2092 const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
2094 &output_data)
2095 {
2096 const ArrayView<const Point<dim>> quadrature_points(
2097 &data.quadrature_points[data_set], quadrature.size());
2098 if (dim > 1 && data.tensor_product_quadrature)
2099 {
2100 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2102 data,
2103 output_data.quadrature_points,
2104 output_data.jacobians,
2105 output_data.inverse_jacobians,
2106 output_data.jacobian_grads);
2107 }
2108 else
2109 {
2113 data,
2114 quadrature_points,
2115 polynomials_1d,
2116 renumber_lexicographic_to_hierarchic,
2117 output_data.quadrature_points,
2118 output_data.jacobians,
2119 output_data.inverse_jacobians);
2120 maybe_update_jacobian_grads<dim, spacedim>(
2122 data,
2123 quadrature_points,
2124 polynomials_1d,
2125 renumber_lexicographic_to_hierarchic,
2126 output_data.jacobian_grads);
2127 }
2128 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2130 data,
2131 quadrature_points,
2132 polynomials_1d,
2133 renumber_lexicographic_to_hierarchic,
2134 output_data.jacobian_pushed_forward_grads);
2135 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2137 data,
2138 quadrature_points,
2139 polynomials_1d,
2140 renumber_lexicographic_to_hierarchic,
2141 output_data.jacobian_2nd_derivatives);
2142 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2144 data,
2145 quadrature_points,
2146 polynomials_1d,
2147 renumber_lexicographic_to_hierarchic,
2149 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2151 data,
2152 quadrature_points,
2153 polynomials_1d,
2154 renumber_lexicographic_to_hierarchic,
2155 output_data.jacobian_3rd_derivatives);
2156 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2158 data,
2159 quadrature_points,
2160 polynomials_1d,
2161 renumber_lexicographic_to_hierarchic,
2163
2165 cell,
2166 face_no,
2167 subface_no,
2168 quadrature.size(),
2169 quadrature.get_weights(),
2170 data,
2171 output_data);
2172 }
2173
2174
2175
2179 template <int dim, int spacedim, int rank>
2180 inline void
2182 const ArrayView<const Tensor<rank, dim>> & input,
2183 const MappingKind mapping_kind,
2184 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2185 const ArrayView<Tensor<rank, spacedim>> & output)
2186 {
2187 AssertDimension(input.size(), output.size());
2188 Assert((dynamic_cast<
2189 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2190 &mapping_data) != nullptr),
2192 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2193 static_cast<
2194 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2195 mapping_data);
2196
2197 switch (mapping_kind)
2198 {
2200 {
2203 "update_contravariant_transformation"));
2204
2205 for (unsigned int i = 0; i < output.size(); ++i)
2206 output[i] = apply_transformation(data.output_data->jacobians[i],
2207 input[i]);
2208
2209 return;
2210 }
2211
2212 case mapping_piola:
2213 {
2216 "update_contravariant_transformation"));
2217 Assert(data.update_each & update_volume_elements,
2219 "update_volume_elements"));
2220 Assert(rank == 1, ExcMessage("Only for rank 1"));
2221 if (rank != 1)
2222 return;
2223
2224 for (unsigned int i = 0; i < output.size(); ++i)
2225 {
2226 output[i] =
2227 apply_transformation(data.output_data->jacobians[i],
2228 input[i]);
2229 output[i] /= data.volume_elements[i];
2230 }
2231 return;
2232 }
2233 // We still allow this operation as in the
2234 // reference cell Derivatives are Tensor
2235 // rather than DerivativeForm
2236 case mapping_covariant:
2237 {
2240 "update_covariant_transformation"));
2241
2242 for (unsigned int i = 0; i < output.size(); ++i)
2243 output[i] = apply_transformation(
2244 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2245
2246 return;
2247 }
2248
2249 default:
2250 Assert(false, ExcNotImplemented());
2251 }
2252 }
2253
2254
2255
2259 template <int dim, int spacedim, int rank>
2260 inline void
2262 const ArrayView<const Tensor<rank, dim>> & input,
2263 const MappingKind mapping_kind,
2264 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2265 const ArrayView<Tensor<rank, spacedim>> & output)
2266 {
2267 AssertDimension(input.size(), output.size());
2268 Assert((dynamic_cast<
2269 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2270 &mapping_data) != nullptr),
2272 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2273 static_cast<
2274 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2275 mapping_data);
2276
2277 switch (mapping_kind)
2278 {
2280 {
2281 Assert(data.update_each & update_covariant_transformation,
2283 "update_covariant_transformation"));
2286 "update_contravariant_transformation"));
2287 Assert(rank == 2, ExcMessage("Only for rank 2"));
2288
2289 for (unsigned int i = 0; i < output.size(); ++i)
2290 {
2292 apply_transformation(data.output_data->jacobians[i],
2293 transpose(input[i]));
2294 output[i] = apply_transformation(
2295 data.output_data->inverse_jacobians[i].transpose(),
2296 A.transpose());
2297 }
2298
2299 return;
2300 }
2301
2303 {
2304 Assert(data.update_each & update_covariant_transformation,
2306 "update_covariant_transformation"));
2307 Assert(rank == 2, ExcMessage("Only for rank 2"));
2308
2309 for (unsigned int i = 0; i < output.size(); ++i)
2310 {
2311 const DerivativeForm<1, dim, spacedim> covariant =
2312 data.output_data->inverse_jacobians[i].transpose();
2314 apply_transformation(covariant, transpose(input[i]));
2315 output[i] = apply_transformation(covariant, A.transpose());
2316 }
2317
2318 return;
2319 }
2320
2322 {
2323 Assert(data.update_each & update_covariant_transformation,
2325 "update_covariant_transformation"));
2328 "update_contravariant_transformation"));
2329 Assert(data.update_each & update_volume_elements,
2331 "update_volume_elements"));
2332 Assert(rank == 2, ExcMessage("Only for rank 2"));
2333
2334 for (unsigned int i = 0; i < output.size(); ++i)
2335 {
2336 const DerivativeForm<1, dim, spacedim> covariant =
2337 data.output_data->inverse_jacobians[i].transpose();
2339 apply_transformation(covariant, input[i]);
2340 const Tensor<2, spacedim> T =
2341 apply_transformation(data.output_data->jacobians[i],
2342 A.transpose());
2343
2344 output[i] = transpose(T);
2345 output[i] /= data.volume_elements[i];
2346 }
2347
2348 return;
2349 }
2350
2351 default:
2352 Assert(false, ExcNotImplemented());
2353 }
2354 }
2355
2356
2357
2361 template <int dim, int spacedim>
2362 inline void
2364 const ArrayView<const Tensor<3, dim>> & input,
2365 const MappingKind mapping_kind,
2366 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2367 const ArrayView<Tensor<3, spacedim>> & output)
2368 {
2369 AssertDimension(input.size(), output.size());
2370 Assert((dynamic_cast<
2371 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2372 &mapping_data) != nullptr),
2374 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2375 static_cast<
2376 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2377 mapping_data);
2378
2379 switch (mapping_kind)
2380 {
2382 {
2383 Assert(data.update_each & update_covariant_transformation,
2385 "update_covariant_transformation"));
2388 "update_contravariant_transformation"));
2389
2390 for (unsigned int q = 0; q < output.size(); ++q)
2391 {
2392 const DerivativeForm<1, dim, spacedim> covariant =
2393 data.output_data->inverse_jacobians[q].transpose();
2394 const DerivativeForm<1, dim, spacedim> contravariant =
2395 data.output_data->jacobians[q];
2396
2397 for (unsigned int i = 0; i < spacedim; ++i)
2398 {
2399 double tmp1[dim][dim];
2400 for (unsigned int J = 0; J < dim; ++J)
2401 for (unsigned int K = 0; K < dim; ++K)
2402 {
2403 tmp1[J][K] =
2404 contravariant[i][0] * input[q][0][J][K];
2405 for (unsigned int I = 1; I < dim; ++I)
2406 tmp1[J][K] +=
2407 contravariant[i][I] * input[q][I][J][K];
2408 }
2409 for (unsigned int j = 0; j < spacedim; ++j)
2410 {
2411 double tmp2[dim];
2412 for (unsigned int K = 0; K < dim; ++K)
2413 {
2414 tmp2[K] = covariant[j][0] * tmp1[0][K];
2415 for (unsigned int J = 1; J < dim; ++J)
2416 tmp2[K] += covariant[j][J] * tmp1[J][K];
2417 }
2418 for (unsigned int k = 0; k < spacedim; ++k)
2419 {
2420 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2421 for (unsigned int K = 1; K < dim; ++K)
2422 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2423 }
2424 }
2425 }
2426 }
2427 return;
2428 }
2429
2431 {
2432 Assert(data.update_each & update_covariant_transformation,
2434 "update_covariant_transformation"));
2435
2436 for (unsigned int q = 0; q < output.size(); ++q)
2437 {
2438 const DerivativeForm<1, dim, spacedim> covariant =
2439 data.output_data->inverse_jacobians[q].transpose();
2440
2441 for (unsigned int i = 0; i < spacedim; ++i)
2442 {
2443 double tmp1[dim][dim];
2444 for (unsigned int J = 0; J < dim; ++J)
2445 for (unsigned int K = 0; K < dim; ++K)
2446 {
2447 tmp1[J][K] = covariant[i][0] * input[q][0][J][K];
2448 for (unsigned int I = 1; I < dim; ++I)
2449 tmp1[J][K] += covariant[i][I] * input[q][I][J][K];
2450 }
2451 for (unsigned int j = 0; j < spacedim; ++j)
2452 {
2453 double tmp2[dim];
2454 for (unsigned int K = 0; K < dim; ++K)
2455 {
2456 tmp2[K] = covariant[j][0] * tmp1[0][K];
2457 for (unsigned int J = 1; J < dim; ++J)
2458 tmp2[K] += covariant[j][J] * tmp1[J][K];
2459 }
2460 for (unsigned int k = 0; k < spacedim; ++k)
2461 {
2462 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2463 for (unsigned int K = 1; K < dim; ++K)
2464 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2465 }
2466 }
2467 }
2468 }
2469
2470 return;
2471 }
2472
2474 {
2475 Assert(data.update_each & update_covariant_transformation,
2477 "update_covariant_transformation"));
2480 "update_contravariant_transformation"));
2481 Assert(data.update_each & update_volume_elements,
2483 "update_volume_elements"));
2484
2485 for (unsigned int q = 0; q < output.size(); ++q)
2486 {
2487 const DerivativeForm<1, dim, spacedim> covariant =
2488 data.output_data->inverse_jacobians[q].transpose();
2489 const DerivativeForm<1, dim, spacedim> contravariant =
2490 data.output_data->jacobians[q];
2491 for (unsigned int i = 0; i < spacedim; ++i)
2492 {
2493 double factor[dim];
2494 for (unsigned int I = 0; I < dim; ++I)
2495 factor[I] =
2496 contravariant[i][I] * (1. / data.volume_elements[q]);
2497 double tmp1[dim][dim];
2498 for (unsigned int J = 0; J < dim; ++J)
2499 for (unsigned int K = 0; K < dim; ++K)
2500 {
2501 tmp1[J][K] = factor[0] * input[q][0][J][K];
2502 for (unsigned int I = 1; I < dim; ++I)
2503 tmp1[J][K] += factor[I] * input[q][I][J][K];
2504 }
2505 for (unsigned int j = 0; j < spacedim; ++j)
2506 {
2507 double tmp2[dim];
2508 for (unsigned int K = 0; K < dim; ++K)
2509 {
2510 tmp2[K] = covariant[j][0] * tmp1[0][K];
2511 for (unsigned int J = 1; J < dim; ++J)
2512 tmp2[K] += covariant[j][J] * tmp1[J][K];
2513 }
2514 for (unsigned int k = 0; k < spacedim; ++k)
2515 {
2516 output[q][i][j][k] = covariant[k][0] * tmp2[0];
2517 for (unsigned int K = 1; K < dim; ++K)
2518 output[q][i][j][k] += covariant[k][K] * tmp2[K];
2519 }
2520 }
2521 }
2522 }
2523
2524 return;
2525 }
2526
2527 default:
2528 Assert(false, ExcNotImplemented());
2529 }
2530 }
2531
2532
2533
2538 template <int dim, int spacedim, int rank>
2539 inline void
2542 const MappingKind mapping_kind,
2543 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_data,
2544 const ArrayView<Tensor<rank + 1, spacedim>> & output)
2545 {
2546 AssertDimension(input.size(), output.size());
2547 Assert((dynamic_cast<
2548 const typename ::MappingQ<dim, spacedim>::InternalData *>(
2549 &mapping_data) != nullptr),
2551 const typename ::MappingQ<dim, spacedim>::InternalData &data =
2552 static_cast<
2553 const typename ::MappingQ<dim, spacedim>::InternalData &>(
2554 mapping_data);
2555
2556 switch (mapping_kind)
2557 {
2558 case mapping_covariant:
2559 {
2562 "update_covariant_transformation"));
2563
2564 for (unsigned int i = 0; i < output.size(); ++i)
2565 output[i] = apply_transformation(
2566 data.output_data->inverse_jacobians[i].transpose(), input[i]);
2567
2568 return;
2569 }
2570 default:
2571 Assert(false, ExcNotImplemented());
2572 }
2573 }
2574 } // namespace MappingQImplementation
2575} // namespace internal
2576
2578
2579#endif
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
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:317
Definition point.h:112
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
const Point< dim > & point(const unsigned int i) 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
InverseQuadraticApproximation(const InverseQuadraticApproximation &)=default
InverseQuadraticApproximation(const std::vector< Point< spacedim > > &real_support_points, const std::vector< Point< dim > > &unit_support_points)
Point< dim, Number > compute(const Point< spacedim, Number > &p) const
#define DEAL_II_ALWAYS_INLINE
Definition config.h:106
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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:4616
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:78
@ mapping_piola
Definition mapping.h:113
@ mapping_covariant_gradient
Definition mapping.h:99
@ mapping_covariant
Definition mapping.h:88
@ mapping_contravariant
Definition mapping.h:93
@ mapping_contravariant_hessian
Definition mapping.h:155
@ mapping_covariant_hessian
Definition mapping.h:149
@ mapping_contravariant_gradient
Definition mapping.h:105
@ mapping_piola_gradient
Definition mapping.h:119
@ mapping_piola_hessian
Definition mapping.h:161
LogStream deallog
Definition logstream.cc:37
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
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)
void store_vectorized_tensor(const unsigned int n_points, const unsigned int cur_index, const DerivativeForm< 1, dim, spacedim, VectorizedArray< double > > &derivative, std::vector< DerivativeForm< 1, dim, spacedim > > &result_array)
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)
Point< dim > do_transform_real_to_unit_cell_internal_codim1(const Point< dim+1 > &p, const Point< dim > &initial_p_unit, const std::vector< Point< dim+1 > > &points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber)
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)
Point< dim, Number > do_transform_real_to_unit_cell_internal(const Point< spacedim, Number > &p, const Point< dim, Number > &initial_p_unit, const std::vector< Point< spacedim > > &points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber, const bool print_iterations_to_deallog=false)
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)
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)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, 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 std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition types.h:213
::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:108
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 > &)
void vectorized_transpose_and_store(const bool add_into, const unsigned int n_entries, const VectorizedArray< Number, width > *in, const unsigned int *offsets, Number *out)