Reference documentation for deal.II version Git 170b4c9308 2021-10-26 16:43:28 -0600
\(\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 
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 MappingQImplementation
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 ::MappingQ<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 ::MappingQ<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 ::MappingQ<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 ::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
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:451
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.
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
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:4587
inline ::Table< 2, double > compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
MappingKind
Definition: mapping.h:64
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:2650
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:303
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:698
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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={})
Definition: tensor.h:506
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
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()
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 > &)
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:2625
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
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()