Reference documentation for deal.II version 9.3.3
\(\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\}}\)
mapping_q_internal.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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
40
41#include <array>
42
43
45
46namespace internal
47{
53 namespace MappingQ1
54 {
55 // These are left as templates on the spatial dimension (even though dim
56 // == spacedim must be true for them to make sense) because templates are
57 // expanded before the compiler eliminates code due to the 'if (dim ==
58 // spacedim)' statement (see the body of the general
59 // transform_real_to_unit_cell).
60 template <int spacedim>
61 inline Point<1>
64 & vertices,
65 const Point<spacedim> &p)
66 {
67 Assert(spacedim == 1, ExcInternalError());
68 return Point<1>((p[0] - vertices[0](0)) /
69 (vertices[1](0) - vertices[0](0)));
70 }
71
72
73
74 template <int spacedim>
75 inline Point<2>
78 & vertices,
79 const Point<spacedim> &p)
80 {
81 Assert(spacedim == 2, ExcInternalError());
82
83 // For accuracy reasons, we do all arithmetic in extended precision
84 // (long double). This has a noticeable effect on the hit rate for
85 // borderline cases and thus makes the algorithm more robust.
86 const long double x = p(0);
87 const long double y = p(1);
88
89 const long double x0 = vertices[0](0);
90 const long double x1 = vertices[1](0);
91 const long double x2 = vertices[2](0);
92 const long double x3 = vertices[3](0);
93
94 const long double y0 = vertices[0](1);
95 const long double y1 = vertices[1](1);
96 const long double y2 = vertices[2](1);
97 const long double y3 = vertices[3](1);
98
99 const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3);
100 const long double b = -(x0 - x1 - x2 + x3) * y + (x - 2 * x1 + x3) * y0 -
101 (x - 2 * x0 + x2) * y1 - (x - x1) * y2 +
102 (x - x0) * y3;
103 const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1;
104
105 const long double discriminant = b * b - 4 * a * c;
106 // exit if the point is not in the cell (this is the only case where the
107 // discriminant is negative)
109 discriminant > 0.0,
111
112 long double eta1;
113 long double eta2;
114 const long double sqrt_discriminant = std::sqrt(discriminant);
115 // special case #1: if a is near-zero to make the discriminant exactly
116 // equal b, then use the linear formula
117 if (b != 0.0 && std::abs(b) == sqrt_discriminant)
118 {
119 eta1 = -c / b;
120 eta2 = -c / b;
121 }
122 // special case #2: a is zero for parallelograms and very small for
123 // near-parallelograms:
124 else if (std::abs(a) < 1e-8 * std::abs(b))
125 {
126 // if both a and c are very small then the root should be near
127 // zero: this first case will capture that
128 eta1 = 2 * c / (-b - sqrt_discriminant);
129 eta2 = 2 * c / (-b + sqrt_discriminant);
130 }
131 // finally, use the plain version:
132 else
133 {
134 eta1 = (-b - sqrt_discriminant) / (2 * a);
135 eta2 = (-b + sqrt_discriminant) / (2 * a);
136 }
137 // pick the one closer to the center of the cell.
138 const long double eta =
139 (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
140
141 /*
142 * There are two ways to compute xi from eta, but either one may have a
143 * zero denominator.
144 */
145 const long double subexpr0 = -eta * x2 + x0 * (eta - 1);
146 const long double xi_denominator0 = eta * x3 - x1 * (eta - 1) + subexpr0;
147 const long double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
148 std::max(std::abs(x2), std::abs(x3)));
149
150 if (std::abs(xi_denominator0) > 1e-10 * max_x)
151 {
152 const double xi = (x + subexpr0) / xi_denominator0;
153 return {xi, static_cast<double>(eta)};
154 }
155 else
156 {
157 const long double max_y =
159 std::max(std::abs(y2), std::abs(y3)));
160 const long double subexpr1 = -eta * y2 + y0 * (eta - 1);
161 const long double xi_denominator1 =
162 eta * y3 - y1 * (eta - 1) + subexpr1;
163 if (std::abs(xi_denominator1) > 1e-10 * max_y)
164 {
165 const double xi = (subexpr1 + y) / xi_denominator1;
166 return {xi, static_cast<double>(eta)};
167 }
168 else // give up and try Newton iteration
169 {
171 false,
172 (typename Mapping<spacedim,
173 spacedim>::ExcTransformationFailed()));
174 }
175 }
176 // bogus return to placate compiler. It should not be possible to get
177 // here.
178 Assert(false, ExcInternalError());
179 return {std::numeric_limits<double>::quiet_NaN(),
180 std::numeric_limits<double>::quiet_NaN()};
181 }
182
183
184
185 template <int spacedim>
186 inline Point<3>
189 & /*vertices*/,
190 const Point<spacedim> & /*p*/)
191 {
192 // It should not be possible to get here
193 Assert(false, ExcInternalError());
194 return {std::numeric_limits<double>::quiet_NaN(),
195 std::numeric_limits<double>::quiet_NaN(),
196 std::numeric_limits<double>::quiet_NaN()};
197 }
198 } // namespace MappingQ1
199
200
201
207 namespace MappingQGenericImplementation
208 {
213 template <int dim>
214 std::vector<Point<dim>>
215 unit_support_points(const std::vector<Point<1>> & line_support_points,
216 const std::vector<unsigned int> &renumbering)
217 {
218 AssertDimension(Utilities::pow(line_support_points.size(), dim),
219 renumbering.size());
220 std::vector<Point<dim>> points(renumbering.size());
221 const unsigned int n1 = line_support_points.size();
222 for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2)
223 for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1)
224 for (unsigned int q0 = 0; q0 < n1; ++q0, ++q)
225 {
226 points[renumbering[q]][0] = line_support_points[q0][0];
227 if (dim > 1)
228 points[renumbering[q]][1] = line_support_points[q1][0];
229 if (dim > 2)
230 points[renumbering[q]][2] = line_support_points[q2][0];
231 }
232 return points;
233 }
234
235
236
244 inline ::Table<2, double>
245 compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
246 {
247 ::Table<2, double> loqvs;
248
249 // we are asked to compute weights for interior support points, but
250 // there are no interior points if degree==1
251 if (polynomial_degree == 1)
252 return loqvs;
253
254 const unsigned int M = polynomial_degree - 1;
255 const unsigned int n_inner_2d = M * M;
256 const unsigned int n_outer_2d = 4 + 4 * M;
257
258 // set the weights of transfinite interpolation
259 loqvs.reinit(n_inner_2d, n_outer_2d);
260 QGaussLobatto<2> gl(polynomial_degree + 1);
261 for (unsigned int i = 0; i < M; ++i)
262 for (unsigned int j = 0; j < M; ++j)
263 {
264 const Point<2> &p =
265 gl.point((i + 1) * (polynomial_degree + 1) + (j + 1));
266 const unsigned int index_table = i * M + j;
267 for (unsigned int v = 0; v < 4; ++v)
268 loqvs(index_table, v) =
270 loqvs(index_table, 4 + i) = 1. - p[0];
271 loqvs(index_table, 4 + i + M) = p[0];
272 loqvs(index_table, 4 + j + 2 * M) = 1. - p[1];
273 loqvs(index_table, 4 + j + 3 * M) = p[1];
274 }
275
276 // the sum of weights of the points at the outer rim should be one.
277 // check this
278 for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
279 Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
280 loqvs[unit_point].end(),
281 0.) -
282 1) < 1e-13 * polynomial_degree,
284
285 return loqvs;
286 }
287
288
289
296 inline ::Table<2, double>
297 compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
298 {
299 ::Table<2, double> lohvs;
300
301 // we are asked to compute weights for interior support points, but
302 // there are no interior points if degree==1
303 if (polynomial_degree == 1)
304 return lohvs;
305
306 const unsigned int M = polynomial_degree - 1;
307
308 const unsigned int n_inner = Utilities::fixed_power<3>(M);
309 const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
310
311 // set the weights of transfinite interpolation
312 lohvs.reinit(n_inner, n_outer);
313 QGaussLobatto<3> gl(polynomial_degree + 1);
314 for (unsigned int i = 0; i < M; ++i)
315 for (unsigned int j = 0; j < M; ++j)
316 for (unsigned int k = 0; k < M; ++k)
317 {
318 const Point<3> & p = gl.point((i + 1) * (M + 2) * (M + 2) +
319 (j + 1) * (M + 2) + (k + 1));
320 const unsigned int index_table = i * M * M + j * M + k;
321
322 // vertices
323 for (unsigned int v = 0; v < 8; ++v)
324 lohvs(index_table, v) =
326
327 // lines
328 {
329 constexpr std::array<unsigned int, 4> line_coordinates_y(
330 {{0, 1, 4, 5}});
331 const Point<2> py(p[0], p[2]);
332 for (unsigned int l = 0; l < 4; ++l)
333 lohvs(index_table, 8 + line_coordinates_y[l] * M + j) =
335 }
336
337 {
338 constexpr std::array<unsigned int, 4> line_coordinates_x(
339 {{2, 3, 6, 7}});
340 const Point<2> px(p[1], p[2]);
341 for (unsigned int l = 0; l < 4; ++l)
342 lohvs(index_table, 8 + line_coordinates_x[l] * M + k) =
344 }
345
346 {
347 constexpr std::array<unsigned int, 4> line_coordinates_z(
348 {{8, 9, 10, 11}});
349 const Point<2> pz(p[0], p[1]);
350 for (unsigned int l = 0; l < 4; ++l)
351 lohvs(index_table, 8 + line_coordinates_z[l] * M + i) =
353 }
354
355 // quads
356 lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
357 1. - p[0];
358 lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0];
359 lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) =
360 1. - p[1];
361 lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1];
362 lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) =
363 1. - p[2];
364 lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
365 }
366
367 // the sum of weights of the points at the outer rim should be one.
368 // check this
369 for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
370 Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
371 lohvs[unit_point].end(),
372 0.) -
373 1) < 1e-13 * polynomial_degree,
375
376 return lohvs;
377 }
378
379
380
385 inline std::vector<::Table<2, double>>
387 const unsigned int polynomial_degree,
388 const unsigned int dim)
389 {
390 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
391 std::vector<::Table<2, double>> output(dim);
392 if (polynomial_degree <= 1)
393 return output;
394
395 // fill the 1D interior weights
396 QGaussLobatto<1> quadrature(polynomial_degree + 1);
397 output[0].reinit(polynomial_degree - 1,
399 for (unsigned int q = 0; q < polynomial_degree - 1; ++q)
400 for (const unsigned int i : GeometryInfo<1>::vertex_indices())
401 output[0](q, i) =
403 i);
404
405 if (dim > 1)
406 output[1] = compute_support_point_weights_on_quad(polynomial_degree);
407
408 if (dim > 2)
409 output[2] = compute_support_point_weights_on_hex(polynomial_degree);
410
411 return output;
412 }
413
414
415
419 template <int dim>
420 inline ::Table<2, double>
421 compute_support_point_weights_cell(const unsigned int polynomial_degree)
422 {
423 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
424 if (polynomial_degree <= 1)
425 return ::Table<2, double>();
426
427 QGaussLobatto<dim> quadrature(polynomial_degree + 1);
428 const std::vector<unsigned int> h2l =
429 FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
430
431 ::Table<2, double> output(quadrature.size() -
434 for (unsigned int q = 0; q < output.size(0); ++q)
435 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
437 quadrature.point(h2l[q + GeometryInfo<dim>::vertices_per_cell]), i);
438
439 return output;
440 }
441
442
443
451 template <int dim, int spacedim>
452 inline Point<spacedim>
454 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data)
455 {
456 AssertDimension(data.shape_values.size(),
457 data.mapping_support_points.size());
458
459 // use now the InternalData to compute the point in real space.
460 Point<spacedim> p_real;
461 for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i)
462 p_real += data.mapping_support_points[i] * data.shape(0, i);
463
464 return p_real;
465 }
466
467
468
473 template <int dim, int spacedim, typename Number>
474 inline Point<dim, Number>
476 const Point<spacedim, Number> & p,
477 const Point<dim, Number> & initial_p_unit,
478 const std::vector<Point<spacedim>> & points,
479 const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
480 const std::vector<unsigned int> & renumber,
481 const bool print_iterations_to_deallog = false)
482 {
483 AssertDimension(points.size(),
484 Utilities::pow(polynomials_1d.size(), dim));
485
486 // Newton iteration to solve
487 // f(x)=p(x)-p=0
488 // where we are looking for 'x' and p(x) is the forward transformation
489 // from unit to real cell. We solve this using a Newton iteration
490 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
491 // The start value is set to be the linear approximation to the cell
492
493 // The shape values and derivatives of the mapping at this point are
494 // previously computed.
495
496 Point<dim, Number> p_unit = initial_p_unit;
498 polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber);
499
500 Tensor<1, spacedim, Number> f = p_real.first - p;
501
502 // early out if we already have our point in all SIMD lanes, i.e.,
503 // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable
504 // this step also for VectorizedArray where we do not have operator <,
505 // compare the result to zero which is defined for SIMD types
506 if (std::max(Number(0.),
507 f.norm_square() - 1e-24 * p_real.second[0].norm_square()) ==
508 Number(0.))
509 return p_unit;
510
511 // we need to compare the position of the computed p(x) against the
512 // given point 'p'. We will terminate the iteration and return 'x' if
513 // they are less than eps apart. The question is how to choose eps --
514 // or, put maybe more generally: in which norm we want these 'p' and
515 // 'p(x)' to be eps apart.
516 //
517 // the question is difficult since we may have to deal with very
518 // elongated cells where we may achieve 1e-12*h for the distance of
519 // these two points in the 'long' direction, but achieving this
520 // tolerance in the 'short' direction of the cell may not be possible
521 //
522 // what we do instead is then to terminate iterations if
523 // \| p(x) - p \|_A < eps
524 // where the A-norm is somehow induced by the transformation of the
525 // cell. in particular, we want to measure distances relative to the
526 // sizes of the cell in its principal directions.
527 //
528 // to define what exactly A should be, note that to first order we have
529 // the following (assuming that x* is the solution of the problem, i.e.,
530 // p(x*)=p):
531 // p(x) - p = p(x) - p(x*)
532 // = -grad p(x) * (x*-x) + higher order terms
533 // This suggest to measure with a norm that corresponds to
534 // A = {[grad p(x)]^T [grad p(x)]}^{-1}
535 // because then
536 // \| p(x) - p \|_A \approx \| x - x* \|
537 // Consequently, we will try to enforce that
538 // \| p(x) - p \|_A = \| f \| <= eps
539 //
540 // Note that using this norm is a bit dangerous since the norm changes
541 // in every iteration (A isn't fixed by depending on xk). However, if
542 // the cell is not too deformed (it may be stretched, but not twisted)
543 // then the mapping is almost linear and A is indeed constant or
544 // nearly so.
545 const double eps = 1.e-11;
546 const unsigned int newton_iteration_limit = 20;
547
548 Point<dim, Number> invalid_point;
549 invalid_point[0] = std::numeric_limits<double>::infinity();
550 bool try_project_to_unit_cell = false;
551
552 unsigned int newton_iteration = 0;
553 Number f_weighted_norm_square = 1.;
554 Number last_f_weighted_norm_square = 1.;
555
556 do
557 {
558 if (print_iterations_to_deallog)
559 deallog << "Newton iteration " << newton_iteration
560 << " for unit point guess " << p_unit << std::endl;
561
562 // f'(x)
564 for (unsigned int d = 0; d < spacedim; ++d)
565 for (unsigned int e = 0; e < dim; ++e)
566 df[d][e] = p_real.second[e][d];
567
568 // check determinand(df) > 0 on all SIMD lanes
569 if (!(std::min(determinant(df),
572 {
573 // We allow to enter this function with an initial guess
574 // outside the unit cell. We might have run into invalid
575 // Jacobians due to that, so we should at least try once to go
576 // back to the unit cell and go on with the Newton iteration
577 // from there. Since the outside case is unlikely, we can
578 // afford spending the extra effort at this place.
579 if (try_project_to_unit_cell == false)
580 {
583 polynomials_1d,
584 points,
585 p_unit,
586 polynomials_1d.size() == 2,
587 renumber);
588 f = p_real.first - p;
589 f_weighted_norm_square = 1.;
590 last_f_weighted_norm_square = 1;
591 try_project_to_unit_cell = true;
592 continue;
593 }
594 else
595 return invalid_point;
596 }
597
598 // Solve [f'(x)]d=f(x)
599 const Tensor<2, spacedim, Number> df_inverse = invert(df);
600 const Tensor<1, spacedim, Number> delta = df_inverse * f;
601 last_f_weighted_norm_square = delta.norm_square();
602
603 if (print_iterations_to_deallog)
604 deallog << " delta=" << delta << std::endl;
605
606 // do a line search
607 double step_length = 1;
608 do
609 {
610 // update of p_unit. The spacedim-th component of transformed
611 // point is simply ignored in codimension one case. When this
612 // component is not zero, then we are projecting the point to
613 // the surface or curve identified by the cell.
614 Point<dim, Number> p_unit_trial = p_unit;
615 for (unsigned int i = 0; i < dim; ++i)
616 p_unit_trial[i] -= step_length * delta[i];
617
618 // shape values and derivatives at new p_unit point
619 const auto p_real_trial =
621 polynomials_1d,
622 points,
623 p_unit_trial,
624 polynomials_1d.size() == 2,
625 renumber);
626 const Tensor<1, spacedim, Number> f_trial =
627 p_real_trial.first - p;
628 f_weighted_norm_square = (df_inverse * f_trial).norm_square();
629
630 if (print_iterations_to_deallog)
631 deallog << " step_length=" << step_length << std::endl
632 << " ||f || =" << f.norm() << std::endl
633 << " ||f*|| =" << f_trial.norm() << std::endl
634 << " ||f*||_A ="
635 << std::sqrt(f_weighted_norm_square) << std::endl;
636
637 // See if we are making progress with the current step length
638 // and if not, reduce it by a factor of two and try again.
639 //
640 // Strictly speaking, we should probably use the same norm as we
641 // use for the outer algorithm. In practice, line search is just
642 // a crutch to find a "reasonable" step length, and so using the
643 // l2 norm is probably just fine.
644 //
645 // check f_trial.norm() < f.norm() in SIMD form. This is a bit
646 // more complicated because some SIMD lanes might not be doing
647 // any progress any more as they have already reached roundoff
648 // accuracy. We define that as the case of updates less than
649 // 1e-6. The tolerance might seem coarse but since we are
650 // dealing with a Newton iteration of a polynomial function we
651 // either converge quadratically or not any more. Thus, our
652 // selection is to terminate if either the norm of f is
653 // decreasing or that threshold of 1e-6 is reached.
654 if (std::max(f_weighted_norm_square - 1e-6 * 1e-6, Number(0.)) *
655 std::max(f_trial.norm_square() - f.norm_square(),
656 Number(0.)) ==
657 Number(0.))
658 {
659 p_real = p_real_trial;
660 p_unit = p_unit_trial;
661 f = f_trial;
662 break;
663 }
664 else if (step_length > 0.05)
665 step_length *= 0.5;
666 else
667 break;
668 }
669 while (true);
670
671 // In case we terminated the line search due to the step becoming
672 // too small, we give the iteration another try with the
673 // projection of the initial guess to the unit cell before we give
674 // up, just like for the negative determinant case.
675 if (step_length <= 0.05 && try_project_to_unit_cell == false)
676 {
679 polynomials_1d,
680 points,
681 p_unit,
682 polynomials_1d.size() == 2,
683 renumber);
684 f = p_real.first - p;
685 f_weighted_norm_square = 1.;
686 last_f_weighted_norm_square = 1;
687 try_project_to_unit_cell = true;
688 continue;
689 }
690 else if (step_length <= 0.05)
691 return invalid_point;
692
693 ++newton_iteration;
694 if (newton_iteration > newton_iteration_limit)
695 return invalid_point;
696 }
697 // Stop if f_weighted_norm_square <= eps^2 on all SIMD lanes or if the
698 // weighted norm is less than 1e-6 and the improvement against the
699 // previous step was less than a factor of 10 (in that regime, we
700 // either have quadratic convergence or roundoff errors due to a bad
701 // mapping)
702 while (
703 !(std::max(f_weighted_norm_square - eps * eps, Number(0.)) *
704 std::max(last_f_weighted_norm_square -
705 std::min(f_weighted_norm_square, Number(1e-6 * 1e-6)) *
706 100.,
707 Number(0.)) ==
708 Number(0.)));
709
710 if (print_iterations_to_deallog)
711 deallog << "Iteration converged for p_unit = [ " << p_unit
712 << " ] and iteration error "
713 << std::sqrt(f_weighted_norm_square) << std::endl;
714
715 return p_unit;
716 }
717
718
719
723 template <int dim>
724 inline Point<dim>
726 const typename ::Triangulation<dim, dim + 1>::cell_iterator &cell,
727 const Point<dim + 1> & p,
728 const Point<dim> &initial_p_unit,
729 typename ::MappingQGeneric<dim, dim + 1>::InternalData &mdata)
730 {
731 const unsigned int spacedim = dim + 1;
732
733 const unsigned int n_shapes = mdata.shape_values.size();
734 (void)n_shapes;
735 Assert(n_shapes != 0, ExcInternalError());
736 Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
737 Assert(mdata.shape_second_derivatives.size() == n_shapes,
739
740 std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
741 Assert(points.size() == n_shapes, ExcInternalError());
742
743 Point<spacedim> p_minus_F;
744
745 Tensor<1, spacedim> DF[dim];
746 Tensor<1, spacedim> D2F[dim][dim];
747
748 Point<dim> p_unit = initial_p_unit;
749 Point<dim> f;
751
752 // Evaluate first and second derivatives
753 mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
754
755 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
756 {
757 const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
758 const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
759 const Point<spacedim> &point_k = points[k];
760
761 for (unsigned int j = 0; j < dim; ++j)
762 {
763 DF[j] += grad_phi_k[j] * point_k;
764 for (unsigned int l = 0; l < dim; ++l)
765 D2F[j][l] += hessian_k[j][l] * point_k;
766 }
767 }
768
769 p_minus_F = p;
770 p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
771
772
773 for (unsigned int j = 0; j < dim; ++j)
774 f[j] = DF[j] * p_minus_F;
775
776 for (unsigned int j = 0; j < dim; ++j)
777 {
778 f[j] = DF[j] * p_minus_F;
779 for (unsigned int l = 0; l < dim; ++l)
780 df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
781 }
782
783
784 const double eps = 1.e-12 * cell->diameter();
785 const unsigned int loop_limit = 10;
786
787 unsigned int loop = 0;
788
789 while (f.norm() > eps && loop++ < loop_limit)
790 {
791 // Solve [df(x)]d=f(x)
792 const Tensor<1, dim> d =
793 invert(df) * static_cast<const Tensor<1, dim> &>(f);
794 p_unit -= d;
795
796 for (unsigned int j = 0; j < dim; ++j)
797 {
798 DF[j].clear();
799 for (unsigned int l = 0; l < dim; ++l)
800 D2F[j][l].clear();
801 }
802
803 mdata.compute_shape_function_values(
804 std::vector<Point<dim>>(1, p_unit));
805
806 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
807 {
808 const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
809 const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
810 const Point<spacedim> &point_k = points[k];
811
812 for (unsigned int j = 0; j < dim; ++j)
813 {
814 DF[j] += grad_phi_k[j] * point_k;
815 for (unsigned int l = 0; l < dim; ++l)
816 D2F[j][l] += hessian_k[j][l] * point_k;
817 }
818 }
819
820 // TODO: implement a line search here in much the same way as for
821 // the corresponding function above that does so for dim==spacedim
822 p_minus_F = p;
823 p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
824
825 for (unsigned int j = 0; j < dim; ++j)
826 {
827 f[j] = DF[j] * p_minus_F;
828 for (unsigned int l = 0; l < dim; ++l)
829 df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
830 }
831 }
832
833
834 // Here we check that in the last execution of while the first
835 // condition was already wrong, meaning the residual was below
836 // eps. Only if the first condition failed, loop will have been
837 // increased and tested, and thus have reached the limit.
838 AssertThrow(loop < loop_limit,
840
841 return p_unit;
842 }
843
844
845
863 template <int dim, int spacedim>
865 {
866 public:
870 static constexpr unsigned int n_functions =
871 (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
872
885 const std::vector<Point<spacedim>> &real_support_points,
886 const std::vector<Point<dim>> & unit_support_points)
887 : normalization_shift(real_support_points[0])
889 1. / real_support_points[0].distance(real_support_points[1]))
890 , is_affine(true)
891 {
892 AssertDimension(real_support_points.size(), unit_support_points.size());
893
894 // For the bi-/trilinear approximation, we cannot build a quadratic
895 // polynomial due to a lack of points (interpolation matrix would get
896 // singular), so pick the affine approximation. Similarly, it is not
897 // entirely clear how to gather enough information for the case dim <
898 // spacedim
899 if (real_support_points.size() ==
901 dim < spacedim)
902 {
903 const auto affine = GridTools::affine_cell_approximation<dim>(
904 make_array_view(real_support_points));
906 affine.first.covariant_form().transpose();
907 coefficients[0] = apply_transformation(A_inv, affine.second);
908 for (unsigned int d = 0; d < spacedim; ++d)
909 for (unsigned int e = 0; e < dim; ++e)
910 coefficients[1 + d][e] = A_inv[e][d];
911 is_affine = true;
912 return;
913 }
914
916 std::array<double, n_functions> shape_values;
917 for (unsigned int q = 0; q < unit_support_points.size(); ++q)
918 {
919 // Evaluate quadratic shape functions in point, with the
920 // normalization applied in order to avoid roundoff issues with
921 // scaling far away from 1.
922 shape_values[0] = 1.;
923 const Tensor<1, spacedim> p_scaled =
925 (real_support_points[q] - normalization_shift);
926 for (unsigned int d = 0; d < spacedim; ++d)
927 shape_values[1 + d] = p_scaled[d];
928 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
929 for (unsigned int e = 0; e <= d; ++e, ++c)
930 shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
931
932 // Build lower diagonal of least squares matrix and rhs, the
933 // essential part being that we construct the matrix with the
934 // real points and the right hand side by comparing to the
935 // reference point positions which sets up an inverse
936 // interpolation.
937 for (unsigned int i = 0; i < n_functions; ++i)
938 for (unsigned int j = 0; j < n_functions; ++j)
939 matrix[i][j] += shape_values[i] * shape_values[j];
940 for (unsigned int i = 0; i < n_functions; ++i)
941 coefficients[i] += shape_values[i] * unit_support_points[q];
942 }
943
944 // Factorize the matrix A = L * L^T in-place with the
945 // Cholesky-Banachiewicz algorithm. The implementation is similar to
946 // FullMatrix::cholesky() but re-implemented to avoid memory
947 // allocations and some unnecessary divisions which we can do here as
948 // we only need to solve with dim right hand sides.
949 for (unsigned int i = 0; i < n_functions; ++i)
950 {
951 double Lij_sum = 0;
952 for (unsigned int j = 0; j < i; ++j)
953 {
954 double Lik_Ljk_sum = 0;
955 for (unsigned int k = 0; k < j; ++k)
956 Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
957 matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
958 Lij_sum += matrix[i][j] * matrix[i][j];
959 }
960 AssertThrow(matrix[i][i] - Lij_sum >= 0,
961 ExcMessage("Matrix not positive definite"));
962
963 // Store the inverse in the diagonal since that is the quantity
964 // needed later in the factorization as well as the forward and
965 // backward substitution, minimizing the number of divisions.
966 matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
967 }
968
969 // Solve lower triangular part, L * y = rhs.
970 for (unsigned int i = 0; i < n_functions; ++i)
971 {
973 for (unsigned int j = 0; j < i; ++j)
974 sum -= matrix[i][j] * coefficients[j];
975 coefficients[i] = sum * matrix[i][i];
976 }
977
978 // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
979 for (unsigned int i = n_functions; i > 0;)
980 {
981 --i;
983 for (unsigned int j = i + 1; j < n_functions; ++j)
984 sum -= matrix[j][i] * coefficients[j];
985 coefficients[i] = sum * matrix[i][i];
986 }
987
988 // Check whether the approximation is indeed affine, allowing to
989 // skip the quadratic terms.
990 is_affine = true;
991 for (unsigned int i = dim + 1; i < n_functions; ++i)
992 if (coefficients[i].norm_square() > 1e-20)
993 {
994 is_affine = false;
995 break;
996 }
997 }
998
1003 default;
1004
1008 template <typename Number>
1011 {
1012 Point<dim, Number> result;
1013 for (unsigned int d = 0; d < dim; ++d)
1014 result[d] = coefficients[0][d];
1015
1016 // Apply the normalization to ensure a good conditioning. Since Number
1017 // might be a vectorized array whereas the normalization is a point of
1018 // doubles, we cannot use the overload of operator- and must instead
1019 // loop over the components of the point.
1020 Point<spacedim, Number> p_scaled;
1021 for (unsigned int d = 0; d < spacedim; ++d)
1022 p_scaled[d] = (p[d] - normalization_shift[d]) * normalization_length;
1023
1024 for (unsigned int d = 0; d < spacedim; ++d)
1025 result += coefficients[1 + d] * p_scaled[d];
1026
1027 if (!is_affine)
1028 {
1029 for (unsigned int d = 0, c = 0; d < spacedim; ++d)
1030 for (unsigned int e = 0; e <= d; ++e, ++c)
1031 result +=
1032 coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
1033 }
1034 return result;
1035 }
1036
1037 private:
1046
1051
1055 std::array<Point<dim>, n_functions> coefficients;
1056
1063 };
1064
1065
1066
1072 template <int dim, int spacedim>
1073 inline void
1075 const CellSimilarity::Similarity cell_similarity,
1076 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1077 std::vector<Point<spacedim>> & quadrature_points,
1078 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1079 {
1080 const UpdateFlags update_flags = data.update_each;
1081
1082 const unsigned int n_shape_values = data.n_shape_functions;
1083 const unsigned int n_q_points = data.shape_info.n_q_points;
1084 constexpr unsigned int n_lanes = VectorizedArray<double>::size();
1085 constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1086 constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1087
1088 EvaluationFlags::EvaluationFlags evaluation_flag =
1091 ((cell_similarity != CellSimilarity::translation) &&
1092 (update_flags & update_contravariant_transformation) ?
1095 ((cell_similarity != CellSimilarity::translation) &&
1096 (update_flags & update_jacobian_grads) ?
1099
1100 Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1102 Assert(!(evaluation_flag & EvaluationFlags::values) ||
1103 n_q_points == quadrature_points.size(),
1104 ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1105 Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1106 data.n_shape_functions > 0,
1108 Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1109 n_q_points == data.contravariant.size(),
1110 ExcDimensionMismatch(n_q_points, data.contravariant.size()));
1111 Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
1112 n_q_points == jacobian_grads.size(),
1113 ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1114
1115 // shortcut in case we have an identity interpolation and only request
1116 // the quadrature points
1117 if (evaluation_flag == EvaluationFlags::values &&
1118 data.shape_info.element_type ==
1120 {
1121 for (unsigned int q = 0; q < n_q_points; ++q)
1123 data.mapping_support_points[data.shape_info
1124 .lexicographic_numbering[q]];
1125 return;
1126 }
1127
1128 // prepare arrays
1129 if (evaluation_flag != EvaluationFlags::nothing)
1130 {
1131 data.values_dofs.resize(n_comp * n_shape_values);
1132 data.values_quad.resize(n_comp * n_q_points);
1133 data.gradients_quad.resize(n_comp * n_q_points * dim);
1134 data.scratch.resize(2 * std::max(n_q_points, n_shape_values));
1135
1136 if (evaluation_flag & EvaluationFlags::hessians)
1137 data.hessians_quad.resize(n_comp * n_q_points * n_hessians);
1138
1139 const std::vector<unsigned int> &renumber_to_lexicographic =
1140 data.shape_info.lexicographic_numbering;
1141 for (unsigned int i = 0; i < n_shape_values; ++i)
1142 for (unsigned int d = 0; d < spacedim; ++d)
1143 {
1144 const unsigned int in_comp = d % n_lanes;
1145 const unsigned int out_comp = d / n_lanes;
1146 data.values_dofs[out_comp * n_shape_values + i][in_comp] =
1147 data.mapping_support_points[renumber_to_lexicographic[i]][d];
1148 }
1149
1150 // do the actual tensorized evaluation
1152 evaluate(n_comp,
1153 evaluation_flag,
1154 data.shape_info,
1155 data.values_dofs.begin(),
1156 data.values_quad.begin(),
1157 data.gradients_quad.begin(),
1158 data.hessians_quad.begin(),
1159 data.scratch.begin());
1160 }
1161
1162 // do the postprocessing
1163 if (evaluation_flag & EvaluationFlags::values)
1164 {
1165 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1166 for (unsigned int i = 0; i < n_q_points; ++i)
1167 for (unsigned int in_comp = 0;
1168 in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1169 ++in_comp)
1170 quadrature_points[i][out_comp * n_lanes + in_comp] =
1171 data.values_quad[out_comp * n_q_points + i][in_comp];
1172 }
1173
1174 if (evaluation_flag & EvaluationFlags::gradients)
1175 {
1176 std::fill(data.contravariant.begin(),
1177 data.contravariant.end(),
1179 // We need to reinterpret the data after evaluate has been applied.
1180 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1181 for (unsigned int point = 0; point < n_q_points; ++point)
1182 for (unsigned int j = 0; j < dim; ++j)
1183 for (unsigned int in_comp = 0;
1184 in_comp < n_lanes &&
1185 in_comp < spacedim - out_comp * n_lanes;
1186 ++in_comp)
1187 {
1188 const unsigned int total_number = point * dim + j;
1189 const unsigned int new_comp = total_number / n_q_points;
1190 const unsigned int new_point = total_number % n_q_points;
1191 data.contravariant[new_point][out_comp * n_lanes + in_comp]
1192 [new_comp] =
1193 data
1194 .gradients_quad[(out_comp * n_q_points + point) * dim +
1195 j][in_comp];
1196 }
1197 }
1198 if (update_flags & update_covariant_transformation)
1199 if (cell_similarity != CellSimilarity::translation)
1200 for (unsigned int point = 0; point < n_q_points; ++point)
1201 data.covariant[point] =
1202 (data.contravariant[point]).covariant_form();
1203
1204 if (update_flags & update_volume_elements)
1205 if (cell_similarity != CellSimilarity::translation)
1206 for (unsigned int point = 0; point < n_q_points; ++point)
1207 data.volume_elements[point] =
1208 data.contravariant[point].determinant();
1209
1210 if (evaluation_flag & EvaluationFlags::hessians)
1211 {
1212 constexpr int desymmetrize_3d[6][2] = {
1213 {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1214 constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1215
1216 // We need to reinterpret the data after evaluate has been applied.
1217 for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1218 for (unsigned int point = 0; point < n_q_points; ++point)
1219 for (unsigned int j = 0; j < n_hessians; ++j)
1220 for (unsigned int in_comp = 0;
1221 in_comp < n_lanes &&
1222 in_comp < spacedim - out_comp * n_lanes;
1223 ++in_comp)
1224 {
1225 const unsigned int total_number = point * n_hessians + j;
1226 const unsigned int new_point = total_number % n_q_points;
1227 const unsigned int new_hessian_comp =
1228 total_number / n_q_points;
1229 const unsigned int new_hessian_comp_i =
1230 dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1231 desymmetrize_3d[new_hessian_comp][0];
1232 const unsigned int new_hessian_comp_j =
1233 dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1234 desymmetrize_3d[new_hessian_comp][1];
1235 const double value =
1236 data.hessians_quad[(out_comp * n_q_points + point) *
1237 n_hessians +
1238 j][in_comp];
1239 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1240 [new_hessian_comp_i][new_hessian_comp_j] =
1241 value;
1242 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1243 [new_hessian_comp_j][new_hessian_comp_i] =
1244 value;
1245 }
1246 }
1247 }
1248
1249
1256 template <int dim, int spacedim>
1257 inline void
1259 const typename QProjector<dim>::DataSetDescriptor data_set,
1260 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1261 std::vector<Point<spacedim>> &quadrature_points)
1262 {
1263 const UpdateFlags update_flags = data.update_each;
1264
1265 if (update_flags & update_quadrature_points)
1266 for (unsigned int point = 0; point < quadrature_points.size(); ++point)
1267 {
1268 const double * shape = &data.shape(point + data_set, 0);
1269 Point<spacedim> result =
1270 (shape[0] * data.mapping_support_points[0]);
1271 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1272 for (unsigned int i = 0; i < spacedim; ++i)
1273 result[i] += shape[k] * data.mapping_support_points[k][i];
1274 quadrature_points[point] = result;
1275 }
1276 }
1277
1278
1279
1288 template <int dim, int spacedim>
1289 inline void
1291 const CellSimilarity::Similarity cell_similarity,
1292 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1293 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data)
1294 {
1295 const UpdateFlags update_flags = data.update_each;
1296
1297 if (update_flags & update_contravariant_transformation)
1298 // if the current cell is just a
1299 // translation of the previous one, no
1300 // need to recompute jacobians...
1301 if (cell_similarity != CellSimilarity::translation)
1302 {
1303 const unsigned int n_q_points = data.contravariant.size();
1304
1305 std::fill(data.contravariant.begin(),
1306 data.contravariant.end(),
1308
1309 Assert(data.n_shape_functions > 0, ExcInternalError());
1310
1311 for (unsigned int point = 0; point < n_q_points; ++point)
1312 {
1313 double result[spacedim][dim];
1314
1315 // peel away part of sum to avoid zeroing the
1316 // entries and adding for the first time
1317 for (unsigned int i = 0; i < spacedim; ++i)
1318 for (unsigned int j = 0; j < dim; ++j)
1319 result[i][j] = data.derivative(point + data_set, 0)[j] *
1320 data.mapping_support_points[0][i];
1321 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1322 for (unsigned int i = 0; i < spacedim; ++i)
1323 for (unsigned int j = 0; j < dim; ++j)
1324 result[i][j] += data.derivative(point + data_set, k)[j] *
1325 data.mapping_support_points[k][i];
1326
1327 // write result into contravariant data. for
1328 // j=dim in the case dim<spacedim, there will
1329 // never be any nonzero data that arrives in
1330 // here, so it is ok anyway because it was
1331 // initialized to zero at the initialization
1332 for (unsigned int i = 0; i < spacedim; ++i)
1333 for (unsigned int j = 0; j < dim; ++j)
1334 data.contravariant[point][i][j] = result[i][j];
1335 }
1336 }
1337
1338 if (update_flags & update_covariant_transformation)
1339 if (cell_similarity != CellSimilarity::translation)
1340 {
1341 const unsigned int n_q_points = data.contravariant.size();
1342 for (unsigned int point = 0; point < n_q_points; ++point)
1343 {
1344 data.covariant[point] =
1345 (data.contravariant[point]).covariant_form();
1346 }
1347 }
1348
1349 if (update_flags & update_volume_elements)
1350 if (cell_similarity != CellSimilarity::translation)
1351 {
1352 const unsigned int n_q_points = data.contravariant.size();
1353 for (unsigned int point = 0; point < n_q_points; ++point)
1354 data.volume_elements[point] =
1355 data.contravariant[point].determinant();
1356 }
1357 }
1358
1359
1360
1367 template <int dim, int spacedim>
1368 inline void
1370 const CellSimilarity::Similarity cell_similarity,
1371 const typename QProjector<dim>::DataSetDescriptor data_set,
1372 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1373 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1374 {
1375 const UpdateFlags update_flags = data.update_each;
1376 if (update_flags & update_jacobian_grads)
1377 {
1378 const unsigned int n_q_points = jacobian_grads.size();
1379
1380 if (cell_similarity != CellSimilarity::translation)
1381 for (unsigned int point = 0; point < n_q_points; ++point)
1382 {
1383 const Tensor<2, dim> *second =
1384 &data.second_derivative(point + data_set, 0);
1385 double result[spacedim][dim][dim];
1386 for (unsigned int i = 0; i < spacedim; ++i)
1387 for (unsigned int j = 0; j < dim; ++j)
1388 for (unsigned int l = 0; l < dim; ++l)
1389 result[i][j][l] =
1390 (second[0][j][l] * data.mapping_support_points[0][i]);
1391 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1392 for (unsigned int i = 0; i < spacedim; ++i)
1393 for (unsigned int j = 0; j < dim; ++j)
1394 for (unsigned int l = 0; l < dim; ++l)
1395 result[i][j][l] +=
1396 (second[k][j][l] * data.mapping_support_points[k][i]);
1397
1398 for (unsigned int i = 0; i < spacedim; ++i)
1399 for (unsigned int j = 0; j < dim; ++j)
1400 for (unsigned int l = 0; l < dim; ++l)
1401 jacobian_grads[point][i][j][l] = result[i][j][l];
1402 }
1403 }
1404 }
1405
1406
1407
1414 template <int dim, int spacedim>
1415 inline void
1417 const CellSimilarity::Similarity cell_similarity,
1418 const typename QProjector<dim>::DataSetDescriptor data_set,
1419 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1420 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1421 {
1422 const UpdateFlags update_flags = data.update_each;
1423 if (update_flags & update_jacobian_pushed_forward_grads)
1424 {
1425 const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
1426
1427 if (cell_similarity != CellSimilarity::translation)
1428 {
1429 double tmp[spacedim][spacedim][spacedim];
1430 for (unsigned int point = 0; point < n_q_points; ++point)
1431 {
1432 const Tensor<2, dim> *second =
1433 &data.second_derivative(point + data_set, 0);
1434 double result[spacedim][dim][dim];
1435 for (unsigned int i = 0; i < spacedim; ++i)
1436 for (unsigned int j = 0; j < dim; ++j)
1437 for (unsigned int l = 0; l < dim; ++l)
1438 result[i][j][l] =
1439 (second[0][j][l] * data.mapping_support_points[0][i]);
1440 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1441 for (unsigned int i = 0; i < spacedim; ++i)
1442 for (unsigned int j = 0; j < dim; ++j)
1443 for (unsigned int l = 0; l < dim; ++l)
1444 result[i][j][l] +=
1445 (second[k][j][l] *
1446 data.mapping_support_points[k][i]);
1447
1448 // first push forward the j-components
1449 for (unsigned int i = 0; i < spacedim; ++i)
1450 for (unsigned int j = 0; j < spacedim; ++j)
1451 for (unsigned int l = 0; l < dim; ++l)
1452 {
1453 tmp[i][j][l] =
1454 result[i][0][l] * data.covariant[point][j][0];
1455 for (unsigned int jr = 1; jr < dim; ++jr)
1456 {
1457 tmp[i][j][l] +=
1458 result[i][jr][l] * data.covariant[point][j][jr];
1459 }
1460 }
1461
1462 // now, pushing forward the l-components
1463 for (unsigned int i = 0; i < spacedim; ++i)
1464 for (unsigned int j = 0; j < spacedim; ++j)
1465 for (unsigned int l = 0; l < spacedim; ++l)
1466 {
1467 jacobian_pushed_forward_grads[point][i][j][l] =
1468 tmp[i][j][0] * data.covariant[point][l][0];
1469 for (unsigned int lr = 1; lr < dim; ++lr)
1470 {
1471 jacobian_pushed_forward_grads[point][i][j][l] +=
1472 tmp[i][j][lr] * data.covariant[point][l][lr];
1473 }
1474 }
1475 }
1476 }
1477 }
1478 }
1479
1480
1481
1488 template <int dim, int spacedim>
1489 inline void
1491 const CellSimilarity::Similarity cell_similarity,
1492 const typename QProjector<dim>::DataSetDescriptor data_set,
1493 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1494 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1495 {
1496 const UpdateFlags update_flags = data.update_each;
1497 if (update_flags & update_jacobian_2nd_derivatives)
1498 {
1499 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1500
1501 if (cell_similarity != CellSimilarity::translation)
1502 {
1503 for (unsigned int point = 0; point < n_q_points; ++point)
1504 {
1505 const Tensor<3, dim> *third =
1506 &data.third_derivative(point + data_set, 0);
1507 double result[spacedim][dim][dim][dim];
1508 for (unsigned int i = 0; i < spacedim; ++i)
1509 for (unsigned int j = 0; j < dim; ++j)
1510 for (unsigned int l = 0; l < dim; ++l)
1511 for (unsigned int m = 0; m < dim; ++m)
1512 result[i][j][l][m] =
1513 (third[0][j][l][m] *
1514 data.mapping_support_points[0][i]);
1515 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1516 for (unsigned int i = 0; i < spacedim; ++i)
1517 for (unsigned int j = 0; j < dim; ++j)
1518 for (unsigned int l = 0; l < dim; ++l)
1519 for (unsigned int m = 0; m < dim; ++m)
1520 result[i][j][l][m] +=
1521 (third[k][j][l][m] *
1522 data.mapping_support_points[k][i]);
1523
1524 for (unsigned int i = 0; i < spacedim; ++i)
1525 for (unsigned int j = 0; j < dim; ++j)
1526 for (unsigned int l = 0; l < dim; ++l)
1527 for (unsigned int m = 0; m < dim; ++m)
1528 jacobian_2nd_derivatives[point][i][j][l][m] =
1529 result[i][j][l][m];
1530 }
1531 }
1532 }
1533 }
1534
1535
1536
1544 template <int dim, int spacedim>
1545 inline void
1547 const CellSimilarity::Similarity cell_similarity,
1548 const typename QProjector<dim>::DataSetDescriptor data_set,
1549 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1550 std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
1551 {
1552 const UpdateFlags update_flags = data.update_each;
1554 {
1555 const unsigned int n_q_points =
1556 jacobian_pushed_forward_2nd_derivatives.size();
1557
1558 if (cell_similarity != CellSimilarity::translation)
1559 {
1560 double tmp[spacedim][spacedim][spacedim][spacedim];
1561 for (unsigned int point = 0; point < n_q_points; ++point)
1562 {
1563 const Tensor<3, dim> *third =
1564 &data.third_derivative(point + data_set, 0);
1565 double result[spacedim][dim][dim][dim];
1566 for (unsigned int i = 0; i < spacedim; ++i)
1567 for (unsigned int j = 0; j < dim; ++j)
1568 for (unsigned int l = 0; l < dim; ++l)
1569 for (unsigned int m = 0; m < dim; ++m)
1570 result[i][j][l][m] =
1571 (third[0][j][l][m] *
1572 data.mapping_support_points[0][i]);
1573 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1574 for (unsigned int i = 0; i < spacedim; ++i)
1575 for (unsigned int j = 0; j < dim; ++j)
1576 for (unsigned int l = 0; l < dim; ++l)
1577 for (unsigned int m = 0; m < dim; ++m)
1578 result[i][j][l][m] +=
1579 (third[k][j][l][m] *
1580 data.mapping_support_points[k][i]);
1581
1582 // push forward the j-coordinate
1583 for (unsigned int i = 0; i < spacedim; ++i)
1584 for (unsigned int j = 0; j < spacedim; ++j)
1585 for (unsigned int l = 0; l < dim; ++l)
1586 for (unsigned int m = 0; m < dim; ++m)
1587 {
1588 jacobian_pushed_forward_2nd_derivatives
1589 [point][i][j][l][m] = result[i][0][l][m] *
1590 data.covariant[point][j][0];
1591 for (unsigned int jr = 1; jr < dim; ++jr)
1592 jacobian_pushed_forward_2nd_derivatives[point][i]
1593 [j][l]
1594 [m] +=
1595 result[i][jr][l][m] *
1596 data.covariant[point][j][jr];
1597 }
1598
1599 // push forward the l-coordinate
1600 for (unsigned int i = 0; i < spacedim; ++i)
1601 for (unsigned int j = 0; j < spacedim; ++j)
1602 for (unsigned int l = 0; l < spacedim; ++l)
1603 for (unsigned int m = 0; m < dim; ++m)
1604 {
1605 tmp[i][j][l][m] =
1606 jacobian_pushed_forward_2nd_derivatives[point][i]
1607 [j][0][m] *
1608 data.covariant[point][l][0];
1609 for (unsigned int lr = 1; lr < dim; ++lr)
1610 tmp[i][j][l][m] +=
1611 jacobian_pushed_forward_2nd_derivatives[point]
1612 [i][j]
1613 [lr][m] *
1614 data.covariant[point][l][lr];
1615 }
1616
1617 // push forward the m-coordinate
1618 for (unsigned int i = 0; i < spacedim; ++i)
1619 for (unsigned int j = 0; j < spacedim; ++j)
1620 for (unsigned int l = 0; l < spacedim; ++l)
1621 for (unsigned int m = 0; m < spacedim; ++m)
1622 {
1623 jacobian_pushed_forward_2nd_derivatives
1624 [point][i][j][l][m] =
1625 tmp[i][j][l][0] * data.covariant[point][m][0];
1626 for (unsigned int mr = 1; mr < dim; ++mr)
1627 jacobian_pushed_forward_2nd_derivatives[point][i]
1628 [j][l]
1629 [m] +=
1630 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1631 }
1632 }
1633 }
1634 }
1635 }
1636
1637
1638
1645 template <int dim, int spacedim>
1646 inline void
1648 const CellSimilarity::Similarity cell_similarity,
1649 const typename QProjector<dim>::DataSetDescriptor data_set,
1650 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1651 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1652 {
1653 const UpdateFlags update_flags = data.update_each;
1654 if (update_flags & update_jacobian_3rd_derivatives)
1655 {
1656 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1657
1658 if (cell_similarity != CellSimilarity::translation)
1659 {
1660 for (unsigned int point = 0; point < n_q_points; ++point)
1661 {
1662 const Tensor<4, dim> *fourth =
1663 &data.fourth_derivative(point + data_set, 0);
1664 double result[spacedim][dim][dim][dim][dim];
1665 for (unsigned int i = 0; i < spacedim; ++i)
1666 for (unsigned int j = 0; j < dim; ++j)
1667 for (unsigned int l = 0; l < dim; ++l)
1668 for (unsigned int m = 0; m < dim; ++m)
1669 for (unsigned int n = 0; n < dim; ++n)
1670 result[i][j][l][m][n] =
1671 (fourth[0][j][l][m][n] *
1672 data.mapping_support_points[0][i]);
1673 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1674 for (unsigned int i = 0; i < spacedim; ++i)
1675 for (unsigned int j = 0; j < dim; ++j)
1676 for (unsigned int l = 0; l < dim; ++l)
1677 for (unsigned int m = 0; m < dim; ++m)
1678 for (unsigned int n = 0; n < dim; ++n)
1679 result[i][j][l][m][n] +=
1680 (fourth[k][j][l][m][n] *
1681 data.mapping_support_points[k][i]);
1682
1683 for (unsigned int i = 0; i < spacedim; ++i)
1684 for (unsigned int j = 0; j < dim; ++j)
1685 for (unsigned int l = 0; l < dim; ++l)
1686 for (unsigned int m = 0; m < dim; ++m)
1687 for (unsigned int n = 0; n < dim; ++n)
1688 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1689 result[i][j][l][m][n];
1690 }
1691 }
1692 }
1693 }
1694
1695
1696
1704 template <int dim, int spacedim>
1705 inline void
1707 const CellSimilarity::Similarity cell_similarity,
1708 const typename QProjector<dim>::DataSetDescriptor data_set,
1709 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1710 std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
1711 {
1712 const UpdateFlags update_flags = data.update_each;
1714 {
1715 const unsigned int n_q_points =
1716 jacobian_pushed_forward_3rd_derivatives.size();
1717
1718 if (cell_similarity != CellSimilarity::translation)
1719 {
1720 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1721 for (unsigned int point = 0; point < n_q_points; ++point)
1722 {
1723 const Tensor<4, dim> *fourth =
1724 &data.fourth_derivative(point + data_set, 0);
1725 double result[spacedim][dim][dim][dim][dim];
1726 for (unsigned int i = 0; i < spacedim; ++i)
1727 for (unsigned int j = 0; j < dim; ++j)
1728 for (unsigned int l = 0; l < dim; ++l)
1729 for (unsigned int m = 0; m < dim; ++m)
1730 for (unsigned int n = 0; n < dim; ++n)
1731 result[i][j][l][m][n] =
1732 (fourth[0][j][l][m][n] *
1733 data.mapping_support_points[0][i]);
1734 for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1735 for (unsigned int i = 0; i < spacedim; ++i)
1736 for (unsigned int j = 0; j < dim; ++j)
1737 for (unsigned int l = 0; l < dim; ++l)
1738 for (unsigned int m = 0; m < dim; ++m)
1739 for (unsigned int n = 0; n < dim; ++n)
1740 result[i][j][l][m][n] +=
1741 (fourth[k][j][l][m][n] *
1742 data.mapping_support_points[k][i]);
1743
1744 // push-forward the j-coordinate
1745 for (unsigned int i = 0; i < spacedim; ++i)
1746 for (unsigned int j = 0; j < spacedim; ++j)
1747 for (unsigned int l = 0; l < dim; ++l)
1748 for (unsigned int m = 0; m < dim; ++m)
1749 for (unsigned int n = 0; n < dim; ++n)
1750 {
1751 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1752 data.covariant[point][j][0];
1753 for (unsigned int jr = 1; jr < dim; ++jr)
1754 tmp[i][j][l][m][n] +=
1755 result[i][jr][l][m][n] *
1756 data.covariant[point][j][jr];
1757 }
1758
1759 // push-forward the l-coordinate
1760 for (unsigned int i = 0; i < spacedim; ++i)
1761 for (unsigned int j = 0; j < spacedim; ++j)
1762 for (unsigned int l = 0; l < spacedim; ++l)
1763 for (unsigned int m = 0; m < dim; ++m)
1764 for (unsigned int n = 0; n < dim; ++n)
1765 {
1766 jacobian_pushed_forward_3rd_derivatives
1767 [point][i][j][l][m][n] =
1768 tmp[i][j][0][m][n] *
1769 data.covariant[point][l][0];
1770 for (unsigned int lr = 1; lr < dim; ++lr)
1771 jacobian_pushed_forward_3rd_derivatives[point]
1772 [i][j][l]
1773 [m][n] +=
1774 tmp[i][j][lr][m][n] *
1775 data.covariant[point][l][lr];
1776 }
1777
1778 // push-forward the m-coordinate
1779 for (unsigned int i = 0; i < spacedim; ++i)
1780 for (unsigned int j = 0; j < spacedim; ++j)
1781 for (unsigned int l = 0; l < spacedim; ++l)
1782 for (unsigned int m = 0; m < spacedim; ++m)
1783 for (unsigned int n = 0; n < dim; ++n)
1784 {
1785 tmp[i][j][l][m][n] =
1786 jacobian_pushed_forward_3rd_derivatives[point]
1787 [i][j][l]
1788 [0][n] *
1789 data.covariant[point][m][0];
1790 for (unsigned int mr = 1; mr < dim; ++mr)
1791 tmp[i][j][l][m][n] +=
1792 jacobian_pushed_forward_3rd_derivatives
1793 [point][i][j][l][mr][n] *
1794 data.covariant[point][m][mr];
1795 }
1796
1797 // push-forward the n-coordinate
1798 for (unsigned int i = 0; i < spacedim; ++i)
1799 for (unsigned int j = 0; j < spacedim; ++j)
1800 for (unsigned int l = 0; l < spacedim; ++l)
1801 for (unsigned int m = 0; m < spacedim; ++m)
1802 for (unsigned int n = 0; n < spacedim; ++n)
1803 {
1804 jacobian_pushed_forward_3rd_derivatives
1805 [point][i][j][l][m][n] =
1806 tmp[i][j][l][m][0] *
1807 data.covariant[point][n][0];
1808 for (unsigned int nr = 1; nr < dim; ++nr)
1809 jacobian_pushed_forward_3rd_derivatives[point]
1810 [i][j][l]
1811 [m][n] +=
1812 tmp[i][j][l][m][nr] *
1813 data.covariant[point][n][nr];
1814 }
1815 }
1816 }
1817 }
1818 }
1819
1820
1821
1831 template <int dim, int spacedim>
1832 inline void
1834 const ::MappingQGeneric<dim, spacedim> &mapping,
1835 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
1836 const unsigned int face_no,
1837 const unsigned int subface_no,
1838 const unsigned int n_q_points,
1839 const std::vector<double> &weights,
1840 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
1842 &output_data)
1843 {
1844 const UpdateFlags update_flags = data.update_each;
1845
1846 if (update_flags &
1849 {
1850 if (update_flags & update_boundary_forms)
1851 AssertDimension(output_data.boundary_forms.size(), n_q_points);
1852 if (update_flags & update_normal_vectors)
1853 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1854 if (update_flags & update_JxW_values)
1855 AssertDimension(output_data.JxW_values.size(), n_q_points);
1856
1857 Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1858
1859 // first compute some common data that is used for evaluating
1860 // all of the flags below
1861
1862 // map the unit tangentials to the real cell. checking for d!=dim-1
1863 // eliminates compiler warnings regarding unsigned int expressions <
1864 // 0.
1865 for (unsigned int d = 0; d != dim - 1; ++d)
1866 {
1868 data.unit_tangentials.size(),
1870 Assert(
1871 data.aux[d].size() <=
1872 data
1873 .unit_tangentials[face_no +
1875 .size(),
1877
1878 mapping.transform(
1880 data.unit_tangentials[face_no +
1883 data,
1884 make_array_view(data.aux[d]));
1885 }
1886
1887 if (update_flags & update_boundary_forms)
1888 {
1889 // if dim==spacedim, we can use the unit tangentials to compute
1890 // the boundary form by simply taking the cross product
1891 if (dim == spacedim)
1892 {
1893 for (unsigned int i = 0; i < n_q_points; ++i)
1894 switch (dim)
1895 {
1896 case 1:
1897 // in 1d, we don't have access to any of the
1898 // data.aux fields (because it has only dim-1
1899 // components), but we can still compute the
1900 // boundary form by simply looking at the number of
1901 // the face
1902 output_data.boundary_forms[i][0] =
1903 (face_no == 0 ? -1 : +1);
1904 break;
1905 case 2:
1906 output_data.boundary_forms[i] =
1907 cross_product_2d(data.aux[0][i]);
1908 break;
1909 case 3:
1910 output_data.boundary_forms[i] =
1911 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1912 break;
1913 default:
1914 Assert(false, ExcNotImplemented());
1915 }
1916 }
1917 else //(dim < spacedim)
1918 {
1919 // in the codim-one case, the boundary form results from the
1920 // cross product of all the face tangential vectors and the
1921 // cell normal vector
1922 //
1923 // to compute the cell normal, use the same method used in
1924 // fill_fe_values for cells above
1925 AssertDimension(data.contravariant.size(), n_q_points);
1926
1927 for (unsigned int point = 0; point < n_q_points; ++point)
1928 {
1929 if (dim == 1)
1930 {
1931 // J is a tangent vector
1932 output_data.boundary_forms[point] =
1933 data.contravariant[point].transpose()[0];
1934 output_data.boundary_forms[point] /=
1935 (face_no == 0 ? -1. : +1.) *
1936 output_data.boundary_forms[point].norm();
1937 }
1938
1939 if (dim == 2)
1940 {
1942 data.contravariant[point].transpose();
1943
1944 Tensor<1, spacedim> cell_normal =
1945 cross_product_3d(DX_t[0], DX_t[1]);
1946 cell_normal /= cell_normal.norm();
1947
1948 // then compute the face normal from the face
1949 // tangent and the cell normal:
1950 output_data.boundary_forms[point] =
1951 cross_product_3d(data.aux[0][point], cell_normal);
1952 }
1953 }
1954 }
1955 }
1956
1957 if (update_flags & update_JxW_values)
1958 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
1959 {
1960 output_data.JxW_values[i] =
1961 output_data.boundary_forms[i].norm() * weights[i];
1962
1963 if (subface_no != numbers::invalid_unsigned_int)
1964 {
1965 const double area_ratio = GeometryInfo<dim>::subface_ratio(
1966 cell->subface_case(face_no), subface_no);
1967 output_data.JxW_values[i] *= area_ratio;
1968 }
1969 }
1970
1971 if (update_flags & update_normal_vectors)
1972 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
1973 output_data.normal_vectors[i] =
1974 Point<spacedim>(output_data.boundary_forms[i] /
1975 output_data.boundary_forms[i].norm());
1976
1977 if (update_flags & update_jacobians)
1978 for (unsigned int point = 0; point < n_q_points; ++point)
1979 output_data.jacobians[point] = data.contravariant[point];
1980
1981 if (update_flags & update_inverse_jacobians)
1982 for (unsigned int point = 0; point < n_q_points; ++point)
1983 output_data.inverse_jacobians[point] =
1984 data.covariant[point].transpose();
1985 }
1986 }
1987
1988
1995 template <int dim, int spacedim>
1996 inline void
1998 const ::MappingQGeneric<dim, spacedim> &mapping,
1999 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
2000 const unsigned int face_no,
2001 const unsigned int subface_no,
2002 const typename QProjector<dim>::DataSetDescriptor data_set,
2003 const Quadrature<dim - 1> & quadrature,
2004 const typename ::MappingQGeneric<dim, spacedim>::InternalData &data,
2006 &output_data)
2007 {
2008 if (dim > 1 && data.tensor_product_quadrature)
2009 {
2010 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2012 data,
2013 output_data.quadrature_points,
2014 output_data.jacobian_grads);
2015 }
2016 else
2017 {
2018 maybe_compute_q_points<dim, spacedim>(data_set,
2019 data,
2020 output_data.quadrature_points);
2021 maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
2022 data_set,
2023 data);
2024 maybe_update_jacobian_grads<dim, spacedim>(
2025 CellSimilarity::none, data_set, data, output_data.jacobian_grads);
2026 }
2027 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2029 data_set,
2030 data,
2031 output_data.jacobian_pushed_forward_grads);
2032 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2034 data_set,
2035 data,
2036 output_data.jacobian_2nd_derivatives);
2037 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2039 data_set,
2040 data,
2042 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2044 data_set,
2045 data,
2046 output_data.jacobian_3rd_derivatives);
2047 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2049 data_set,
2050 data,
2052
2054 cell,
2055 face_no,
2056 subface_no,
2057 quadrature.size(),
2058 quadrature.get_weights(),
2059 data,
2060 output_data);
2061 }
2062
2063
2064
2068 template <int dim, int spacedim, int rank>
2069 inline void
2071 const ArrayView<const Tensor<rank, dim>> & input,
2072 const MappingKind mapping_kind,
2073 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2074 const ArrayView<Tensor<rank, spacedim>> & output)
2075 {
2076 AssertDimension(input.size(), output.size());
2077 Assert((dynamic_cast<const typename ::
2078 MappingQGeneric<dim, spacedim>::InternalData *>(
2079 &mapping_data) != nullptr),
2081 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2082 &data =
2083 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
2084 InternalData &>(mapping_data);
2085
2086 switch (mapping_kind)
2087 {
2089 {
2092 "update_contravariant_transformation"));
2093
2094 for (unsigned int i = 0; i < output.size(); ++i)
2095 output[i] =
2096 apply_transformation(data.contravariant[i], input[i]);
2097
2098 return;
2099 }
2100
2101 case mapping_piola:
2102 {
2105 "update_contravariant_transformation"));
2106 Assert(data.update_each & update_volume_elements,
2108 "update_volume_elements"));
2109 Assert(rank == 1, ExcMessage("Only for rank 1"));
2110 if (rank != 1)
2111 return;
2112
2113 for (unsigned int i = 0; i < output.size(); ++i)
2114 {
2115 output[i] =
2116 apply_transformation(data.contravariant[i], input[i]);
2117 output[i] /= data.volume_elements[i];
2118 }
2119 return;
2120 }
2121 // We still allow this operation as in the
2122 // reference cell Derivatives are Tensor
2123 // rather than DerivativeForm
2124 case mapping_covariant:
2125 {
2128 "update_covariant_transformation"));
2129
2130 for (unsigned int i = 0; i < output.size(); ++i)
2131 output[i] = apply_transformation(data.covariant[i], input[i]);
2132
2133 return;
2134 }
2135
2136 default:
2137 Assert(false, ExcNotImplemented());
2138 }
2139 }
2140
2141
2142
2146 template <int dim, int spacedim, int rank>
2147 inline void
2149 const ArrayView<const Tensor<rank, dim>> & input,
2150 const MappingKind mapping_kind,
2151 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2152 const ArrayView<Tensor<rank, spacedim>> & output)
2153 {
2154 AssertDimension(input.size(), output.size());
2155 Assert((dynamic_cast<const typename ::
2156 MappingQGeneric<dim, spacedim>::InternalData *>(
2157 &mapping_data) != nullptr),
2159 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2160 &data =
2161 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
2162 InternalData &>(mapping_data);
2163
2164 switch (mapping_kind)
2165 {
2167 {
2168 Assert(data.update_each & update_covariant_transformation,
2170 "update_covariant_transformation"));
2173 "update_contravariant_transformation"));
2174 Assert(rank == 2, ExcMessage("Only for rank 2"));
2175
2176 for (unsigned int i = 0; i < output.size(); ++i)
2177 {
2179 apply_transformation(data.contravariant[i],
2180 transpose(input[i]));
2181 output[i] =
2182 apply_transformation(data.covariant[i], A.transpose());
2183 }
2184
2185 return;
2186 }
2187
2189 {
2190 Assert(data.update_each & update_covariant_transformation,
2192 "update_covariant_transformation"));
2193 Assert(rank == 2, ExcMessage("Only for rank 2"));
2194
2195 for (unsigned int i = 0; i < output.size(); ++i)
2196 {
2198 apply_transformation(data.covariant[i],
2199 transpose(input[i]));
2200 output[i] =
2201 apply_transformation(data.covariant[i], A.transpose());
2202 }
2203
2204 return;
2205 }
2206
2208 {
2209 Assert(data.update_each & update_covariant_transformation,
2211 "update_covariant_transformation"));
2214 "update_contravariant_transformation"));
2215 Assert(data.update_each & update_volume_elements,
2217 "update_volume_elements"));
2218 Assert(rank == 2, ExcMessage("Only for rank 2"));
2219
2220 for (unsigned int i = 0; i < output.size(); ++i)
2221 {
2223 apply_transformation(data.covariant[i], input[i]);
2224 const Tensor<2, spacedim> T =
2225 apply_transformation(data.contravariant[i], A.transpose());
2226
2227 output[i] = transpose(T);
2228 output[i] /= data.volume_elements[i];
2229 }
2230
2231 return;
2232 }
2233
2234 default:
2235 Assert(false, ExcNotImplemented());
2236 }
2237 }
2238
2239
2240
2244 template <int dim, int spacedim>
2245 inline void
2247 const ArrayView<const Tensor<3, dim>> & input,
2248 const MappingKind mapping_kind,
2249 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2250 const ArrayView<Tensor<3, spacedim>> & output)
2251 {
2252 AssertDimension(input.size(), output.size());
2253 Assert((dynamic_cast<const typename ::
2254 MappingQGeneric<dim, spacedim>::InternalData *>(
2255 &mapping_data) != nullptr),
2257 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2258 &data =
2259 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
2260 InternalData &>(mapping_data);
2261
2262 switch (mapping_kind)
2263 {
2265 {
2266 Assert(data.update_each & update_covariant_transformation,
2268 "update_covariant_transformation"));
2271 "update_contravariant_transformation"));
2272
2273 for (unsigned int q = 0; q < output.size(); ++q)
2274 for (unsigned int i = 0; i < spacedim; ++i)
2275 {
2276 double tmp1[dim][dim];
2277 for (unsigned int J = 0; J < dim; ++J)
2278 for (unsigned int K = 0; K < dim; ++K)
2279 {
2280 tmp1[J][K] =
2281 data.contravariant[q][i][0] * input[q][0][J][K];
2282 for (unsigned int I = 1; I < dim; ++I)
2283 tmp1[J][K] +=
2284 data.contravariant[q][i][I] * input[q][I][J][K];
2285 }
2286 for (unsigned int j = 0; j < spacedim; ++j)
2287 {
2288 double tmp2[dim];
2289 for (unsigned int K = 0; K < dim; ++K)
2290 {
2291 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2292 for (unsigned int J = 1; J < dim; ++J)
2293 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2294 }
2295 for (unsigned int k = 0; k < spacedim; ++k)
2296 {
2297 output[q][i][j][k] =
2298 data.covariant[q][k][0] * tmp2[0];
2299 for (unsigned int K = 1; K < dim; ++K)
2300 output[q][i][j][k] +=
2301 data.covariant[q][k][K] * tmp2[K];
2302 }
2303 }
2304 }
2305 return;
2306 }
2307
2309 {
2310 Assert(data.update_each & update_covariant_transformation,
2312 "update_covariant_transformation"));
2313
2314 for (unsigned int q = 0; q < output.size(); ++q)
2315 for (unsigned int i = 0; i < spacedim; ++i)
2316 {
2317 double tmp1[dim][dim];
2318 for (unsigned int J = 0; J < dim; ++J)
2319 for (unsigned int K = 0; K < dim; ++K)
2320 {
2321 tmp1[J][K] =
2322 data.covariant[q][i][0] * input[q][0][J][K];
2323 for (unsigned int I = 1; I < dim; ++I)
2324 tmp1[J][K] +=
2325 data.covariant[q][i][I] * input[q][I][J][K];
2326 }
2327 for (unsigned int j = 0; j < spacedim; ++j)
2328 {
2329 double tmp2[dim];
2330 for (unsigned int K = 0; K < dim; ++K)
2331 {
2332 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2333 for (unsigned int J = 1; J < dim; ++J)
2334 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2335 }
2336 for (unsigned int k = 0; k < spacedim; ++k)
2337 {
2338 output[q][i][j][k] =
2339 data.covariant[q][k][0] * tmp2[0];
2340 for (unsigned int K = 1; K < dim; ++K)
2341 output[q][i][j][k] +=
2342 data.covariant[q][k][K] * tmp2[K];
2343 }
2344 }
2345 }
2346
2347 return;
2348 }
2349
2351 {
2352 Assert(data.update_each & update_covariant_transformation,
2354 "update_covariant_transformation"));
2357 "update_contravariant_transformation"));
2358 Assert(data.update_each & update_volume_elements,
2360 "update_volume_elements"));
2361
2362 for (unsigned int q = 0; q < output.size(); ++q)
2363 for (unsigned int i = 0; i < spacedim; ++i)
2364 {
2365 double factor[dim];
2366 for (unsigned int I = 0; I < dim; ++I)
2367 factor[I] =
2368 data.contravariant[q][i][I] / data.volume_elements[q];
2369 double tmp1[dim][dim];
2370 for (unsigned int J = 0; J < dim; ++J)
2371 for (unsigned int K = 0; K < dim; ++K)
2372 {
2373 tmp1[J][K] = factor[0] * input[q][0][J][K];
2374 for (unsigned int I = 1; I < dim; ++I)
2375 tmp1[J][K] += factor[I] * input[q][I][J][K];
2376 }
2377 for (unsigned int j = 0; j < spacedim; ++j)
2378 {
2379 double tmp2[dim];
2380 for (unsigned int K = 0; K < dim; ++K)
2381 {
2382 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2383 for (unsigned int J = 1; J < dim; ++J)
2384 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2385 }
2386 for (unsigned int k = 0; k < spacedim; ++k)
2387 {
2388 output[q][i][j][k] =
2389 data.covariant[q][k][0] * tmp2[0];
2390 for (unsigned int K = 1; K < dim; ++K)
2391 output[q][i][j][k] +=
2392 data.covariant[q][k][K] * tmp2[K];
2393 }
2394 }
2395 }
2396
2397 return;
2398 }
2399
2400 default:
2401 Assert(false, ExcNotImplemented());
2402 }
2403 }
2404
2405
2406
2411 template <int dim, int spacedim, int rank>
2412 inline void
2415 const MappingKind mapping_kind,
2416 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_data,
2417 const ArrayView<Tensor<rank + 1, spacedim>> & output)
2418 {
2419 AssertDimension(input.size(), output.size());
2420 Assert((dynamic_cast<const typename ::
2421 MappingQGeneric<dim, spacedim>::InternalData *>(
2422 &mapping_data) != nullptr),
2424 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2425 &data =
2426 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
2427 InternalData &>(mapping_data);
2428
2429 switch (mapping_kind)
2430 {
2431 case mapping_covariant:
2432 {
2435 "update_covariant_transformation"));
2436
2437 for (unsigned int i = 0; i < output.size(); ++i)
2438 output[i] = apply_transformation(data.covariant[i], input[i]);
2439
2440 return;
2441 }
2442 default:
2443 Assert(false, ExcNotImplemented());
2444 }
2445 }
2446 } // namespace MappingQGenericImplementation
2447} // namespace internal
2448
2450
2451#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:697
DerivativeForm< 1, spacedim, dim, Number > transpose() const
Abstract base class for mapping classes.
Definition: mapping.h:304
const Point< dim > & point(const unsigned int i) const
const std::vector< double > & get_weights() const
unsigned int size() const
Definition: tensor.h:472
numbers::NumberTraits< Number >::real_type norm() const
constexpr void clear()
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
std::vector< Tensor< 1, spacedim > > normal_vectors
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< Point< spacedim > > quadrature_points
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< Tensor< 1, spacedim > > boundary_forms
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
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
InverseQuadraticApproximation(const InverseQuadraticApproximation &)=default
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
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)
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_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
Point< 2 > second
Definition: grid_out.cc:4588
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:439
MappingKind
Definition: mapping.h:65
@ mapping_piola
Definition: mapping.h:100
@ mapping_covariant_gradient
Definition: mapping.h:86
@ mapping_covariant
Definition: mapping.h:75
@ mapping_contravariant
Definition: mapping.h:80
@ mapping_contravariant_hessian
Definition: mapping.h:142
@ mapping_covariant_hessian
Definition: mapping.h:136
@ mapping_contravariant_gradient
Definition: mapping.h:92
@ mapping_piola_gradient
Definition: mapping.h:106
@ mapping_piola_hessian
Definition: mapping.h:148
LogStream deallog
Definition: logstream.cc:37
Expression fabs(const Expression &x)
EvaluationFlags
The EvaluationFlags enum.
@ matrix
Contents is actually a matrix.
static const char A
static const char T
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
Definition: generators.cc:451
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
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 QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, 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 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)
Point< spacedim > compute_mapped_location_of_point(const typename ::MappingQGeneric< dim, spacedim >::InternalData &data)
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)
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_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim > > &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< Tensor< 5, spacedim > > &jacobian_pushed_forward_3rd_derivatives)
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data)
inline ::Table< 2, double > compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
Point< dim > do_transform_real_to_unit_cell_internal_codim1(const typename ::Triangulation< dim, dim+1 >::cell_iterator &cell, const Point< dim+1 > &p, const Point< dim > &initial_p_unit, typename ::MappingQGeneric< dim, dim+1 >::InternalData &mdata)
void do_fill_fe_face_values(const ::MappingQGeneric< 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 ::MappingQGeneric< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
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 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_q_points_Jacobians_and_grads_tensor(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< Point< spacedim > > &quadrature_points, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
void maybe_compute_face_data(const ::MappingQGeneric< 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 ::MappingQGeneric< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
inline ::Table< 2, double > compute_support_point_weights_cell(const unsigned int polynomial_degree)
void maybe_compute_q_points(const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< Point< spacedim > > &quadrature_points)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim > > &jacobian_2nd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim > > &jacobian_pushed_forward_grads)
inline ::Table< 2, double > compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
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={})
static const unsigned int invalid_unsigned_int
Definition: types.h:196
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2539
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2564