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