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