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