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