Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55: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 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_mapping_q_internal_h
17 #define dealii_mapping_q_internal_h
18 
19 #include <deal.II/base/config.h>
20 
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 typename ::Triangulation<dim, dim + 1>::cell_iterator &cell,
737  const Point<dim + 1> & p,
738  const Point<dim> & initial_p_unit,
739  typename ::MappingQ<dim, dim + 1>::InternalData &mdata)
740  {
741  const unsigned int spacedim = dim + 1;
742 
743  const unsigned int n_shapes = mdata.shape_values.size();
744  (void)n_shapes;
745  Assert(n_shapes != 0, ExcInternalError());
746  Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
747  Assert(mdata.shape_second_derivatives.size() == n_shapes,
748  ExcInternalError());
749 
750  std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
751  Assert(points.size() == n_shapes, ExcInternalError());
752 
753  Point<spacedim> p_minus_F;
754 
755  Tensor<1, spacedim> DF[dim];
756  Tensor<1, spacedim> D2F[dim][dim];
757 
758  Point<dim> p_unit = initial_p_unit;
759  Point<dim> f;
760  Tensor<2, dim> df;
761 
762  // Evaluate first and second derivatives
763  mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
764 
765  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
766  {
767  const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
768  const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
769  const Point<spacedim> &point_k = points[k];
770 
771  for (unsigned int j = 0; j < dim; ++j)
772  {
773  DF[j] += grad_phi_k[j] * point_k;
774  for (unsigned int l = 0; l < dim; ++l)
775  D2F[j][l] += hessian_k[j][l] * point_k;
776  }
777  }
778 
779  p_minus_F = p;
780  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
781 
782 
783  for (unsigned int j = 0; j < dim; ++j)
784  f[j] = DF[j] * p_minus_F;
785 
786  for (unsigned int j = 0; j < dim; ++j)
787  {
788  f[j] = DF[j] * p_minus_F;
789  for (unsigned int l = 0; l < dim; ++l)
790  df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
791  }
792 
793 
794  const double eps = 1.e-12 * cell->diameter();
795  const unsigned int loop_limit = 10;
796 
797  unsigned int loop = 0;
798 
799  while (f.norm() > eps && loop++ < loop_limit)
800  {
801  // Solve [df(x)]d=f(x)
802  const Tensor<1, dim> d =
803  invert(df) * static_cast<const Tensor<1, dim> &>(f);
804  p_unit -= d;
805 
806  for (unsigned int j = 0; j < dim; ++j)
807  {
808  DF[j].clear();
809  for (unsigned int l = 0; l < dim; ++l)
810  D2F[j][l].clear();
811  }
812 
813  mdata.compute_shape_function_values(
814  std::vector<Point<dim>>(1, p_unit));
815 
816  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
817  {
818  const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
819  const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
820  const Point<spacedim> &point_k = points[k];
821 
822  for (unsigned int j = 0; j < dim; ++j)
823  {
824  DF[j] += grad_phi_k[j] * point_k;
825  for (unsigned int l = 0; l < dim; ++l)
826  D2F[j][l] += hessian_k[j][l] * point_k;
827  }
828  }
829 
830  // TODO: implement a line search here in much the same way as for
831  // the corresponding function above that does so for dim==spacedim
832  p_minus_F = p;
833  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
834 
835  for (unsigned int j = 0; j < dim; ++j)
836  {
837  f[j] = DF[j] * p_minus_F;
838  for (unsigned int l = 0; l < dim; ++l)
839  df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
840  }
841  }
842 
843 
844  // Here we check that in the last execution of while the first
845  // condition was already wrong, meaning the residual was below
846  // eps. Only if the first condition failed, loop will have been
847  // increased and tested, and thus have reached the limit.
848  AssertThrow(loop < loop_limit,
850 
851  return p_unit;
852  }
853 
854 
855 
873  template <int dim, int spacedim>
875  {
876  public:
880  static constexpr unsigned int n_functions =
881  (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
882 
895  const std::vector<Point<spacedim>> &real_support_points,
896  const std::vector<Point<dim>> & unit_support_points)
897  : normalization_shift(real_support_points[0])
899  1. / real_support_points[0].distance(real_support_points[1]))
900  , is_affine(true)
901  {
902  AssertDimension(real_support_points.size(), unit_support_points.size());
903 
904  // For the bi-/trilinear approximation, we cannot build a quadratic
905  // polynomial due to a lack of points (interpolation matrix would
906  // get singular). Similarly, it is not entirely clear how to gather
907  // enough information for the case dim < spacedim.
908  //
909  // In both cases we require the vector real_support_points to
910  // contain the vertex positions and fall back to an affine
911  // approximation:
912  Assert(dim == spacedim || real_support_points.size() ==
914  ExcInternalError());
915  if (real_support_points.size() == GeometryInfo<dim>::vertices_per_cell)
916  {
917  const auto affine = GridTools::affine_cell_approximation<dim>(
918  make_array_view(real_support_points));
920  affine.first.covariant_form().transpose();
921 
922  // The code for evaluation assumes an additional transformation of
923  // the form (x - normalization_shift) * normalization_length --
924  // account for this in the definition of the coefficients.
925  coefficients[0] =
927  for (unsigned int d = 0; d < spacedim; ++d)
928  for (unsigned int e = 0; e < dim; ++e)
929  coefficients[1 + d][e] =
930  A_inv[e][d] * (1.0 / normalization_length);
931  is_affine = true;
932  return;
933  }
934 
936  std::array<double, n_functions> shape_values;
937  for (unsigned int q = 0; q < unit_support_points.size(); ++q)
938  {
939  // Evaluate quadratic shape functions in point, with the
940  // normalization applied in order to avoid roundoff issues with
941  // scaling far away from 1.
942  shape_values[0] = 1.;
943  const Tensor<1, spacedim> p_scaled =
945  (real_support_points[q] - normalization_shift);
946  for (unsigned int d = 0; d < spacedim; ++d)
947  shape_values[1 + d] = p_scaled[d];
948  for (unsigned int d = 0, c = 0; d < spacedim; ++d)
949  for (unsigned int e = 0; e <= d; ++e, ++c)
950  shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
951 
952  // Build lower diagonal of least squares matrix and rhs, the
953  // essential part being that we construct the matrix with the
954  // real points and the right hand side by comparing to the
955  // reference point positions which sets up an inverse
956  // interpolation.
957  for (unsigned int i = 0; i < n_functions; ++i)
958  for (unsigned int j = 0; j < n_functions; ++j)
959  matrix[i][j] += shape_values[i] * shape_values[j];
960  for (unsigned int i = 0; i < n_functions; ++i)
961  coefficients[i] += shape_values[i] * unit_support_points[q];
962  }
963 
964  // Factorize the matrix A = L * L^T in-place with the
965  // Cholesky-Banachiewicz algorithm. The implementation is similar to
966  // FullMatrix::cholesky() but re-implemented to avoid memory
967  // allocations and some unnecessary divisions which we can do here as
968  // we only need to solve with dim right hand sides.
969  for (unsigned int i = 0; i < n_functions; ++i)
970  {
971  double Lij_sum = 0;
972  for (unsigned int j = 0; j < i; ++j)
973  {
974  double Lik_Ljk_sum = 0;
975  for (unsigned int k = 0; k < j; ++k)
976  Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
977  matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
978  Lij_sum += matrix[i][j] * matrix[i][j];
979  }
980  AssertThrow(matrix[i][i] - Lij_sum >= 0,
981  ExcMessage("Matrix of normal equations not positive "
982  "definite"));
983 
984  // Store the inverse in the diagonal since that is the quantity
985  // needed later in the factorization as well as the forward and
986  // backward substitution, minimizing the number of divisions.
987  matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
988  }
989 
990  // Solve lower triangular part, L * y = rhs.
991  for (unsigned int i = 0; i < n_functions; ++i)
992  {
994  for (unsigned int j = 0; j < i; ++j)
995  sum -= matrix[i][j] * coefficients[j];
996  coefficients[i] = sum * matrix[i][i];
997  }
998 
999  // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
1000  for (unsigned int i = n_functions; i > 0;)
1001  {
1002  --i;
1004  for (unsigned int j = i + 1; j < n_functions; ++j)
1005  sum -= matrix[j][i] * coefficients[j];
1006  coefficients[i] = sum * matrix[i][i];
1007  }
1008 
1009  // Check whether the approximation is indeed affine, allowing to
1010  // skip the quadratic terms.
1011  is_affine = true;
1012  for (unsigned int i = dim + 1; i < n_functions; ++i)
1013  if (coefficients[i].norm_square() > 1e-20)
1014  {
1015  is_affine = false;
1016  break;
1017  }
1018  }
1019 
1024  default;
1025 
1029  template <typename Number>
1032  {
1033  Point<dim, Number> result;
1034  for (unsigned int d = 0; d < dim; ++d)
1035  result[d] = coefficients[0][d];
1036 
1037  // Apply the normalization to ensure a good conditioning. Since Number
1038  // might be a vectorized array whereas the normalization is a point of
1039  // doubles, we cannot use the overload of operator- and must instead
1040  // loop over the components of the point.
1041  Point<spacedim, Number> p_scaled;
1042  for (unsigned int d = 0; d < spacedim; ++d)
1043  p_scaled[d] = (p[d] - normalization_shift[d]) * normalization_length;
1044 
1045  for (unsigned int d = 0; d < spacedim; ++d)
1046  result += coefficients[1 + d] * p_scaled[d];
1047 
1048  if (!is_affine)
1049  {
1050  Point<dim, Number> result_affine = result;
1051  for (unsigned int d = 0, c = 0; d < spacedim; ++d)
1052  for (unsigned int e = 0; e <= d; ++e, ++c)
1053  result +=
1054  coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
1055 
1056  // Check if the quadratic approximation ends up considerably
1057  // farther outside the unit cell on some or all SIMD lanes than
1058  // the affine approximation - in that case, we switch those
1059  // components back to the affine approximation. Note that the
1060  // quadratic approximation will grow more quickly away from the
1061  // unit cell. We make the selection for each SIMD lane with a
1062  // ternary operation.
1063  const Number distance_to_unit_cell = result.distance_square(
1065  const Number affine_distance_to_unit_cell =
1066  result_affine.distance_square(
1068  for (unsigned int d = 0; d < dim; ++d)
1069  result[d] = compare_and_apply_mask<SIMDComparison::greater_than>(
1070  distance_to_unit_cell,
1071  affine_distance_to_unit_cell + 0.5,
1072  result_affine[d],
1073  result[d]);
1074  }
1075  return result;
1076  }
1077 
1078  private:
1087 
1091  const double normalization_length;
1092 
1096  std::array<Point<dim>, n_functions> coefficients;
1097 
1104  };
1105 
1106 
1107 
1113  template <int dim, int spacedim>
1114  inline void
1116  const CellSimilarity::Similarity cell_similarity,
1117  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1118  std::vector<Point<spacedim>> & quadrature_points,
1119  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1120  {
1121  const UpdateFlags update_flags = data.update_each;
1122 
1123  using VectorizedArrayType =
1124  typename ::MappingQ<dim,
1125  spacedim>::InternalData::VectorizedArrayType;
1126  const unsigned int n_shape_values = data.n_shape_functions;
1127  const unsigned int n_q_points = data.shape_info.n_q_points;
1128  constexpr unsigned int n_lanes = VectorizedArrayType::size();
1129  constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1130  constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1131 
1132  EvaluationFlags::EvaluationFlags evaluation_flag =
1135  ((cell_similarity != CellSimilarity::translation) &&
1136  (update_flags & update_contravariant_transformation) ?
1139  ((cell_similarity != CellSimilarity::translation) &&
1140  (update_flags & update_jacobian_grads) ?
1143 
1144  Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1145  ExcInternalError());
1146  Assert(!(evaluation_flag & EvaluationFlags::values) ||
1147  n_q_points == quadrature_points.size(),
1148  ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1149  Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1150  data.n_shape_functions > 0,
1151  ExcInternalError());
1152  Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1153  n_q_points == data.contravariant.size(),
1154  ExcDimensionMismatch(n_q_points, data.contravariant.size()));
1155  Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
1156  n_q_points == jacobian_grads.size(),
1157  ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1158 
1159  // shortcut in case we have an identity interpolation and only request
1160  // the quadrature points
1161  if (evaluation_flag == EvaluationFlags::values &&
1162  data.shape_info.element_type ==
1164  {
1165  for (unsigned int q = 0; q < n_q_points; ++q)
1166  quadrature_points[q] =
1167  data.mapping_support_points[data.shape_info
1168  .lexicographic_numbering[q]];
1169  return;
1170  }
1171 
1173 
1174  // prepare arrays
1175  if (evaluation_flag != EvaluationFlags::nothing)
1176  {
1177  eval.set_data_pointers(&data.scratch, n_comp);
1178 
1179  // make sure to initialize on all lanes also when some are unused in
1180  // the code below
1181  for (unsigned int i = 0; i < n_shape_values * n_comp; ++i)
1182  eval.begin_dof_values()[i] = VectorizedArrayType();
1183 
1184  const std::vector<unsigned int> &renumber_to_lexicographic =
1185  data.shape_info.lexicographic_numbering;
1186  for (unsigned int i = 0; i < n_shape_values; ++i)
1187  for (unsigned int d = 0; d < spacedim; ++d)
1188  {
1189  const unsigned int in_comp = d % n_lanes;
1190  const unsigned int out_comp = d / n_lanes;
1191  eval
1192  .begin_dof_values()[out_comp * n_shape_values + i][in_comp] =
1193  data.mapping_support_points[renumber_to_lexicographic[i]][d];
1194  }
1195 
1196  // do the actual tensorized evaluation
1198  n_comp, evaluation_flag, eval.begin_dof_values(), eval);
1199  }
1200 
1201  // do the postprocessing
1202  if (evaluation_flag & EvaluationFlags::values)
1203  {
1204  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1205  for (unsigned int i = 0; i < n_q_points; ++i)
1206  for (unsigned int in_comp = 0;
1207  in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
1208  ++in_comp)
1209  quadrature_points[i][out_comp * n_lanes + in_comp] =
1210  eval.begin_values()[out_comp * n_q_points + i][in_comp];
1211  }
1212 
1213  if (evaluation_flag & EvaluationFlags::gradients)
1214  {
1215  std::fill(data.contravariant.begin(),
1216  data.contravariant.end(),
1218  // We need to reinterpret the data after evaluate has been applied.
1219  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1220  for (unsigned int point = 0; point < n_q_points; ++point)
1221  for (unsigned int j = 0; j < dim; ++j)
1222  for (unsigned int in_comp = 0;
1223  in_comp < n_lanes &&
1224  in_comp < spacedim - out_comp * n_lanes;
1225  ++in_comp)
1226  {
1227  const unsigned int total_number = point * dim + j;
1228  const unsigned int new_comp = total_number / n_q_points;
1229  const unsigned int new_point = total_number % n_q_points;
1230  data.contravariant[new_point][out_comp * n_lanes + in_comp]
1231  [new_comp] =
1232  eval.begin_gradients()[(out_comp * n_q_points + point) *
1233  dim +
1234  j][in_comp];
1235  }
1236  }
1237  if (update_flags & update_covariant_transformation)
1238  if (cell_similarity != CellSimilarity::translation)
1239  for (unsigned int point = 0; point < n_q_points; ++point)
1240  data.covariant[point] =
1241  (data.contravariant[point]).covariant_form();
1242 
1243  if (update_flags & update_volume_elements)
1244  if (cell_similarity != CellSimilarity::translation)
1245  for (unsigned int point = 0; point < n_q_points; ++point)
1246  data.volume_elements[point] =
1247  data.contravariant[point].determinant();
1248 
1249  if (evaluation_flag & EvaluationFlags::hessians)
1250  {
1251  constexpr int desymmetrize_3d[6][2] = {
1252  {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1253  constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1254 
1255  // We need to reinterpret the data after evaluate has been applied.
1256  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1257  for (unsigned int point = 0; point < n_q_points; ++point)
1258  for (unsigned int j = 0; j < n_hessians; ++j)
1259  for (unsigned int in_comp = 0;
1260  in_comp < n_lanes &&
1261  in_comp < spacedim - out_comp * n_lanes;
1262  ++in_comp)
1263  {
1264  const unsigned int total_number = point * n_hessians + j;
1265  const unsigned int new_point = total_number % n_q_points;
1266  const unsigned int new_hessian_comp =
1267  total_number / n_q_points;
1268  const unsigned int new_hessian_comp_i =
1269  dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1270  desymmetrize_3d[new_hessian_comp][0];
1271  const unsigned int new_hessian_comp_j =
1272  dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1273  desymmetrize_3d[new_hessian_comp][1];
1274  const double value =
1275  eval.begin_hessians()[(out_comp * n_q_points + point) *
1276  n_hessians +
1277  j][in_comp];
1278  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1279  [new_hessian_comp_i][new_hessian_comp_j] =
1280  value;
1281  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1282  [new_hessian_comp_j][new_hessian_comp_i] =
1283  value;
1284  }
1285  }
1286  }
1287 
1288 
1295  template <int dim, int spacedim>
1296  inline void
1298  const typename QProjector<dim>::DataSetDescriptor data_set,
1299  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1300  std::vector<Point<spacedim>> &quadrature_points)
1301  {
1302  const UpdateFlags update_flags = data.update_each;
1303 
1304  if (update_flags & update_quadrature_points)
1305  for (unsigned int point = 0; point < quadrature_points.size(); ++point)
1306  {
1307  const double * shape = &data.shape(point + data_set, 0);
1308  Point<spacedim> result =
1309  (shape[0] * data.mapping_support_points[0]);
1310  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1311  for (unsigned int i = 0; i < spacedim; ++i)
1312  result[i] += shape[k] * data.mapping_support_points[k][i];
1313  quadrature_points[point] = result;
1314  }
1315  }
1316 
1317 
1318 
1327  template <int dim, int spacedim>
1328  inline void
1330  const CellSimilarity::Similarity cell_similarity,
1331  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1332  const typename ::MappingQ<dim, spacedim>::InternalData &data)
1333  {
1334  const UpdateFlags update_flags = data.update_each;
1335 
1336  if (update_flags & update_contravariant_transformation)
1337  // if the current cell is just a
1338  // translation of the previous one, no
1339  // need to recompute jacobians...
1340  if (cell_similarity != CellSimilarity::translation)
1341  {
1342  const unsigned int n_q_points = data.contravariant.size();
1343 
1344  std::fill(data.contravariant.begin(),
1345  data.contravariant.end(),
1347 
1348  Assert(data.n_shape_functions > 0, ExcInternalError());
1349 
1350  for (unsigned int point = 0; point < n_q_points; ++point)
1351  {
1352  double result[spacedim][dim];
1353 
1354  // peel away part of sum to avoid zeroing the
1355  // entries and adding for the first time
1356  for (unsigned int i = 0; i < spacedim; ++i)
1357  for (unsigned int j = 0; j < dim; ++j)
1358  result[i][j] = data.derivative(point + data_set, 0)[j] *
1359  data.mapping_support_points[0][i];
1360  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1361  for (unsigned int i = 0; i < spacedim; ++i)
1362  for (unsigned int j = 0; j < dim; ++j)
1363  result[i][j] += data.derivative(point + data_set, k)[j] *
1364  data.mapping_support_points[k][i];
1365 
1366  // write result into contravariant data. for
1367  // j=dim in the case dim<spacedim, there will
1368  // never be any nonzero data that arrives in
1369  // here, so it is ok anyway because it was
1370  // initialized to zero at the initialization
1371  for (unsigned int i = 0; i < spacedim; ++i)
1372  for (unsigned int j = 0; j < dim; ++j)
1373  data.contravariant[point][i][j] = result[i][j];
1374  }
1375  }
1376 
1377  if (update_flags & update_covariant_transformation)
1378  if (cell_similarity != CellSimilarity::translation)
1379  {
1380  const unsigned int n_q_points = data.contravariant.size();
1381  for (unsigned int point = 0; point < n_q_points; ++point)
1382  {
1383  data.covariant[point] =
1384  (data.contravariant[point]).covariant_form();
1385  }
1386  }
1387 
1388  if (update_flags & update_volume_elements)
1389  if (cell_similarity != CellSimilarity::translation)
1390  {
1391  const unsigned int n_q_points = data.contravariant.size();
1392  for (unsigned int point = 0; point < n_q_points; ++point)
1393  data.volume_elements[point] =
1394  data.contravariant[point].determinant();
1395  }
1396  }
1397 
1398 
1399 
1406  template <int dim, int spacedim>
1407  inline void
1409  const CellSimilarity::Similarity cell_similarity,
1410  const typename QProjector<dim>::DataSetDescriptor data_set,
1411  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1412  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1413  {
1414  const UpdateFlags update_flags = data.update_each;
1415  if (update_flags & update_jacobian_grads)
1416  {
1417  const unsigned int n_q_points = jacobian_grads.size();
1418 
1419  if (cell_similarity != CellSimilarity::translation)
1420  for (unsigned int point = 0; point < n_q_points; ++point)
1421  {
1422  const Tensor<2, dim> *second =
1423  &data.second_derivative(point + data_set, 0);
1424  double result[spacedim][dim][dim];
1425  for (unsigned int i = 0; i < spacedim; ++i)
1426  for (unsigned int j = 0; j < dim; ++j)
1427  for (unsigned int l = 0; l < dim; ++l)
1428  result[i][j][l] =
1429  (second[0][j][l] * data.mapping_support_points[0][i]);
1430  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1431  for (unsigned int i = 0; i < spacedim; ++i)
1432  for (unsigned int j = 0; j < dim; ++j)
1433  for (unsigned int l = 0; l < dim; ++l)
1434  result[i][j][l] +=
1435  (second[k][j][l] * data.mapping_support_points[k][i]);
1436 
1437  for (unsigned int i = 0; i < spacedim; ++i)
1438  for (unsigned int j = 0; j < dim; ++j)
1439  for (unsigned int l = 0; l < dim; ++l)
1440  jacobian_grads[point][i][j][l] = result[i][j][l];
1441  }
1442  }
1443  }
1444 
1445 
1446 
1453  template <int dim, int spacedim>
1454  inline void
1456  const CellSimilarity::Similarity cell_similarity,
1457  const typename QProjector<dim>::DataSetDescriptor data_set,
1458  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1459  std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1460  {
1461  const UpdateFlags update_flags = data.update_each;
1462  if (update_flags & update_jacobian_pushed_forward_grads)
1463  {
1464  const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
1465 
1466  if (cell_similarity != CellSimilarity::translation)
1467  {
1468  double tmp[spacedim][spacedim][spacedim];
1469  for (unsigned int point = 0; point < n_q_points; ++point)
1470  {
1471  const Tensor<2, dim> *second =
1472  &data.second_derivative(point + data_set, 0);
1473  double result[spacedim][dim][dim];
1474  for (unsigned int i = 0; i < spacedim; ++i)
1475  for (unsigned int j = 0; j < dim; ++j)
1476  for (unsigned int l = 0; l < dim; ++l)
1477  result[i][j][l] =
1478  (second[0][j][l] * data.mapping_support_points[0][i]);
1479  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1480  for (unsigned int i = 0; i < spacedim; ++i)
1481  for (unsigned int j = 0; j < dim; ++j)
1482  for (unsigned int l = 0; l < dim; ++l)
1483  result[i][j][l] +=
1484  (second[k][j][l] *
1485  data.mapping_support_points[k][i]);
1486 
1487  // first push forward the j-components
1488  for (unsigned int i = 0; i < spacedim; ++i)
1489  for (unsigned int j = 0; j < spacedim; ++j)
1490  for (unsigned int l = 0; l < dim; ++l)
1491  {
1492  tmp[i][j][l] =
1493  result[i][0][l] * data.covariant[point][j][0];
1494  for (unsigned int jr = 1; jr < dim; ++jr)
1495  {
1496  tmp[i][j][l] +=
1497  result[i][jr][l] * data.covariant[point][j][jr];
1498  }
1499  }
1500 
1501  // now, pushing forward the l-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 < spacedim; ++l)
1505  {
1506  jacobian_pushed_forward_grads[point][i][j][l] =
1507  tmp[i][j][0] * data.covariant[point][l][0];
1508  for (unsigned int lr = 1; lr < dim; ++lr)
1509  {
1510  jacobian_pushed_forward_grads[point][i][j][l] +=
1511  tmp[i][j][lr] * data.covariant[point][l][lr];
1512  }
1513  }
1514  }
1515  }
1516  }
1517  }
1518 
1519 
1520 
1527  template <int dim, int spacedim>
1528  inline void
1530  const CellSimilarity::Similarity cell_similarity,
1531  const typename QProjector<dim>::DataSetDescriptor data_set,
1532  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1533  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1534  {
1535  const UpdateFlags update_flags = data.update_each;
1536  if (update_flags & update_jacobian_2nd_derivatives)
1537  {
1538  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1539 
1540  if (cell_similarity != CellSimilarity::translation)
1541  {
1542  for (unsigned int point = 0; point < n_q_points; ++point)
1543  {
1544  const Tensor<3, dim> *third =
1545  &data.third_derivative(point + data_set, 0);
1546  double result[spacedim][dim][dim][dim];
1547  for (unsigned int i = 0; i < spacedim; ++i)
1548  for (unsigned int j = 0; j < dim; ++j)
1549  for (unsigned int l = 0; l < dim; ++l)
1550  for (unsigned int m = 0; m < dim; ++m)
1551  result[i][j][l][m] =
1552  (third[0][j][l][m] *
1553  data.mapping_support_points[0][i]);
1554  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1555  for (unsigned int i = 0; i < spacedim; ++i)
1556  for (unsigned int j = 0; j < dim; ++j)
1557  for (unsigned int l = 0; l < dim; ++l)
1558  for (unsigned int m = 0; m < dim; ++m)
1559  result[i][j][l][m] +=
1560  (third[k][j][l][m] *
1561  data.mapping_support_points[k][i]);
1562 
1563  for (unsigned int i = 0; i < spacedim; ++i)
1564  for (unsigned int j = 0; j < dim; ++j)
1565  for (unsigned int l = 0; l < dim; ++l)
1566  for (unsigned int m = 0; m < dim; ++m)
1567  jacobian_2nd_derivatives[point][i][j][l][m] =
1568  result[i][j][l][m];
1569  }
1570  }
1571  }
1572  }
1573 
1574 
1575 
1583  template <int dim, int spacedim>
1584  inline void
1586  const CellSimilarity::Similarity cell_similarity,
1587  const typename QProjector<dim>::DataSetDescriptor data_set,
1588  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1589  std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
1590  {
1591  const UpdateFlags update_flags = data.update_each;
1593  {
1594  const unsigned int n_q_points =
1595  jacobian_pushed_forward_2nd_derivatives.size();
1596 
1597  if (cell_similarity != CellSimilarity::translation)
1598  {
1599  double tmp[spacedim][spacedim][spacedim][spacedim];
1600  for (unsigned int point = 0; point < n_q_points; ++point)
1601  {
1602  const Tensor<3, dim> *third =
1603  &data.third_derivative(point + data_set, 0);
1604  double result[spacedim][dim][dim][dim];
1605  for (unsigned int i = 0; i < spacedim; ++i)
1606  for (unsigned int j = 0; j < dim; ++j)
1607  for (unsigned int l = 0; l < dim; ++l)
1608  for (unsigned int m = 0; m < dim; ++m)
1609  result[i][j][l][m] =
1610  (third[0][j][l][m] *
1611  data.mapping_support_points[0][i]);
1612  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1613  for (unsigned int i = 0; i < spacedim; ++i)
1614  for (unsigned int j = 0; j < dim; ++j)
1615  for (unsigned int l = 0; l < dim; ++l)
1616  for (unsigned int m = 0; m < dim; ++m)
1617  result[i][j][l][m] +=
1618  (third[k][j][l][m] *
1619  data.mapping_support_points[k][i]);
1620 
1621  // push forward the j-coordinate
1622  for (unsigned int i = 0; i < spacedim; ++i)
1623  for (unsigned int j = 0; j < spacedim; ++j)
1624  for (unsigned int l = 0; l < dim; ++l)
1625  for (unsigned int m = 0; m < dim; ++m)
1626  {
1627  jacobian_pushed_forward_2nd_derivatives
1628  [point][i][j][l][m] = result[i][0][l][m] *
1629  data.covariant[point][j][0];
1630  for (unsigned int jr = 1; jr < dim; ++jr)
1631  jacobian_pushed_forward_2nd_derivatives[point][i]
1632  [j][l]
1633  [m] +=
1634  result[i][jr][l][m] *
1635  data.covariant[point][j][jr];
1636  }
1637 
1638  // push forward the l-coordinate
1639  for (unsigned int i = 0; i < spacedim; ++i)
1640  for (unsigned int j = 0; j < spacedim; ++j)
1641  for (unsigned int l = 0; l < spacedim; ++l)
1642  for (unsigned int m = 0; m < dim; ++m)
1643  {
1644  tmp[i][j][l][m] =
1645  jacobian_pushed_forward_2nd_derivatives[point][i]
1646  [j][0][m] *
1647  data.covariant[point][l][0];
1648  for (unsigned int lr = 1; lr < dim; ++lr)
1649  tmp[i][j][l][m] +=
1650  jacobian_pushed_forward_2nd_derivatives[point]
1651  [i][j]
1652  [lr][m] *
1653  data.covariant[point][l][lr];
1654  }
1655 
1656  // push forward the m-coordinate
1657  for (unsigned int i = 0; i < spacedim; ++i)
1658  for (unsigned int j = 0; j < spacedim; ++j)
1659  for (unsigned int l = 0; l < spacedim; ++l)
1660  for (unsigned int m = 0; m < spacedim; ++m)
1661  {
1662  jacobian_pushed_forward_2nd_derivatives
1663  [point][i][j][l][m] =
1664  tmp[i][j][l][0] * data.covariant[point][m][0];
1665  for (unsigned int mr = 1; mr < dim; ++mr)
1666  jacobian_pushed_forward_2nd_derivatives[point][i]
1667  [j][l]
1668  [m] +=
1669  tmp[i][j][l][mr] * data.covariant[point][m][mr];
1670  }
1671  }
1672  }
1673  }
1674  }
1675 
1676 
1677 
1684  template <int dim, int spacedim>
1685  inline void
1687  const CellSimilarity::Similarity cell_similarity,
1688  const typename QProjector<dim>::DataSetDescriptor data_set,
1689  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1690  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1691  {
1692  const UpdateFlags update_flags = data.update_each;
1693  if (update_flags & update_jacobian_3rd_derivatives)
1694  {
1695  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1696 
1697  if (cell_similarity != CellSimilarity::translation)
1698  {
1699  for (unsigned int point = 0; point < n_q_points; ++point)
1700  {
1701  const Tensor<4, dim> *fourth =
1702  &data.fourth_derivative(point + data_set, 0);
1703  double result[spacedim][dim][dim][dim][dim];
1704  for (unsigned int i = 0; i < spacedim; ++i)
1705  for (unsigned int j = 0; j < dim; ++j)
1706  for (unsigned int l = 0; l < dim; ++l)
1707  for (unsigned int m = 0; m < dim; ++m)
1708  for (unsigned int n = 0; n < dim; ++n)
1709  result[i][j][l][m][n] =
1710  (fourth[0][j][l][m][n] *
1711  data.mapping_support_points[0][i]);
1712  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1713  for (unsigned int i = 0; i < spacedim; ++i)
1714  for (unsigned int j = 0; j < dim; ++j)
1715  for (unsigned int l = 0; l < dim; ++l)
1716  for (unsigned int m = 0; m < dim; ++m)
1717  for (unsigned int n = 0; n < dim; ++n)
1718  result[i][j][l][m][n] +=
1719  (fourth[k][j][l][m][n] *
1720  data.mapping_support_points[k][i]);
1721 
1722  for (unsigned int i = 0; i < spacedim; ++i)
1723  for (unsigned int j = 0; j < dim; ++j)
1724  for (unsigned int l = 0; l < dim; ++l)
1725  for (unsigned int m = 0; m < dim; ++m)
1726  for (unsigned int n = 0; n < dim; ++n)
1727  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1728  result[i][j][l][m][n];
1729  }
1730  }
1731  }
1732  }
1733 
1734 
1735 
1743  template <int dim, int spacedim>
1744  inline void
1746  const CellSimilarity::Similarity cell_similarity,
1747  const typename QProjector<dim>::DataSetDescriptor data_set,
1748  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1749  std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
1750  {
1751  const UpdateFlags update_flags = data.update_each;
1753  {
1754  const unsigned int n_q_points =
1755  jacobian_pushed_forward_3rd_derivatives.size();
1756 
1757  if (cell_similarity != CellSimilarity::translation)
1758  {
1759  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1760  for (unsigned int point = 0; point < n_q_points; ++point)
1761  {
1762  const Tensor<4, dim> *fourth =
1763  &data.fourth_derivative(point + data_set, 0);
1764  double result[spacedim][dim][dim][dim][dim];
1765  for (unsigned int i = 0; i < spacedim; ++i)
1766  for (unsigned int j = 0; j < dim; ++j)
1767  for (unsigned int l = 0; l < dim; ++l)
1768  for (unsigned int m = 0; m < dim; ++m)
1769  for (unsigned int n = 0; n < dim; ++n)
1770  result[i][j][l][m][n] =
1771  (fourth[0][j][l][m][n] *
1772  data.mapping_support_points[0][i]);
1773  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1774  for (unsigned int i = 0; i < spacedim; ++i)
1775  for (unsigned int j = 0; j < dim; ++j)
1776  for (unsigned int l = 0; l < dim; ++l)
1777  for (unsigned int m = 0; m < dim; ++m)
1778  for (unsigned int n = 0; n < dim; ++n)
1779  result[i][j][l][m][n] +=
1780  (fourth[k][j][l][m][n] *
1781  data.mapping_support_points[k][i]);
1782 
1783  // push-forward the j-coordinate
1784  for (unsigned int i = 0; i < spacedim; ++i)
1785  for (unsigned int j = 0; j < spacedim; ++j)
1786  for (unsigned int l = 0; l < dim; ++l)
1787  for (unsigned int m = 0; m < dim; ++m)
1788  for (unsigned int n = 0; n < dim; ++n)
1789  {
1790  tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1791  data.covariant[point][j][0];
1792  for (unsigned int jr = 1; jr < dim; ++jr)
1793  tmp[i][j][l][m][n] +=
1794  result[i][jr][l][m][n] *
1795  data.covariant[point][j][jr];
1796  }
1797 
1798  // push-forward the l-coordinate
1799  for (unsigned int i = 0; i < spacedim; ++i)
1800  for (unsigned int j = 0; j < spacedim; ++j)
1801  for (unsigned int l = 0; l < spacedim; ++l)
1802  for (unsigned int m = 0; m < dim; ++m)
1803  for (unsigned int n = 0; n < dim; ++n)
1804  {
1805  jacobian_pushed_forward_3rd_derivatives
1806  [point][i][j][l][m][n] =
1807  tmp[i][j][0][m][n] *
1808  data.covariant[point][l][0];
1809  for (unsigned int lr = 1; lr < dim; ++lr)
1810  jacobian_pushed_forward_3rd_derivatives[point]
1811  [i][j][l]
1812  [m][n] +=
1813  tmp[i][j][lr][m][n] *
1814  data.covariant[point][l][lr];
1815  }
1816 
1817  // push-forward the m-coordinate
1818  for (unsigned int i = 0; i < spacedim; ++i)
1819  for (unsigned int j = 0; j < spacedim; ++j)
1820  for (unsigned int l = 0; l < spacedim; ++l)
1821  for (unsigned int m = 0; m < spacedim; ++m)
1822  for (unsigned int n = 0; n < dim; ++n)
1823  {
1824  tmp[i][j][l][m][n] =
1825  jacobian_pushed_forward_3rd_derivatives[point]
1826  [i][j][l]
1827  [0][n] *
1828  data.covariant[point][m][0];
1829  for (unsigned int mr = 1; mr < dim; ++mr)
1830  tmp[i][j][l][m][n] +=
1831  jacobian_pushed_forward_3rd_derivatives
1832  [point][i][j][l][mr][n] *
1833  data.covariant[point][m][mr];
1834  }
1835 
1836  // push-forward the n-coordinate
1837  for (unsigned int i = 0; i < spacedim; ++i)
1838  for (unsigned int j = 0; j < spacedim; ++j)
1839  for (unsigned int l = 0; l < spacedim; ++l)
1840  for (unsigned int m = 0; m < spacedim; ++m)
1841  for (unsigned int n = 0; n < spacedim; ++n)
1842  {
1843  jacobian_pushed_forward_3rd_derivatives
1844  [point][i][j][l][m][n] =
1845  tmp[i][j][l][m][0] *
1846  data.covariant[point][n][0];
1847  for (unsigned int nr = 1; nr < dim; ++nr)
1848  jacobian_pushed_forward_3rd_derivatives[point]
1849  [i][j][l]
1850  [m][n] +=
1851  tmp[i][j][l][m][nr] *
1852  data.covariant[point][n][nr];
1853  }
1854  }
1855  }
1856  }
1857  }
1858 
1859 
1860 
1870  template <int dim, int spacedim>
1871  inline void
1873  const ::MappingQ<dim, spacedim> &mapping,
1874  const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
1875  const unsigned int face_no,
1876  const unsigned int subface_no,
1877  const unsigned int n_q_points,
1878  const std::vector<double> & weights,
1879  const typename ::MappingQ<dim, spacedim>::InternalData &data,
1881  &output_data)
1882  {
1883  const UpdateFlags update_flags = data.update_each;
1884 
1885  if (update_flags &
1888  {
1889  if (update_flags & update_boundary_forms)
1890  AssertDimension(output_data.boundary_forms.size(), n_q_points);
1891  if (update_flags & update_normal_vectors)
1892  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1893  if (update_flags & update_JxW_values)
1894  AssertDimension(output_data.JxW_values.size(), n_q_points);
1895 
1896  Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1897 
1898  // first compute some common data that is used for evaluating
1899  // all of the flags below
1900 
1901  // map the unit tangentials to the real cell. checking for d!=dim-1
1902  // eliminates compiler warnings regarding unsigned int expressions <
1903  // 0.
1904  for (unsigned int d = 0; d != dim - 1; ++d)
1905  {
1907  data.unit_tangentials.size(),
1908  ExcInternalError());
1909  Assert(
1910  data.aux[d].size() <=
1911  data
1912  .unit_tangentials[face_no +
1914  .size(),
1915  ExcInternalError());
1916 
1917  mapping.transform(
1919  data.unit_tangentials[face_no +
1922  data,
1923  make_array_view(data.aux[d].begin(), data.aux[d].end()));
1924  }
1925 
1926  if (update_flags & update_boundary_forms)
1927  {
1928  // if dim==spacedim, we can use the unit tangentials to compute
1929  // the boundary form by simply taking the cross product
1930  if (dim == spacedim)
1931  {
1932  for (unsigned int i = 0; i < n_q_points; ++i)
1933  switch (dim)
1934  {
1935  case 1:
1936  // in 1d, we don't have access to any of the
1937  // data.aux fields (because it has only dim-1
1938  // components), but we can still compute the
1939  // boundary form by simply looking at the number of
1940  // the face
1941  output_data.boundary_forms[i][0] =
1942  (face_no == 0 ? -1 : +1);
1943  break;
1944  case 2:
1945  output_data.boundary_forms[i] =
1946  cross_product_2d(data.aux[0][i]);
1947  break;
1948  case 3:
1949  output_data.boundary_forms[i] =
1950  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1951  break;
1952  default:
1953  Assert(false, ExcNotImplemented());
1954  }
1955  }
1956  else //(dim < spacedim)
1957  {
1958  // in the codim-one case, the boundary form results from the
1959  // cross product of all the face tangential vectors and the
1960  // cell normal vector
1961  //
1962  // to compute the cell normal, use the same method used in
1963  // fill_fe_values for cells above
1964  AssertDimension(data.contravariant.size(), n_q_points);
1965 
1966  for (unsigned int point = 0; point < n_q_points; ++point)
1967  {
1968  if (dim == 1)
1969  {
1970  // J is a tangent vector
1971  output_data.boundary_forms[point] =
1972  data.contravariant[point].transpose()[0];
1973  output_data.boundary_forms[point] /=
1974  (face_no == 0 ? -1. : +1.) *
1975  output_data.boundary_forms[point].norm();
1976  }
1977 
1978  if (dim == 2)
1979  {
1981  data.contravariant[point].transpose();
1982 
1983  Tensor<1, spacedim> cell_normal =
1984  cross_product_3d(DX_t[0], DX_t[1]);
1985  cell_normal /= cell_normal.norm();
1986 
1987  // then compute the face normal from the face
1988  // tangent and the cell normal:
1989  output_data.boundary_forms[point] =
1990  cross_product_3d(data.aux[0][point], cell_normal);
1991  }
1992  }
1993  }
1994  }
1995 
1996  if (update_flags & update_JxW_values)
1997  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
1998  {
1999  output_data.JxW_values[i] =
2000  output_data.boundary_forms[i].norm() * weights[i];
2001 
2002  if (subface_no != numbers::invalid_unsigned_int)
2003  {
2004  const double area_ratio = GeometryInfo<dim>::subface_ratio(
2005  cell->subface_case(face_no), subface_no);
2006  output_data.JxW_values[i] *= area_ratio;
2007  }
2008  }
2009 
2010  if (update_flags & update_normal_vectors)
2011  for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
2012  output_data.normal_vectors[i] =
2013  Point<spacedim>(output_data.boundary_forms[i] /
2014  output_data.boundary_forms[i].norm());
2015 
2016  if (update_flags & update_jacobians)
2017  for (unsigned int point = 0; point < n_q_points; ++point)
2018  output_data.jacobians[point] = data.contravariant[point];
2019 
2020  if (update_flags & update_inverse_jacobians)
2021  for (unsigned int point = 0; point < n_q_points; ++point)
2022  output_data.inverse_jacobians[point] =
2023  data.covariant[point].transpose();
2024  }
2025  }
2026 
2027 
2034  template <int dim, int spacedim>
2035  inline void
2037  const ::MappingQ<dim, spacedim> &mapping,
2038  const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
2039  const unsigned int face_no,
2040  const unsigned int subface_no,
2041  const typename QProjector<dim>::DataSetDescriptor data_set,
2042  const Quadrature<dim - 1> & quadrature,
2043  const typename ::MappingQ<dim, spacedim>::InternalData &data,
2045  &output_data)
2046  {
2047  if (dim > 1 && data.tensor_product_quadrature)
2048  {
2049  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2051  data,
2052  output_data.quadrature_points,
2053  output_data.jacobian_grads);
2054  }
2055  else
2056  {
2057  maybe_compute_q_points<dim, spacedim>(data_set,
2058  data,
2059  output_data.quadrature_points);
2060  maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
2061  data_set,
2062  data);
2063  maybe_update_jacobian_grads<dim, spacedim>(
2064  CellSimilarity::none, data_set, data, output_data.jacobian_grads);
2065  }
2066  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2068  data_set,
2069  data,
2070  output_data.jacobian_pushed_forward_grads);
2071  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2073  data_set,
2074  data,
2075  output_data.jacobian_2nd_derivatives);
2076  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2078  data_set,
2079  data,
2081  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2083  data_set,
2084  data,
2085  output_data.jacobian_3rd_derivatives);
2086  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2088  data_set,
2089  data,
2091 
2092  maybe_compute_face_data(mapping,
2093  cell,
2094  face_no,
2095  subface_no,
2096  quadrature.size(),
2097  quadrature.get_weights(),
2098  data,
2099  output_data);
2100  }
2101 
2102 
2103 
2107  template <int dim, int spacedim, int rank>
2108  inline void
2110  const ArrayView<const Tensor<rank, dim>> & input,
2111  const MappingKind mapping_kind,
2112  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2113  const ArrayView<Tensor<rank, spacedim>> & output)
2114  {
2115  AssertDimension(input.size(), output.size());
2116  Assert((dynamic_cast<
2117  const typename ::MappingQ<dim, spacedim>::InternalData *>(
2118  &mapping_data) != nullptr),
2119  ExcInternalError());
2120  const typename ::MappingQ<dim, spacedim>::InternalData &data =
2121  static_cast<
2122  const typename ::MappingQ<dim, spacedim>::InternalData &>(
2123  mapping_data);
2124 
2125  switch (mapping_kind)
2126  {
2127  case mapping_contravariant:
2128  {
2129  Assert(data.update_each & update_contravariant_transformation,
2131  "update_contravariant_transformation"));
2132 
2133  for (unsigned int i = 0; i < output.size(); ++i)
2134  output[i] =
2135  apply_transformation(data.contravariant[i], input[i]);
2136 
2137  return;
2138  }
2139 
2140  case mapping_piola:
2141  {
2142  Assert(data.update_each & update_contravariant_transformation,
2144  "update_contravariant_transformation"));
2145  Assert(data.update_each & update_volume_elements,
2147  "update_volume_elements"));
2148  Assert(rank == 1, ExcMessage("Only for rank 1"));
2149  if (rank != 1)
2150  return;
2151 
2152  for (unsigned int i = 0; i < output.size(); ++i)
2153  {
2154  output[i] =
2155  apply_transformation(data.contravariant[i], input[i]);
2156  output[i] /= data.volume_elements[i];
2157  }
2158  return;
2159  }
2160  // We still allow this operation as in the
2161  // reference cell Derivatives are Tensor
2162  // rather than DerivativeForm
2163  case mapping_covariant:
2164  {
2165  Assert(data.update_each & update_contravariant_transformation,
2167  "update_covariant_transformation"));
2168 
2169  for (unsigned int i = 0; i < output.size(); ++i)
2170  output[i] = apply_transformation(data.covariant[i], input[i]);
2171 
2172  return;
2173  }
2174 
2175  default:
2176  Assert(false, ExcNotImplemented());
2177  }
2178  }
2179 
2180 
2181 
2185  template <int dim, int spacedim, int rank>
2186  inline void
2188  const ArrayView<const Tensor<rank, dim>> & input,
2189  const MappingKind mapping_kind,
2190  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2191  const ArrayView<Tensor<rank, spacedim>> & output)
2192  {
2193  AssertDimension(input.size(), output.size());
2194  Assert((dynamic_cast<
2195  const typename ::MappingQ<dim, spacedim>::InternalData *>(
2196  &mapping_data) != nullptr),
2197  ExcInternalError());
2198  const typename ::MappingQ<dim, spacedim>::InternalData &data =
2199  static_cast<
2200  const typename ::MappingQ<dim, spacedim>::InternalData &>(
2201  mapping_data);
2202 
2203  switch (mapping_kind)
2204  {
2206  {
2207  Assert(data.update_each & update_covariant_transformation,
2209  "update_covariant_transformation"));
2210  Assert(data.update_each & update_contravariant_transformation,
2212  "update_contravariant_transformation"));
2213  Assert(rank == 2, ExcMessage("Only for rank 2"));
2214 
2215  for (unsigned int i = 0; i < output.size(); ++i)
2216  {
2218  apply_transformation(data.contravariant[i],
2219  transpose(input[i]));
2220  output[i] =
2221  apply_transformation(data.covariant[i], A.transpose());
2222  }
2223 
2224  return;
2225  }
2226 
2228  {
2229  Assert(data.update_each & update_covariant_transformation,
2231  "update_covariant_transformation"));
2232  Assert(rank == 2, ExcMessage("Only for rank 2"));
2233 
2234  for (unsigned int i = 0; i < output.size(); ++i)
2235  {
2237  apply_transformation(data.covariant[i],
2238  transpose(input[i]));
2239  output[i] =
2240  apply_transformation(data.covariant[i], A.transpose());
2241  }
2242 
2243  return;
2244  }
2245 
2247  {
2248  Assert(data.update_each & update_covariant_transformation,
2250  "update_covariant_transformation"));
2251  Assert(data.update_each & update_contravariant_transformation,
2253  "update_contravariant_transformation"));
2254  Assert(data.update_each & update_volume_elements,
2256  "update_volume_elements"));
2257  Assert(rank == 2, ExcMessage("Only for rank 2"));
2258 
2259  for (unsigned int i = 0; i < output.size(); ++i)
2260  {
2262  apply_transformation(data.covariant[i], input[i]);
2263  const Tensor<2, spacedim> T =
2264  apply_transformation(data.contravariant[i], A.transpose());
2265 
2266  output[i] = transpose(T);
2267  output[i] /= data.volume_elements[i];
2268  }
2269 
2270  return;
2271  }
2272 
2273  default:
2274  Assert(false, ExcNotImplemented());
2275  }
2276  }
2277 
2278 
2279 
2283  template <int dim, int spacedim>
2284  inline void
2286  const ArrayView<const Tensor<3, dim>> & input,
2287  const MappingKind mapping_kind,
2288  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2289  const ArrayView<Tensor<3, spacedim>> & output)
2290  {
2291  AssertDimension(input.size(), output.size());
2292  Assert((dynamic_cast<
2293  const typename ::MappingQ<dim, spacedim>::InternalData *>(
2294  &mapping_data) != nullptr),
2295  ExcInternalError());
2296  const typename ::MappingQ<dim, spacedim>::InternalData &data =
2297  static_cast<
2298  const typename ::MappingQ<dim, spacedim>::InternalData &>(
2299  mapping_data);
2300 
2301  switch (mapping_kind)
2302  {
2304  {
2305  Assert(data.update_each & update_covariant_transformation,
2307  "update_covariant_transformation"));
2308  Assert(data.update_each & update_contravariant_transformation,
2310  "update_contravariant_transformation"));
2311 
2312  for (unsigned int q = 0; q < output.size(); ++q)
2313  for (unsigned int i = 0; i < spacedim; ++i)
2314  {
2315  double tmp1[dim][dim];
2316  for (unsigned int J = 0; J < dim; ++J)
2317  for (unsigned int K = 0; K < dim; ++K)
2318  {
2319  tmp1[J][K] =
2320  data.contravariant[q][i][0] * input[q][0][J][K];
2321  for (unsigned int I = 1; I < dim; ++I)
2322  tmp1[J][K] +=
2323  data.contravariant[q][i][I] * input[q][I][J][K];
2324  }
2325  for (unsigned int j = 0; j < spacedim; ++j)
2326  {
2327  double tmp2[dim];
2328  for (unsigned int K = 0; K < dim; ++K)
2329  {
2330  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2331  for (unsigned int J = 1; J < dim; ++J)
2332  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2333  }
2334  for (unsigned int k = 0; k < spacedim; ++k)
2335  {
2336  output[q][i][j][k] =
2337  data.covariant[q][k][0] * tmp2[0];
2338  for (unsigned int K = 1; K < dim; ++K)
2339  output[q][i][j][k] +=
2340  data.covariant[q][k][K] * tmp2[K];
2341  }
2342  }
2343  }
2344  return;
2345  }
2346 
2348  {
2349  Assert(data.update_each & update_covariant_transformation,
2351  "update_covariant_transformation"));
2352 
2353  for (unsigned int q = 0; q < output.size(); ++q)
2354  for (unsigned int i = 0; i < spacedim; ++i)
2355  {
2356  double tmp1[dim][dim];
2357  for (unsigned int J = 0; J < dim; ++J)
2358  for (unsigned int K = 0; K < dim; ++K)
2359  {
2360  tmp1[J][K] =
2361  data.covariant[q][i][0] * input[q][0][J][K];
2362  for (unsigned int I = 1; I < dim; ++I)
2363  tmp1[J][K] +=
2364  data.covariant[q][i][I] * input[q][I][J][K];
2365  }
2366  for (unsigned int j = 0; j < spacedim; ++j)
2367  {
2368  double tmp2[dim];
2369  for (unsigned int K = 0; K < dim; ++K)
2370  {
2371  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2372  for (unsigned int J = 1; J < dim; ++J)
2373  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2374  }
2375  for (unsigned int k = 0; k < spacedim; ++k)
2376  {
2377  output[q][i][j][k] =
2378  data.covariant[q][k][0] * tmp2[0];
2379  for (unsigned int K = 1; K < dim; ++K)
2380  output[q][i][j][k] +=
2381  data.covariant[q][k][K] * tmp2[K];
2382  }
2383  }
2384  }
2385 
2386  return;
2387  }
2388 
2389  case mapping_piola_hessian:
2390  {
2391  Assert(data.update_each & update_covariant_transformation,
2393  "update_covariant_transformation"));
2394  Assert(data.update_each & update_contravariant_transformation,
2396  "update_contravariant_transformation"));
2397  Assert(data.update_each & update_volume_elements,
2399  "update_volume_elements"));
2400 
2401  for (unsigned int q = 0; q < output.size(); ++q)
2402  for (unsigned int i = 0; i < spacedim; ++i)
2403  {
2404  double factor[dim];
2405  for (unsigned int I = 0; I < dim; ++I)
2406  factor[I] =
2407  data.contravariant[q][i][I] / data.volume_elements[q];
2408  double tmp1[dim][dim];
2409  for (unsigned int J = 0; J < dim; ++J)
2410  for (unsigned int K = 0; K < dim; ++K)
2411  {
2412  tmp1[J][K] = factor[0] * input[q][0][J][K];
2413  for (unsigned int I = 1; I < dim; ++I)
2414  tmp1[J][K] += factor[I] * input[q][I][J][K];
2415  }
2416  for (unsigned int j = 0; j < spacedim; ++j)
2417  {
2418  double tmp2[dim];
2419  for (unsigned int K = 0; K < dim; ++K)
2420  {
2421  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2422  for (unsigned int J = 1; J < dim; ++J)
2423  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2424  }
2425  for (unsigned int k = 0; k < spacedim; ++k)
2426  {
2427  output[q][i][j][k] =
2428  data.covariant[q][k][0] * tmp2[0];
2429  for (unsigned int K = 1; K < dim; ++K)
2430  output[q][i][j][k] +=
2431  data.covariant[q][k][K] * tmp2[K];
2432  }
2433  }
2434  }
2435 
2436  return;
2437  }
2438 
2439  default:
2440  Assert(false, ExcNotImplemented());
2441  }
2442  }
2443 
2444 
2445 
2450  template <int dim, int spacedim, int rank>
2451  inline void
2453  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2454  const MappingKind mapping_kind,
2455  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_data,
2456  const ArrayView<Tensor<rank + 1, spacedim>> & output)
2457  {
2458  AssertDimension(input.size(), output.size());
2459  Assert((dynamic_cast<
2460  const typename ::MappingQ<dim, spacedim>::InternalData *>(
2461  &mapping_data) != nullptr),
2462  ExcInternalError());
2463  const typename ::MappingQ<dim, spacedim>::InternalData &data =
2464  static_cast<
2465  const typename ::MappingQ<dim, spacedim>::InternalData &>(
2466  mapping_data);
2467 
2468  switch (mapping_kind)
2469  {
2470  case mapping_covariant:
2471  {
2472  Assert(data.update_each & update_contravariant_transformation,
2474  "update_covariant_transformation"));
2475 
2476  for (unsigned int i = 0; i < output.size(); ++i)
2477  output[i] = apply_transformation(data.covariant[i], input[i]);
2478 
2479  return;
2480  }
2481  default:
2482  Assert(false, ExcNotImplemented());
2483  }
2484  }
2485  } // namespace MappingQImplementation
2486 } // namespace internal
2487 
2489 
2490 #endif
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:704
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:314
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 void clear()
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_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Point< 2 > second
Definition: grid_out.cc:4616
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
MappingKind
Definition: mapping.h:75
@ mapping_piola
Definition: mapping.h:110
@ mapping_covariant_gradient
Definition: mapping.h:96
@ mapping_covariant
Definition: mapping.h:85
@ mapping_contravariant
Definition: mapping.h:90
@ mapping_contravariant_hessian
Definition: mapping.h:152
@ mapping_covariant_hessian
Definition: mapping.h:146
@ mapping_contravariant_gradient
Definition: mapping.h:102
@ mapping_piola_gradient
Definition: mapping.h:116
@ mapping_piola_hessian
Definition: mapping.h:158
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:189
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:493
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)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
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)
void maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 5, spacedim >> &jacobian_pushed_forward_3rd_derivatives)
Point< spacedim > compute_mapped_location_of_point(const typename ::MappingQ< dim, spacedim >::InternalData &data)
void maybe_update_q_points_Jacobians_and_grads_tensor(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Point< spacedim >> &quadrature_points, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
Point< dim > do_transform_real_to_unit_cell_internal_codim1(const typename ::Triangulation< dim, dim+1 >::cell_iterator &cell, const Point< dim+1 > &p, const Point< dim > &initial_p_unit, typename ::MappingQ< dim, dim+1 >::InternalData &mdata)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 4, spacedim >> &jacobian_pushed_forward_2nd_derivatives)
void maybe_compute_q_points(const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Point< spacedim >> &quadrature_points)
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data)
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
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_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
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)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition: types.h:213
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 > &)