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