Reference documentation for deal.II version Git 53084d8c6a 2020-10-21 00:34:33 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mapping_q_generic.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 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 
23 #include <deal.II/base/table.h>
25 
26 #include <deal.II/fe/fe_dgq.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping_q1.h>
31 
34 #include <deal.II/grid/tria.h>
36 
39 
43 
44 #include <boost/container/small_vector.hpp>
45 
46 #include <algorithm>
47 #include <array>
48 #include <cmath>
49 #include <memory>
50 #include <numeric>
51 
52 
54 
55 
56 namespace internal
57 {
58  namespace MappingQ1
59  {
60  namespace
61  {
62  // These are left as templates on the spatial dimension (even though dim
63  // == spacedim must be true for them to make sense) because templates are
64  // expanded before the compiler eliminates code due to the 'if (dim ==
65  // spacedim)' statement (see the body of the general
66  // transform_real_to_unit_cell).
67  template <int spacedim>
68  Point<1>
69  transform_real_to_unit_cell(
71  & vertices,
72  const Point<spacedim> &p)
73  {
74  Assert(spacedim == 1, ExcInternalError());
75  return Point<1>((p[0] - vertices[0](0)) /
76  (vertices[1](0) - vertices[0](0)));
77  }
78 
79 
80 
81  template <int spacedim>
82  Point<2>
83  transform_real_to_unit_cell(
85  & vertices,
86  const Point<spacedim> &p)
87  {
88  Assert(spacedim == 2, ExcInternalError());
89 
90  // For accuracy reasons, we do all arithmetic in extended precision
91  // (long double). This has a noticeable effect on the hit rate for
92  // borderline cases and thus makes the algorithm more robust.
93  const long double x = p(0);
94  const long double y = p(1);
95 
96  const long double x0 = vertices[0](0);
97  const long double x1 = vertices[1](0);
98  const long double x2 = vertices[2](0);
99  const long double x3 = vertices[3](0);
100 
101  const long double y0 = vertices[0](1);
102  const long double y1 = vertices[1](1);
103  const long double y2 = vertices[2](1);
104  const long double y3 = vertices[3](1);
105 
106  const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3);
107  const long double b = -(x0 - x1 - x2 + x3) * y +
108  (x - 2 * x1 + x3) * y0 - (x - 2 * x0 + x2) * y1 -
109  (x - x1) * y2 + (x - x0) * y3;
110  const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1;
111 
112  const long double discriminant = b * b - 4 * a * c;
113  // exit if the point is not in the cell (this is the only case where the
114  // discriminant is negative)
115  AssertThrow(
116  discriminant > 0.0,
118 
119  long double eta1;
120  long double eta2;
121  const long double sqrt_discriminant = std::sqrt(discriminant);
122  // special case #1: if a is near-zero to make the discriminant exactly
123  // equal b, then use the linear formula
124  if (b != 0.0 && std::abs(b) == sqrt_discriminant)
125  {
126  eta1 = -c / b;
127  eta2 = -c / b;
128  }
129  // special case #2: a is zero for parallelograms and very small for
130  // near-parallelograms:
131  else if (std::abs(a) < 1e-8 * std::abs(b))
132  {
133  // if both a and c are very small then the root should be near
134  // zero: this first case will capture that
135  eta1 = 2 * c / (-b - sqrt_discriminant);
136  eta2 = 2 * c / (-b + sqrt_discriminant);
137  }
138  // finally, use the plain version:
139  else
140  {
141  eta1 = (-b - sqrt_discriminant) / (2 * a);
142  eta2 = (-b + sqrt_discriminant) / (2 * a);
143  }
144  // pick the one closer to the center of the cell.
145  const long double eta =
146  (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
147 
148  /*
149  * There are two ways to compute xi from eta, but either one may have a
150  * zero denominator.
151  */
152  const long double subexpr0 = -eta * x2 + x0 * (eta - 1);
153  const long double xi_denominator0 =
154  eta * x3 - x1 * (eta - 1) + subexpr0;
155  const long double max_x =
156  std::max(std::max(std::abs(x0), std::abs(x1)),
157  std::max(std::abs(x2), std::abs(x3)));
158 
159  if (std::abs(xi_denominator0) > 1e-10 * max_x)
160  {
161  const double xi = (x + subexpr0) / xi_denominator0;
162  return {xi, static_cast<double>(eta)};
163  }
164  else
165  {
166  const long double max_y =
167  std::max(std::max(std::abs(y0), std::abs(y1)),
168  std::max(std::abs(y2), std::abs(y3)));
169  const long double subexpr1 = -eta * y2 + y0 * (eta - 1);
170  const long double xi_denominator1 =
171  eta * y3 - y1 * (eta - 1) + subexpr1;
172  if (std::abs(xi_denominator1) > 1e-10 * max_y)
173  {
174  const double xi = (subexpr1 + y) / xi_denominator1;
175  return {xi, static_cast<double>(eta)};
176  }
177  else // give up and try Newton iteration
178  {
179  AssertThrow(
180  false,
181  (typename Mapping<spacedim,
182  spacedim>::ExcTransformationFailed()));
183  }
184  }
185  // bogus return to placate compiler. It should not be possible to get
186  // here.
187  Assert(false, ExcInternalError());
188  return {std::numeric_limits<double>::quiet_NaN(),
189  std::numeric_limits<double>::quiet_NaN()};
190  }
191 
192 
193 
194  template <int spacedim>
195  Point<3>
196  transform_real_to_unit_cell(
198  & /*vertices*/,
199  const Point<spacedim> & /*p*/)
200  {
201  // It should not be possible to get here
202  Assert(false, ExcInternalError());
203  return Point<3>();
204  }
205  } // namespace
206  } // namespace MappingQ1
207 } // namespace internal
208 
209 
210 
211 template <int dim, int spacedim>
213  const unsigned int polynomial_degree)
214  : polynomial_degree(polynomial_degree)
215  , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
216  , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
217  , tensor_product_quadrature(false)
218 {}
219 
220 
221 
222 template <int dim, int spacedim>
223 std::size_t
225 {
226  return (
239 }
240 
241 
242 template <int dim, int spacedim>
243 void
245  const UpdateFlags update_flags,
246  const Quadrature<dim> &q,
247  const unsigned int n_original_q_points)
248 {
249  // store the flags in the internal data object so we can access them
250  // in fill_fe_*_values()
251  this->update_each = update_flags;
252 
253  const unsigned int n_q_points = q.size();
254 
255  const bool needs_higher_order_terms =
256  this->update_each &
261 
263  covariant.resize(n_original_q_points);
264 
266  contravariant.resize(n_original_q_points);
267 
269  volume_elements.resize(n_original_q_points);
270 
272 
273  // use of MatrixFree only for higher order elements and with more than one
274  // point where tensor products do not make sense
275  if (polynomial_degree < 2 || n_q_points == 1)
277 
278  if (dim > 1)
279  {
280  // find out if the one-dimensional formula is the same
281  // in all directions
283  {
284  const std::array<Quadrature<1>, dim> quad_array =
285  q.get_tensor_basis();
286  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
287  {
288  if (quad_array[i - 1].size() != quad_array[i].size())
289  {
290  tensor_product_quadrature = false;
291  break;
292  }
293  else
294  {
295  const std::vector<Point<1>> &points_1 =
296  quad_array[i - 1].get_points();
297  const std::vector<Point<1>> &points_2 =
298  quad_array[i].get_points();
299  const std::vector<double> &weights_1 =
300  quad_array[i - 1].get_weights();
301  const std::vector<double> &weights_2 =
302  quad_array[i].get_weights();
303  for (unsigned int j = 0; j < quad_array[i].size(); ++j)
304  {
305  if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
306  std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
307  {
308  tensor_product_quadrature = false;
309  break;
310  }
311  }
312  }
313  }
314 
315  if (tensor_product_quadrature)
316  {
317  // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
318  // numbering manually (building an FE_Q<dim> is relatively
319  // expensive due to constraints)
320  const FE_DGQ<1> fe(polynomial_degree);
321  shape_info.reinit(q.get_tensor_basis()[0], fe);
323  FETools::lexicographic_to_hierarchic_numbering<dim>(
325  shape_info.n_q_points = q.size();
328  }
329  }
330  }
331 
332  // Only fill the big arrays on demand in case we cannot use the tensor
333  // product quadrature code path
334  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
335  {
336  // see if we need the (transformation) shape function values
337  // and/or gradients and resize the necessary arrays
339  shape_values.resize(n_shape_functions * n_q_points);
340 
341  if (this->update_each &
342  (update_covariant_transformation |
343  update_contravariant_transformation | update_JxW_values |
351  shape_derivatives.resize(n_shape_functions * n_q_points);
352 
353  if (this->update_each &
355  shape_second_derivatives.resize(n_shape_functions * n_q_points);
356 
359  shape_third_derivatives.resize(n_shape_functions * n_q_points);
360 
363  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
364 
365  // now also fill the various fields with their correct values
367  }
368 }
369 
370 
371 
372 template <int dim, int spacedim>
373 void
375  const UpdateFlags update_flags,
376  const Quadrature<dim> &q,
377  const unsigned int n_original_q_points)
378 {
379  initialize(update_flags, q, n_original_q_points);
380 
381  if (dim > 1 && tensor_product_quadrature)
382  {
383  constexpr unsigned int facedim = dim - 1;
384  const FE_DGQ<1> fe(polynomial_degree);
385  shape_info.reinit(q.get_tensor_basis()[0], fe);
387  FETools::lexicographic_to_hierarchic_numbering<facedim>(
389  shape_info.n_q_points = n_original_q_points;
392  }
393 
394  if (dim > 1)
395  {
396  if (this->update_each &
399  {
400  aux.resize(dim - 1,
401  std::vector<Tensor<1, spacedim>>(n_original_q_points));
402 
403  // Compute tangentials to the unit cell.
404  for (const unsigned int i : GeometryInfo<dim>::face_indices())
405  {
406  unit_tangentials[i].resize(n_original_q_points);
407  std::fill(unit_tangentials[i].begin(),
408  unit_tangentials[i].end(),
410  if (dim > 2)
411  {
413  .resize(n_original_q_points);
414  std::fill(
416  .begin(),
418  .end(),
420  }
421  }
422  }
423  }
424 }
425 
426 
427 
428 template <int dim, int spacedim>
429 void
431  const std::vector<Point<dim>> &unit_points)
432 {
433  const unsigned int n_points = unit_points.size();
434 
435  // Construct the tensor product polynomials used as shape functions for
436  // the Qp mapping of cells at the boundary.
437  const TensorProductPolynomials<dim> tensor_pols(
440  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
441 
442  // then also construct the mapping from lexicographic to the Qp shape
443  // function numbering
444  const std::vector<unsigned int> renumber =
445  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
446 
447  std::vector<double> values;
448  std::vector<Tensor<1, dim>> grads;
449  if (shape_values.size() != 0)
450  {
451  Assert(shape_values.size() == n_shape_functions * n_points,
452  ExcInternalError());
453  values.resize(n_shape_functions);
454  }
455  if (shape_derivatives.size() != 0)
456  {
457  Assert(shape_derivatives.size() == n_shape_functions * n_points,
458  ExcInternalError());
459  grads.resize(n_shape_functions);
460  }
461 
462  std::vector<Tensor<2, dim>> grad2;
463  if (shape_second_derivatives.size() != 0)
464  {
466  ExcInternalError());
467  grad2.resize(n_shape_functions);
468  }
469 
470  std::vector<Tensor<3, dim>> grad3;
471  if (shape_third_derivatives.size() != 0)
472  {
473  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
474  ExcInternalError());
475  grad3.resize(n_shape_functions);
476  }
477 
478  std::vector<Tensor<4, dim>> grad4;
479  if (shape_fourth_derivatives.size() != 0)
480  {
482  ExcInternalError());
483  grad4.resize(n_shape_functions);
484  }
485 
486 
487  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
488  shape_second_derivatives.size() != 0 ||
489  shape_third_derivatives.size() != 0 ||
490  shape_fourth_derivatives.size() != 0)
491  for (unsigned int point = 0; point < n_points; ++point)
492  {
493  tensor_pols.evaluate(
494  unit_points[point], values, grads, grad2, grad3, grad4);
495 
496  if (shape_values.size() != 0)
497  for (unsigned int i = 0; i < n_shape_functions; ++i)
498  shape(point, i) = values[renumber[i]];
499 
500  if (shape_derivatives.size() != 0)
501  for (unsigned int i = 0; i < n_shape_functions; ++i)
502  derivative(point, i) = grads[renumber[i]];
503 
504  if (shape_second_derivatives.size() != 0)
505  for (unsigned int i = 0; i < n_shape_functions; ++i)
506  second_derivative(point, i) = grad2[renumber[i]];
507 
508  if (shape_third_derivatives.size() != 0)
509  for (unsigned int i = 0; i < n_shape_functions; ++i)
510  third_derivative(point, i) = grad3[renumber[i]];
511 
512  if (shape_fourth_derivatives.size() != 0)
513  for (unsigned int i = 0; i < n_shape_functions; ++i)
514  fourth_derivative(point, i) = grad4[renumber[i]];
515  }
516 }
517 
518 
519 
520 namespace internal
521 {
522  namespace MappingQGenericImplementation
523  {
524  namespace
525  {
534  compute_support_point_weights_on_quad(
535  const unsigned int polynomial_degree)
536  {
537  ::Table<2, double> loqvs;
538 
539  // we are asked to compute weights for interior support points, but
540  // there are no interior points if degree==1
541  if (polynomial_degree == 1)
542  return loqvs;
543 
544  const unsigned int M = polynomial_degree - 1;
545  const unsigned int n_inner_2d = M * M;
546  const unsigned int n_outer_2d = 4 + 4 * M;
547 
548  // set the weights of transfinite interpolation
549  loqvs.reinit(n_inner_2d, n_outer_2d);
550  QGaussLobatto<2> gl(polynomial_degree + 1);
551  for (unsigned int i = 0; i < M; ++i)
552  for (unsigned int j = 0; j < M; ++j)
553  {
554  const Point<2> p =
555  gl.point((i + 1) * (polynomial_degree + 1) + (j + 1));
556  const unsigned int index_table = i * M + j;
557  for (unsigned int v = 0; v < 4; ++v)
558  loqvs(index_table, v) =
560  loqvs(index_table, 4 + i) = 1. - p[0];
561  loqvs(index_table, 4 + i + M) = p[0];
562  loqvs(index_table, 4 + j + 2 * M) = 1. - p[1];
563  loqvs(index_table, 4 + j + 3 * M) = p[1];
564  }
565 
566  // the sum of weights of the points at the outer rim should be one.
567  // check this
568  for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
569  Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
570  loqvs[unit_point].end(),
571  0.) -
572  1) < 1e-13 * polynomial_degree,
573  ExcInternalError());
574 
575  return loqvs;
576  }
577 
578 
579 
587  compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
588  {
589  ::Table<2, double> lohvs;
590 
591  // we are asked to compute weights for interior support points, but
592  // there are no interior points if degree==1
593  if (polynomial_degree == 1)
594  return lohvs;
595 
596  const unsigned int M = polynomial_degree - 1;
597 
598  const unsigned int n_inner = Utilities::fixed_power<3>(M);
599  const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
600 
601  // set the weights of transfinite interpolation
602  lohvs.reinit(n_inner, n_outer);
603  QGaussLobatto<3> gl(polynomial_degree + 1);
604  for (unsigned int i = 0; i < M; ++i)
605  for (unsigned int j = 0; j < M; ++j)
606  for (unsigned int k = 0; k < M; ++k)
607  {
608  const Point<3> p = gl.point((i + 1) * (M + 2) * (M + 2) +
609  (j + 1) * (M + 2) + (k + 1));
610  const unsigned int index_table = i * M * M + j * M + k;
611 
612  // vertices
613  for (unsigned int v = 0; v < 8; ++v)
614  lohvs(index_table, v) =
616 
617  // lines
618  {
619  constexpr std::array<unsigned int, 4> line_coordinates_y(
620  {{0, 1, 4, 5}});
621  const Point<2> py(p[0], p[2]);
622  for (unsigned int l = 0; l < 4; ++l)
623  lohvs(index_table, 8 + line_coordinates_y[l] * M + j) =
625  }
626 
627  {
628  constexpr std::array<unsigned int, 4> line_coordinates_x(
629  {{2, 3, 6, 7}});
630  const Point<2> px(p[1], p[2]);
631  for (unsigned int l = 0; l < 4; ++l)
632  lohvs(index_table, 8 + line_coordinates_x[l] * M + k) =
634  }
635 
636  {
637  constexpr std::array<unsigned int, 4> line_coordinates_z(
638  {{8, 9, 10, 11}});
639  const Point<2> pz(p[0], p[1]);
640  for (unsigned int l = 0; l < 4; ++l)
641  lohvs(index_table, 8 + line_coordinates_z[l] * M + i) =
643  }
644 
645  // quads
646  lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
647  1. - p[0];
648  lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0];
649  lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) =
650  1. - p[1];
651  lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1];
652  lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) =
653  1. - p[2];
654  lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
655  }
656 
657  // the sum of weights of the points at the outer rim should be one.
658  // check this
659  for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
660  Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
661  lohvs[unit_point].end(),
662  0.) -
663  1) < 1e-13 * polynomial_degree,
664  ExcInternalError());
665 
666  return lohvs;
667  }
668 
669 
670 
675  std::vector<::Table<2, double>>
676  compute_support_point_weights_perimeter_to_interior(
677  const unsigned int polynomial_degree,
678  const unsigned int dim)
679  {
680  Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
681  std::vector<::Table<2, double>> output(dim);
682  if (polynomial_degree <= 1)
683  return output;
684 
685  // fill the 1D interior weights
686  QGaussLobatto<1> quadrature(polynomial_degree + 1);
687  output[0].reinit(polynomial_degree - 1,
689  for (unsigned int q = 0; q < polynomial_degree - 1; ++q)
690  for (const unsigned int i : GeometryInfo<1>::vertex_indices())
691  output[0](q, i) =
693  i);
694 
695  if (dim > 1)
696  output[1] = compute_support_point_weights_on_quad(polynomial_degree);
697 
698  if (dim > 2)
699  output[2] = compute_support_point_weights_on_hex(polynomial_degree);
700 
701  return output;
702  }
703 
707  template <int dim>
709  compute_support_point_weights_cell(const unsigned int polynomial_degree)
710  {
711  Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
712  if (polynomial_degree <= 1)
713  return ::Table<2, double>();
714 
715  QGaussLobatto<dim> quadrature(polynomial_degree + 1);
716  const std::vector<unsigned int> h2l =
717  FETools::hierarchic_to_lexicographic_numbering<dim>(
719 
720  ::Table<2, double> output(quadrature.size() -
723  for (unsigned int q = 0; q < output.size(0); ++q)
724  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
726  quadrature.point(h2l[q + GeometryInfo<dim>::vertices_per_cell]),
727  i);
728 
729  return output;
730  }
731 
732 
733 
741  template <int dim, int spacedim>
743  compute_mapped_location_of_point(
744  const typename ::MappingQGeneric<dim, spacedim>::InternalData
745  &data)
746  {
747  AssertDimension(data.shape_values.size(),
748  data.mapping_support_points.size());
749 
750  // use now the InternalData to compute the point in real space.
751  Point<spacedim> p_real;
752  for (unsigned int i = 0; i < data.mapping_support_points.size(); ++i)
753  p_real += data.mapping_support_points[i] * data.shape(0, i);
754 
755  return p_real;
756  }
757 
758 
763  template <int dim, int spacedim, typename Number>
765  do_transform_real_to_unit_cell_internal(
766  const Point<spacedim, Number> & p,
767  const Point<dim, Number> & initial_p_unit,
768  const std::vector<Point<spacedim>> & points,
770  const std::vector<unsigned int> & renumber)
771  {
772  AssertDimension(points.size(),
773  Utilities::pow(polynomials_1d.size(), dim));
774 
775  // Newton iteration to solve
776  // f(x)=p(x)-p=0
777  // where we are looking for 'x' and p(x) is the forward transformation
778  // from unit to real cell. We solve this using a Newton iteration
779  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
780  // The start value is set to be the linear approximation to the cell
781 
782  // The shape values and derivatives of the mapping at this point are
783  // previously computed.
784 
785  Point<dim, Number> p_unit = initial_p_unit;
787  polynomials_1d, points, p_unit, polynomials_1d.size() == 2, renumber);
788 
789  Tensor<1, spacedim, Number> f = p_real.first - p;
790 
791  // early out if we already have our point in all SIMD lanes, i.e.,
792  // f.norm_square() < 1e-24 * p_real.second[0].norm_square(). To enable
793  // this step also for VectorizedArray where we do not have operator <,
794  // compare the result to zero which is defined for SIMD types
795  if (std::max(Number(0.),
796  f.norm_square() -
797  1e-24 * p_real.second[0].norm_square()) == Number(0.))
798  return p_unit;
799 
800  // we need to compare the position of the computed p(x) against the
801  // given point 'p'. We will terminate the iteration and return 'x' if
802  // they are less than eps apart. The question is how to choose eps --
803  // or, put maybe more generally: in which norm we want these 'p' and
804  // 'p(x)' to be eps apart.
805  //
806  // the question is difficult since we may have to deal with very
807  // elongated cells where we may achieve 1e-12*h for the distance of
808  // these two points in the 'long' direction, but achieving this
809  // tolerance in the 'short' direction of the cell may not be possible
810  //
811  // what we do instead is then to terminate iterations if
812  // \| p(x) - p \|_A < eps
813  // where the A-norm is somehow induced by the transformation of the
814  // cell. in particular, we want to measure distances relative to the
815  // sizes of the cell in its principal directions.
816  //
817  // to define what exactly A should be, note that to first order we have
818  // the following (assuming that x* is the solution of the problem, i.e.,
819  // p(x*)=p):
820  // p(x) - p = p(x) - p(x*)
821  // = -grad p(x) * (x*-x) + higher order terms
822  // This suggest to measure with a norm that corresponds to
823  // A = {[grad p(x]^T [grad p(x)]}^{-1}
824  // because then
825  // \| p(x) - p \|_A \approx \| x - x* \|
826  // Consequently, we will try to enforce that
827  // \| p(x) - p \|_A = \| f \| <= eps
828  //
829  // Note that using this norm is a bit dangerous since the norm changes
830  // in every iteration (A isn't fixed by depending on xk). However, if
831  // the cell is not too deformed (it may be stretched, but not twisted)
832  // then the mapping is almost linear and A is indeed constant or
833  // nearly so.
834  const double eps = 1.e-11;
835  const unsigned int newton_iteration_limit = 20;
836 
837  Point<dim, Number> invalid_point;
838  invalid_point[0] = std::numeric_limits<double>::infinity();
839 
840  unsigned int newton_iteration = 0;
841  Number last_f_weighted_norm_square = 0.;
842  do
843  {
844 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
845  std::cout << "Newton iteration " << newton_iteration << std::endl;
846 #endif
847 
848  // f'(x)
850  for (unsigned int d = 0; d < spacedim; ++d)
851  for (unsigned int e = 0; e < dim; ++e)
852  df[d][e] = p_real.second[e][d];
853 
854  // check if the determinant is positive on all SIMD lanes
855  if (!(std::min(determinant(df),
858  return invalid_point;
859 
860  // Solve [f'(x)]d=f(x)
861  const Tensor<2, spacedim, Number> df_inverse = invert(df);
862  const Tensor<1, spacedim, Number> delta = df_inverse * f;
863  last_f_weighted_norm_square = (df_inverse * f).norm_square();
864 
865 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
866  std::cout << " delta=" << delta << std::endl;
867 #endif
868 
869  // do a line search
870  double step_length = 1;
871  do
872  {
873  // update of p_unit. The spacedim-th component of transformed
874  // point is simply ignored in codimension one case. When this
875  // component is not zero, then we are projecting the point to
876  // the surface or curve identified by the cell.
877  Point<dim, Number> p_unit_trial = p_unit;
878  for (unsigned int i = 0; i < dim; ++i)
879  p_unit_trial[i] -= step_length * delta[i];
880 
881  // shape values and derivatives at new p_unit point
882  const auto p_real_trial =
885  points,
886  p_unit_trial,
887  polynomials_1d.size() == 2,
888  renumber);
889  const Tensor<1, spacedim, Number> f_trial =
890  p_real_trial.first - p;
891 
892 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
893  std::cout << " step_length=" << step_length << std::endl
894  << " ||f || =" << f.norm() << std::endl
895  << " ||f*|| =" << f_trial.norm() << std::endl
896  << " ||f*||_A ="
897  << (df_inverse * f_trial).norm() << std::endl;
898 #endif
899 
900  // see if we are making progress with the current step length
901  // and if not, reduce it by a factor of two and try again
902  //
903  // strictly speaking, we should probably use the same norm as we
904  // use for the outer algorithm. in practice, line search is just
905  // a crutch to find a "reasonable" step length, and so using the
906  // l2 norm is probably just fine
907  //
908  // due to the possible use of VectorizedArray<double>, we must
909  // turn the check f_trial.norm() < f.norm() into a more
910  // complicated statement. We are done if either
911  // last_f_weighted_norm_square is less than the Newton
912  // tolerance (i.e., that particular SIMD lane is already
913  // converged in the previous the Newton iteration, so we might
914  // not be able to decrease the right hand side norm any more)
915  // or if the norm did not increase in the line search
916  if (std::max(last_f_weighted_norm_square - eps * eps,
917  Number(0.)) *
918  std::max(f_trial.norm_square() - f.norm_square(),
919  Number(0.)) ==
920  Number(0.))
921  {
922  p_real = p_real_trial;
923  p_unit = p_unit_trial;
924  f = f_trial;
925  break;
926  }
927  else if (step_length > 0.05)
928  step_length *= 0.5;
929  else
930  return invalid_point;
931  }
932  while (true);
933 
934  ++newton_iteration;
935  if (newton_iteration > newton_iteration_limit)
936  return invalid_point;
937  }
938  while (std::max(eps * eps - last_f_weighted_norm_square, Number(0.)) ==
939  Number(0.));
940 
941  return p_unit;
942  }
943 
944 
945 
949  template <int dim>
950  Point<dim>
951  do_transform_real_to_unit_cell_internal_codim1(
952  const typename ::Triangulation<dim, dim + 1>::cell_iterator &cell,
953  const Point<dim + 1> & p,
954  const Point<dim> &initial_p_unit,
955  typename ::MappingQGeneric<dim, dim + 1>::InternalData &mdata)
956  {
957  const unsigned int spacedim = dim + 1;
958 
959  const unsigned int n_shapes = mdata.shape_values.size();
960  (void)n_shapes;
961  Assert(n_shapes != 0, ExcInternalError());
962  Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
963  Assert(mdata.shape_second_derivatives.size() == n_shapes,
964  ExcInternalError());
965 
966  std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
967  Assert(points.size() == n_shapes, ExcInternalError());
968 
969  Point<spacedim> p_minus_F;
970 
971  Tensor<1, spacedim> DF[dim];
972  Tensor<1, spacedim> D2F[dim][dim];
973 
974  Point<dim> p_unit = initial_p_unit;
975  Point<dim> f;
976  Tensor<2, dim> df;
977 
978  // Evaluate first and second derivatives
979  mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
980 
981  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
982  {
983  const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
984  const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
985  const Point<spacedim> &point_k = points[k];
986 
987  for (unsigned int j = 0; j < dim; ++j)
988  {
989  DF[j] += grad_phi_k[j] * point_k;
990  for (unsigned int l = 0; l < dim; ++l)
991  D2F[j][l] += hessian_k[j][l] * point_k;
992  }
993  }
994 
995  p_minus_F = p;
996  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
997 
998 
999  for (unsigned int j = 0; j < dim; ++j)
1000  f[j] = DF[j] * p_minus_F;
1001 
1002  for (unsigned int j = 0; j < dim; ++j)
1003  {
1004  f[j] = DF[j] * p_minus_F;
1005  for (unsigned int l = 0; l < dim; ++l)
1006  df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
1007  }
1008 
1009 
1010  const double eps = 1.e-12 * cell->diameter();
1011  const unsigned int loop_limit = 10;
1012 
1013  unsigned int loop = 0;
1014 
1015  while (f.norm() > eps && loop++ < loop_limit)
1016  {
1017  // Solve [df(x)]d=f(x)
1018  const Tensor<1, dim> d =
1019  invert(df) * static_cast<const Tensor<1, dim> &>(f);
1020  p_unit -= d;
1021 
1022  for (unsigned int j = 0; j < dim; ++j)
1023  {
1024  DF[j].clear();
1025  for (unsigned int l = 0; l < dim; ++l)
1026  D2F[j][l].clear();
1027  }
1028 
1029  mdata.compute_shape_function_values(
1030  std::vector<Point<dim>>(1, p_unit));
1031 
1032  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1033  {
1034  const Tensor<1, dim> &grad_phi_k = mdata.derivative(0, k);
1035  const Tensor<2, dim> &hessian_k = mdata.second_derivative(0, k);
1036  const Point<spacedim> &point_k = points[k];
1037 
1038  for (unsigned int j = 0; j < dim; ++j)
1039  {
1040  DF[j] += grad_phi_k[j] * point_k;
1041  for (unsigned int l = 0; l < dim; ++l)
1042  D2F[j][l] += hessian_k[j][l] * point_k;
1043  }
1044  }
1045 
1046  // TODO: implement a line search here in much the same way as for
1047  // the corresponding function above that does so for dim==spacedim
1048  p_minus_F = p;
1049  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
1050 
1051  for (unsigned int j = 0; j < dim; ++j)
1052  {
1053  f[j] = DF[j] * p_minus_F;
1054  for (unsigned int l = 0; l < dim; ++l)
1055  df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
1056  }
1057  }
1058 
1059 
1060  // Here we check that in the last execution of while the first
1061  // condition was already wrong, meaning the residual was below
1062  // eps. Only if the first condition failed, loop will have been
1063  // increased and tested, and thus have reached the limit.
1064  AssertThrow(
1065  loop < loop_limit,
1067 
1068  return p_unit;
1069  }
1070 
1076  template <int dim, int spacedim>
1077  void
1078  maybe_update_q_points_Jacobians_and_grads_tensor(
1079  const CellSimilarity::Similarity cell_similarity,
1080  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1081  & data,
1082  std::vector<Point<spacedim>> & quadrature_points,
1083  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1084  {
1085  const UpdateFlags update_flags = data.update_each;
1086 
1087  const unsigned int n_shape_values = data.n_shape_functions;
1088  const unsigned int n_q_points = data.shape_info.n_q_points;
1089  constexpr unsigned int n_lanes = VectorizedArray<double>::size();
1090  constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1091  constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1092 
1093  EvaluationFlags::EvaluationFlags evaluation_flag =
1096  ((cell_similarity != CellSimilarity::translation) &&
1097  (update_flags & update_contravariant_transformation) ?
1100  ((cell_similarity != CellSimilarity::translation) &&
1101  (update_flags & update_jacobian_grads) ?
1104 
1105  Assert(!(evaluation_flag & EvaluationFlags::values) || n_q_points > 0,
1106  ExcInternalError());
1107  Assert(!(evaluation_flag & EvaluationFlags::values) ||
1108  n_q_points == quadrature_points.size(),
1109  ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1110  Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1111  data.n_shape_functions > 0,
1112  ExcInternalError());
1113  Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
1114  n_q_points == data.contravariant.size(),
1115  ExcDimensionMismatch(n_q_points, data.contravariant.size()));
1116  Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
1117  n_q_points == jacobian_grads.size(),
1118  ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1119 
1120  // shortcut in case we have an identity interpolation and only request
1121  // the quadrature points
1122  if (evaluation_flag == EvaluationFlags::values &&
1123  data.shape_info.element_type ==
1125  {
1126  for (unsigned int q = 0; q < n_q_points; ++q)
1127  quadrature_points[q] =
1128  data.mapping_support_points[data.shape_info
1129  .lexicographic_numbering[q]];
1130  return;
1131  }
1132 
1133  // prepare arrays
1134  if (evaluation_flag != EvaluationFlags::nothing)
1135  {
1136  data.values_dofs.resize(n_comp * n_shape_values);
1137  data.values_quad.resize(n_comp * n_q_points);
1138  data.gradients_quad.resize(n_comp * n_q_points * dim);
1139  data.scratch.resize(2 * std::max(n_q_points, n_shape_values));
1140 
1141  if (evaluation_flag & EvaluationFlags::hessians)
1142  data.hessians_quad.resize(n_comp * n_q_points * n_hessians);
1143 
1144  const std::vector<unsigned int> &renumber_to_lexicographic =
1145  data.shape_info.lexicographic_numbering;
1146  for (unsigned int i = 0; i < n_shape_values; ++i)
1147  for (unsigned int d = 0; d < spacedim; ++d)
1148  {
1149  const unsigned int in_comp = d % n_lanes;
1150  const unsigned int out_comp = d / n_lanes;
1151  data.values_dofs[out_comp * n_shape_values + i][in_comp] =
1152  data
1153  .mapping_support_points[renumber_to_lexicographic[i]][d];
1154  }
1155 
1156  // do the actual tensorized evaluation
1158  dim,
1159  double,
1160  VectorizedArray<double>>::evaluate(n_comp,
1161  evaluation_flag,
1162  data.shape_info,
1163  data.values_dofs.begin(),
1164  data.values_quad.begin(),
1165  data.gradients_quad.begin(),
1166  data.hessians_quad.begin(),
1167  data.scratch.begin());
1168  }
1169 
1170  // do the postprocessing
1171  if (evaluation_flag & EvaluationFlags::values)
1172  {
1173  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1174  for (unsigned int i = 0; i < n_q_points; ++i)
1175  for (unsigned int in_comp = 0;
1176  in_comp < n_lanes &&
1177  in_comp < spacedim - out_comp * n_lanes;
1178  ++in_comp)
1179  quadrature_points[i][out_comp * n_lanes + in_comp] =
1180  data.values_quad[out_comp * n_q_points + i][in_comp];
1181  }
1182 
1183  if (evaluation_flag & EvaluationFlags::gradients)
1184  {
1185  std::fill(data.contravariant.begin(),
1186  data.contravariant.end(),
1188  // We need to reinterpret the data after evaluate has been applied.
1189  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1190  for (unsigned int point = 0; point < n_q_points; ++point)
1191  for (unsigned int j = 0; j < dim; ++j)
1192  for (unsigned int in_comp = 0;
1193  in_comp < n_lanes &&
1194  in_comp < spacedim - out_comp * n_lanes;
1195  ++in_comp)
1196  {
1197  const unsigned int total_number = point * dim + j;
1198  const unsigned int new_comp = total_number / n_q_points;
1199  const unsigned int new_point = total_number % n_q_points;
1200  data.contravariant[new_point][out_comp * n_lanes +
1201  in_comp][new_comp] =
1202  data.gradients_quad[(out_comp * n_q_points + point) *
1203  dim +
1204  j][in_comp];
1205  }
1206  }
1207  if (update_flags & update_covariant_transformation)
1208  if (cell_similarity != CellSimilarity::translation)
1209  for (unsigned int point = 0; point < n_q_points; ++point)
1210  data.covariant[point] =
1211  (data.contravariant[point]).covariant_form();
1212 
1213  if (update_flags & update_volume_elements)
1214  if (cell_similarity != CellSimilarity::translation)
1215  for (unsigned int point = 0; point < n_q_points; ++point)
1216  data.volume_elements[point] =
1217  data.contravariant[point].determinant();
1218 
1219  if (evaluation_flag & EvaluationFlags::hessians)
1220  {
1221  constexpr int desymmetrize_3d[6][2] = {
1222  {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1223  constexpr int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1224 
1225  // We need to reinterpret the data after evaluate has been applied.
1226  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1227  for (unsigned int point = 0; point < n_q_points; ++point)
1228  for (unsigned int j = 0; j < n_hessians; ++j)
1229  for (unsigned int in_comp = 0;
1230  in_comp < n_lanes &&
1231  in_comp < spacedim - out_comp * n_lanes;
1232  ++in_comp)
1233  {
1234  const unsigned int total_number = point * n_hessians + j;
1235  const unsigned int new_point = total_number % n_q_points;
1236  const unsigned int new_hessian_comp =
1237  total_number / n_q_points;
1238  const unsigned int new_hessian_comp_i =
1239  dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1240  desymmetrize_3d[new_hessian_comp][0];
1241  const unsigned int new_hessian_comp_j =
1242  dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1243  desymmetrize_3d[new_hessian_comp][1];
1244  const double value =
1245  data.hessians_quad[(out_comp * n_q_points + point) *
1246  n_hessians +
1247  j][in_comp];
1248  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1249  [new_hessian_comp_i][new_hessian_comp_j] =
1250  value;
1251  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1252  [new_hessian_comp_j][new_hessian_comp_i] =
1253  value;
1254  }
1255  }
1256  }
1257 
1258 
1265  template <int dim, int spacedim>
1266  void
1267  maybe_compute_q_points(
1268  const typename QProjector<dim>::DataSetDescriptor data_set,
1269  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1270  & data,
1271  std::vector<Point<spacedim>> &quadrature_points)
1272  {
1273  const UpdateFlags update_flags = data.update_each;
1274 
1275  if (update_flags & update_quadrature_points)
1276  for (unsigned int point = 0; point < quadrature_points.size();
1277  ++point)
1278  {
1279  const double * shape = &data.shape(point + data_set, 0);
1280  Point<spacedim> result =
1281  (shape[0] * data.mapping_support_points[0]);
1282  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1283  for (unsigned int i = 0; i < spacedim; ++i)
1284  result[i] += shape[k] * data.mapping_support_points[k][i];
1285  quadrature_points[point] = result;
1286  }
1287  }
1288 
1289 
1290 
1299  template <int dim, int spacedim>
1300  void
1301  maybe_update_Jacobians(
1302  const CellSimilarity::Similarity cell_similarity,
1303  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1304  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1305  &data)
1306  {
1307  const UpdateFlags update_flags = data.update_each;
1308 
1309  if (update_flags & update_contravariant_transformation)
1310  // if the current cell is just a
1311  // translation of the previous one, no
1312  // need to recompute jacobians...
1313  if (cell_similarity != CellSimilarity::translation)
1314  {
1315  const unsigned int n_q_points = data.contravariant.size();
1316 
1317  std::fill(data.contravariant.begin(),
1318  data.contravariant.end(),
1320 
1321  Assert(data.n_shape_functions > 0, ExcInternalError());
1322 
1323  const Tensor<1, spacedim> *supp_pts =
1324  data.mapping_support_points.data();
1325 
1326  for (unsigned int point = 0; point < n_q_points; ++point)
1327  {
1328  const Tensor<1, dim> *data_derv =
1329  &data.derivative(point + data_set, 0);
1330 
1331  double result[spacedim][dim];
1332 
1333  // peel away part of sum to avoid zeroing the
1334  // entries and adding for the first time
1335  for (unsigned int i = 0; i < spacedim; ++i)
1336  for (unsigned int j = 0; j < dim; ++j)
1337  result[i][j] = data_derv[0][j] * supp_pts[0][i];
1338  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1339  for (unsigned int i = 0; i < spacedim; ++i)
1340  for (unsigned int j = 0; j < dim; ++j)
1341  result[i][j] += data_derv[k][j] * supp_pts[k][i];
1342 
1343  // write result into contravariant data. for
1344  // j=dim in the case dim<spacedim, there will
1345  // never be any nonzero data that arrives in
1346  // here, so it is ok anyway because it was
1347  // initialized to zero at the initialization
1348  for (unsigned int i = 0; i < spacedim; ++i)
1349  for (unsigned int j = 0; j < dim; ++j)
1350  data.contravariant[point][i][j] = result[i][j];
1351  }
1352  }
1353 
1354  if (update_flags & update_covariant_transformation)
1355  if (cell_similarity != CellSimilarity::translation)
1356  {
1357  const unsigned int n_q_points = data.contravariant.size();
1358  for (unsigned int point = 0; point < n_q_points; ++point)
1359  {
1360  data.covariant[point] =
1361  (data.contravariant[point]).covariant_form();
1362  }
1363  }
1364 
1365  if (update_flags & update_volume_elements)
1366  if (cell_similarity != CellSimilarity::translation)
1367  {
1368  const unsigned int n_q_points = data.contravariant.size();
1369  for (unsigned int point = 0; point < n_q_points; ++point)
1370  data.volume_elements[point] =
1371  data.contravariant[point].determinant();
1372  }
1373  }
1374 
1381  template <int dim, int spacedim>
1382  void
1383  maybe_update_jacobian_grads(
1384  const CellSimilarity::Similarity cell_similarity,
1385  const typename QProjector<dim>::DataSetDescriptor data_set,
1386  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1387  & data,
1388  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1389  {
1390  const UpdateFlags update_flags = data.update_each;
1391  if (update_flags & update_jacobian_grads)
1392  {
1393  const unsigned int n_q_points = jacobian_grads.size();
1394 
1395  if (cell_similarity != CellSimilarity::translation)
1396  for (unsigned int point = 0; point < n_q_points; ++point)
1397  {
1398  const Tensor<2, dim> *second =
1399  &data.second_derivative(point + data_set, 0);
1400  double result[spacedim][dim][dim];
1401  for (unsigned int i = 0; i < spacedim; ++i)
1402  for (unsigned int j = 0; j < dim; ++j)
1403  for (unsigned int l = 0; l < dim; ++l)
1404  result[i][j][l] =
1405  (second[0][j][l] * data.mapping_support_points[0][i]);
1406  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1407  for (unsigned int i = 0; i < spacedim; ++i)
1408  for (unsigned int j = 0; j < dim; ++j)
1409  for (unsigned int l = 0; l < dim; ++l)
1410  result[i][j][l] +=
1411  (second[k][j][l] *
1412  data.mapping_support_points[k][i]);
1413 
1414  for (unsigned int i = 0; i < spacedim; ++i)
1415  for (unsigned int j = 0; j < dim; ++j)
1416  for (unsigned int l = 0; l < dim; ++l)
1417  jacobian_grads[point][i][j][l] = result[i][j][l];
1418  }
1419  }
1420  }
1421 
1428  template <int dim, int spacedim>
1429  void
1430  maybe_update_jacobian_pushed_forward_grads(
1431  const CellSimilarity::Similarity cell_similarity,
1432  const typename QProjector<dim>::DataSetDescriptor data_set,
1433  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1434  & data,
1435  std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1436  {
1437  const UpdateFlags update_flags = data.update_each;
1438  if (update_flags & update_jacobian_pushed_forward_grads)
1439  {
1440  const unsigned int n_q_points =
1441  jacobian_pushed_forward_grads.size();
1442 
1443  if (cell_similarity != CellSimilarity::translation)
1444  {
1445  double tmp[spacedim][spacedim][spacedim];
1446  for (unsigned int point = 0; point < n_q_points; ++point)
1447  {
1448  const Tensor<2, dim> *second =
1449  &data.second_derivative(point + data_set, 0);
1450  double result[spacedim][dim][dim];
1451  for (unsigned int i = 0; i < spacedim; ++i)
1452  for (unsigned int j = 0; j < dim; ++j)
1453  for (unsigned int l = 0; l < dim; ++l)
1454  result[i][j][l] = (second[0][j][l] *
1455  data.mapping_support_points[0][i]);
1456  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1457  for (unsigned int i = 0; i < spacedim; ++i)
1458  for (unsigned int j = 0; j < dim; ++j)
1459  for (unsigned int l = 0; l < dim; ++l)
1460  result[i][j][l] +=
1461  (second[k][j][l] *
1462  data.mapping_support_points[k][i]);
1463 
1464  // first push forward the j-components
1465  for (unsigned int i = 0; i < spacedim; ++i)
1466  for (unsigned int j = 0; j < spacedim; ++j)
1467  for (unsigned int l = 0; l < dim; ++l)
1468  {
1469  tmp[i][j][l] =
1470  result[i][0][l] * data.covariant[point][j][0];
1471  for (unsigned int jr = 1; jr < dim; ++jr)
1472  {
1473  tmp[i][j][l] += result[i][jr][l] *
1474  data.covariant[point][j][jr];
1475  }
1476  }
1477 
1478  // now, pushing forward the l-components
1479  for (unsigned int i = 0; i < spacedim; ++i)
1480  for (unsigned int j = 0; j < spacedim; ++j)
1481  for (unsigned int l = 0; l < spacedim; ++l)
1482  {
1483  jacobian_pushed_forward_grads[point][i][j][l] =
1484  tmp[i][j][0] * data.covariant[point][l][0];
1485  for (unsigned int lr = 1; lr < dim; ++lr)
1486  {
1487  jacobian_pushed_forward_grads[point][i][j][l] +=
1488  tmp[i][j][lr] * data.covariant[point][l][lr];
1489  }
1490  }
1491  }
1492  }
1493  }
1494  }
1495 
1502  template <int dim, int spacedim>
1503  void
1504  maybe_update_jacobian_2nd_derivatives(
1505  const CellSimilarity::Similarity cell_similarity,
1506  const typename QProjector<dim>::DataSetDescriptor data_set,
1507  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1508  & data,
1509  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1510  {
1511  const UpdateFlags update_flags = data.update_each;
1512  if (update_flags & update_jacobian_2nd_derivatives)
1513  {
1514  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1515 
1516  if (cell_similarity != CellSimilarity::translation)
1517  {
1518  for (unsigned int point = 0; point < n_q_points; ++point)
1519  {
1520  const Tensor<3, dim> *third =
1521  &data.third_derivative(point + data_set, 0);
1522  double result[spacedim][dim][dim][dim];
1523  for (unsigned int i = 0; i < spacedim; ++i)
1524  for (unsigned int j = 0; j < dim; ++j)
1525  for (unsigned int l = 0; l < dim; ++l)
1526  for (unsigned int m = 0; m < dim; ++m)
1527  result[i][j][l][m] =
1528  (third[0][j][l][m] *
1529  data.mapping_support_points[0][i]);
1530  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1531  for (unsigned int i = 0; i < spacedim; ++i)
1532  for (unsigned int j = 0; j < dim; ++j)
1533  for (unsigned int l = 0; l < dim; ++l)
1534  for (unsigned int m = 0; m < dim; ++m)
1535  result[i][j][l][m] +=
1536  (third[k][j][l][m] *
1537  data.mapping_support_points[k][i]);
1538 
1539  for (unsigned int i = 0; i < spacedim; ++i)
1540  for (unsigned int j = 0; j < dim; ++j)
1541  for (unsigned int l = 0; l < dim; ++l)
1542  for (unsigned int m = 0; m < dim; ++m)
1543  jacobian_2nd_derivatives[point][i][j][l][m] =
1544  result[i][j][l][m];
1545  }
1546  }
1547  }
1548  }
1549 
1557  template <int dim, int spacedim>
1558  void
1559  maybe_update_jacobian_pushed_forward_2nd_derivatives(
1560  const CellSimilarity::Similarity cell_similarity,
1561  const typename QProjector<dim>::DataSetDescriptor data_set,
1562  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1563  &data,
1564  std::vector<Tensor<4, spacedim>>
1565  &jacobian_pushed_forward_2nd_derivatives)
1566  {
1567  const UpdateFlags update_flags = data.update_each;
1569  {
1570  const unsigned int n_q_points =
1571  jacobian_pushed_forward_2nd_derivatives.size();
1572 
1573  if (cell_similarity != CellSimilarity::translation)
1574  {
1575  double tmp[spacedim][spacedim][spacedim][spacedim];
1576  for (unsigned int point = 0; point < n_q_points; ++point)
1577  {
1578  const Tensor<3, dim> *third =
1579  &data.third_derivative(point + data_set, 0);
1580  double result[spacedim][dim][dim][dim];
1581  for (unsigned int i = 0; i < spacedim; ++i)
1582  for (unsigned int j = 0; j < dim; ++j)
1583  for (unsigned int l = 0; l < dim; ++l)
1584  for (unsigned int m = 0; m < dim; ++m)
1585  result[i][j][l][m] =
1586  (third[0][j][l][m] *
1587  data.mapping_support_points[0][i]);
1588  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1589  for (unsigned int i = 0; i < spacedim; ++i)
1590  for (unsigned int j = 0; j < dim; ++j)
1591  for (unsigned int l = 0; l < dim; ++l)
1592  for (unsigned int m = 0; m < dim; ++m)
1593  result[i][j][l][m] +=
1594  (third[k][j][l][m] *
1595  data.mapping_support_points[k][i]);
1596 
1597  // push forward the j-coordinate
1598  for (unsigned int i = 0; i < spacedim; ++i)
1599  for (unsigned int j = 0; j < spacedim; ++j)
1600  for (unsigned int l = 0; l < dim; ++l)
1601  for (unsigned int m = 0; m < dim; ++m)
1602  {
1603  jacobian_pushed_forward_2nd_derivatives
1604  [point][i][j][l][m] =
1605  result[i][0][l][m] *
1606  data.covariant[point][j][0];
1607  for (unsigned int jr = 1; jr < dim; ++jr)
1608  jacobian_pushed_forward_2nd_derivatives[point]
1609  [i][j][l]
1610  [m] +=
1611  result[i][jr][l][m] *
1612  data.covariant[point][j][jr];
1613  }
1614 
1615  // push forward the l-coordinate
1616  for (unsigned int i = 0; i < spacedim; ++i)
1617  for (unsigned int j = 0; j < spacedim; ++j)
1618  for (unsigned int l = 0; l < spacedim; ++l)
1619  for (unsigned int m = 0; m < dim; ++m)
1620  {
1621  tmp[i][j][l][m] =
1622  jacobian_pushed_forward_2nd_derivatives[point]
1623  [i][j][0]
1624  [m] *
1625  data.covariant[point][l][0];
1626  for (unsigned int lr = 1; lr < dim; ++lr)
1627  tmp[i][j][l][m] +=
1628  jacobian_pushed_forward_2nd_derivatives
1629  [point][i][j][lr][m] *
1630  data.covariant[point][l][lr];
1631  }
1632 
1633  // push forward the m-coordinate
1634  for (unsigned int i = 0; i < spacedim; ++i)
1635  for (unsigned int j = 0; j < spacedim; ++j)
1636  for (unsigned int l = 0; l < spacedim; ++l)
1637  for (unsigned int m = 0; m < spacedim; ++m)
1638  {
1639  jacobian_pushed_forward_2nd_derivatives
1640  [point][i][j][l][m] =
1641  tmp[i][j][l][0] * data.covariant[point][m][0];
1642  for (unsigned int mr = 1; mr < dim; ++mr)
1643  jacobian_pushed_forward_2nd_derivatives[point]
1644  [i][j][l]
1645  [m] +=
1646  tmp[i][j][l][mr] *
1647  data.covariant[point][m][mr];
1648  }
1649  }
1650  }
1651  }
1652  }
1653 
1660  template <int dim, int spacedim>
1661  void
1662  maybe_update_jacobian_3rd_derivatives(
1663  const CellSimilarity::Similarity cell_similarity,
1664  const typename QProjector<dim>::DataSetDescriptor data_set,
1665  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1666  & data,
1667  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1668  {
1669  const UpdateFlags update_flags = data.update_each;
1670  if (update_flags & update_jacobian_3rd_derivatives)
1671  {
1672  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1673 
1674  if (cell_similarity != CellSimilarity::translation)
1675  {
1676  for (unsigned int point = 0; point < n_q_points; ++point)
1677  {
1678  const Tensor<4, dim> *fourth =
1679  &data.fourth_derivative(point + data_set, 0);
1680  double result[spacedim][dim][dim][dim][dim];
1681  for (unsigned int i = 0; i < spacedim; ++i)
1682  for (unsigned int j = 0; j < dim; ++j)
1683  for (unsigned int l = 0; l < dim; ++l)
1684  for (unsigned int m = 0; m < dim; ++m)
1685  for (unsigned int n = 0; n < dim; ++n)
1686  result[i][j][l][m][n] =
1687  (fourth[0][j][l][m][n] *
1688  data.mapping_support_points[0][i]);
1689  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1690  for (unsigned int i = 0; i < spacedim; ++i)
1691  for (unsigned int j = 0; j < dim; ++j)
1692  for (unsigned int l = 0; l < dim; ++l)
1693  for (unsigned int m = 0; m < dim; ++m)
1694  for (unsigned int n = 0; n < dim; ++n)
1695  result[i][j][l][m][n] +=
1696  (fourth[k][j][l][m][n] *
1697  data.mapping_support_points[k][i]);
1698 
1699  for (unsigned int i = 0; i < spacedim; ++i)
1700  for (unsigned int j = 0; j < dim; ++j)
1701  for (unsigned int l = 0; l < dim; ++l)
1702  for (unsigned int m = 0; m < dim; ++m)
1703  for (unsigned int n = 0; n < dim; ++n)
1704  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1705  result[i][j][l][m][n];
1706  }
1707  }
1708  }
1709  }
1710 
1718  template <int dim, int spacedim>
1719  void
1720  maybe_update_jacobian_pushed_forward_3rd_derivatives(
1721  const CellSimilarity::Similarity cell_similarity,
1722  const typename QProjector<dim>::DataSetDescriptor data_set,
1723  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1724  &data,
1725  std::vector<Tensor<5, spacedim>>
1726  &jacobian_pushed_forward_3rd_derivatives)
1727  {
1728  const UpdateFlags update_flags = data.update_each;
1730  {
1731  const unsigned int n_q_points =
1732  jacobian_pushed_forward_3rd_derivatives.size();
1733 
1734  if (cell_similarity != CellSimilarity::translation)
1735  {
1736  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1737  for (unsigned int point = 0; point < n_q_points; ++point)
1738  {
1739  const Tensor<4, dim> *fourth =
1740  &data.fourth_derivative(point + data_set, 0);
1741  double result[spacedim][dim][dim][dim][dim];
1742  for (unsigned int i = 0; i < spacedim; ++i)
1743  for (unsigned int j = 0; j < dim; ++j)
1744  for (unsigned int l = 0; l < dim; ++l)
1745  for (unsigned int m = 0; m < dim; ++m)
1746  for (unsigned int n = 0; n < dim; ++n)
1747  result[i][j][l][m][n] =
1748  (fourth[0][j][l][m][n] *
1749  data.mapping_support_points[0][i]);
1750  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
1751  for (unsigned int i = 0; i < spacedim; ++i)
1752  for (unsigned int j = 0; j < dim; ++j)
1753  for (unsigned int l = 0; l < dim; ++l)
1754  for (unsigned int m = 0; m < dim; ++m)
1755  for (unsigned int n = 0; n < dim; ++n)
1756  result[i][j][l][m][n] +=
1757  (fourth[k][j][l][m][n] *
1758  data.mapping_support_points[k][i]);
1759 
1760  // push-forward the j-coordinate
1761  for (unsigned int i = 0; i < spacedim; ++i)
1762  for (unsigned int j = 0; j < spacedim; ++j)
1763  for (unsigned int l = 0; l < dim; ++l)
1764  for (unsigned int m = 0; m < dim; ++m)
1765  for (unsigned int n = 0; n < dim; ++n)
1766  {
1767  tmp[i][j][l][m][n] =
1768  result[i][0][l][m][n] *
1769  data.covariant[point][j][0];
1770  for (unsigned int jr = 1; jr < dim; ++jr)
1771  tmp[i][j][l][m][n] +=
1772  result[i][jr][l][m][n] *
1773  data.covariant[point][j][jr];
1774  }
1775 
1776  // push-forward the l-coordinate
1777  for (unsigned int i = 0; i < spacedim; ++i)
1778  for (unsigned int j = 0; j < spacedim; ++j)
1779  for (unsigned int l = 0; l < spacedim; ++l)
1780  for (unsigned int m = 0; m < dim; ++m)
1781  for (unsigned int n = 0; n < dim; ++n)
1782  {
1783  jacobian_pushed_forward_3rd_derivatives
1784  [point][i][j][l][m][n] =
1785  tmp[i][j][0][m][n] *
1786  data.covariant[point][l][0];
1787  for (unsigned int lr = 1; lr < dim; ++lr)
1788  jacobian_pushed_forward_3rd_derivatives
1789  [point][i][j][l][m][n] +=
1790  tmp[i][j][lr][m][n] *
1791  data.covariant[point][l][lr];
1792  }
1793 
1794  // push-forward the m-coordinate
1795  for (unsigned int i = 0; i < spacedim; ++i)
1796  for (unsigned int j = 0; j < spacedim; ++j)
1797  for (unsigned int l = 0; l < spacedim; ++l)
1798  for (unsigned int m = 0; m < spacedim; ++m)
1799  for (unsigned int n = 0; n < dim; ++n)
1800  {
1801  tmp[i][j][l][m][n] =
1802  jacobian_pushed_forward_3rd_derivatives
1803  [point][i][j][l][0][n] *
1804  data.covariant[point][m][0];
1805  for (unsigned int mr = 1; mr < dim; ++mr)
1806  tmp[i][j][l][m][n] +=
1807  jacobian_pushed_forward_3rd_derivatives
1808  [point][i][j][l][mr][n] *
1809  data.covariant[point][m][mr];
1810  }
1811 
1812  // push-forward the n-coordinate
1813  for (unsigned int i = 0; i < spacedim; ++i)
1814  for (unsigned int j = 0; j < spacedim; ++j)
1815  for (unsigned int l = 0; l < spacedim; ++l)
1816  for (unsigned int m = 0; m < spacedim; ++m)
1817  for (unsigned int n = 0; n < spacedim; ++n)
1818  {
1819  jacobian_pushed_forward_3rd_derivatives
1820  [point][i][j][l][m][n] =
1821  tmp[i][j][l][m][0] *
1822  data.covariant[point][n][0];
1823  for (unsigned int nr = 1; nr < dim; ++nr)
1824  jacobian_pushed_forward_3rd_derivatives
1825  [point][i][j][l][m][n] +=
1826  tmp[i][j][l][m][nr] *
1827  data.covariant[point][n][nr];
1828  }
1829  }
1830  }
1831  }
1832  }
1833  } // namespace
1834  } // namespace MappingQGenericImplementation
1835 } // namespace internal
1836 
1837 
1838 
1839 template <int dim, int spacedim>
1841  : polynomial_degree(p)
1843  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
1844  , polynomials_1d(
1849  internal::MappingQGenericImplementation::
1850  compute_support_point_weights_perimeter_to_interior(
1851  this->polynomial_degree,
1852  dim))
1854  internal::MappingQGenericImplementation::
1855  compute_support_point_weights_cell<dim>(this->polynomial_degree))
1856 {
1857  Assert(p >= 1,
1858  ExcMessage("It only makes sense to create polynomial mappings "
1859  "with a polynomial degree greater or equal to one."));
1860 }
1861 
1862 
1863 
1864 template <int dim, int spacedim>
1866  const MappingQGeneric<dim, spacedim> &mapping)
1869  , polynomials_1d(mapping.polynomials_1d)
1875 {}
1876 
1877 
1878 
1879 template <int dim, int spacedim>
1880 std::unique_ptr<Mapping<dim, spacedim>>
1882 {
1883  return std::make_unique<MappingQGeneric<dim, spacedim>>(*this);
1884 }
1885 
1886 
1887 
1888 template <int dim, int spacedim>
1889 unsigned int
1891 {
1892  return polynomial_degree;
1893 }
1894 
1895 
1896 
1897 template <int dim, int spacedim>
1900  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1901  const Point<dim> & p) const
1902 {
1905  this->compute_mapping_support_points(cell),
1906  p,
1907  polynomials_1d.size() == 2,
1909  .first);
1910 }
1911 
1912 
1913 // In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
1914 // seeing which of the overloaded versions of
1915 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
1916 // error messages and, generally, nothing very good. Avoid this by ensuring
1917 // that this class exists, but does not have an inner InternalData
1918 // type, thereby ruling out the codim-1 version of the function
1919 // below when doing overload resolution.
1920 template <>
1921 class MappingQGeneric<3, 4>
1922 {};
1923 
1924 
1925 
1926 // visual studio freaks out when trying to determine if
1927 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
1928 // candidate. So instead of letting the compiler pick the correct overload, we
1929 // use template specialization to make sure we pick up the right function to
1930 // call:
1931 
1932 template <int dim, int spacedim>
1933 Point<dim>
1936  const Point<spacedim> &,
1937  const Point<dim> &) const
1938 {
1939  // default implementation (should never be called)
1940  Assert(false, ExcInternalError());
1941  return Point<dim>();
1942 }
1943 
1944 
1945 
1946 template <>
1947 Point<1>
1950  const Point<1> & p,
1951  const Point<1> & initial_p_unit) const
1952 {
1953  // dispatch to the various specializations for spacedim=dim,
1954  // spacedim=dim+1, etc
1955  return internal::MappingQGenericImplementation::
1956  do_transform_real_to_unit_cell_internal<1>(
1957  p,
1958  initial_p_unit,
1959  this->compute_mapping_support_points(cell),
1962 }
1963 
1964 
1965 
1966 template <>
1967 Point<2>
1970  const Point<2> & p,
1971  const Point<2> & initial_p_unit) const
1972 {
1973  return internal::MappingQGenericImplementation::
1974  do_transform_real_to_unit_cell_internal<2>(
1975  p,
1976  initial_p_unit,
1977  this->compute_mapping_support_points(cell),
1980 }
1981 
1982 
1983 
1984 template <>
1985 Point<3>
1988  const Point<3> & p,
1989  const Point<3> & initial_p_unit) const
1990 {
1991  return internal::MappingQGenericImplementation::
1992  do_transform_real_to_unit_cell_internal<3>(
1993  p,
1994  initial_p_unit,
1995  this->compute_mapping_support_points(cell),
1998 }
1999 
2000 
2001 
2002 template <>
2003 Point<1>
2006  const Point<2> & p,
2007  const Point<1> & initial_p_unit) const
2008 {
2009  const int dim = 1;
2010  const int spacedim = 2;
2011 
2012  const Quadrature<dim> point_quadrature(initial_p_unit);
2013 
2015  if (spacedim > dim)
2016  update_flags |= update_jacobian_grads;
2018  get_data(update_flags, point_quadrature));
2019 
2021 
2022  // dispatch to the various specializations for spacedim=dim,
2023  // spacedim=dim+1, etc
2024  return internal::MappingQGenericImplementation::
2025  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
2026  p,
2027  initial_p_unit,
2028  *mdata);
2029 }
2030 
2031 
2032 
2033 template <>
2034 Point<2>
2037  const Point<3> & p,
2038  const Point<2> & initial_p_unit) const
2039 {
2040  const int dim = 2;
2041  const int spacedim = 3;
2042 
2043  const Quadrature<dim> point_quadrature(initial_p_unit);
2044 
2046  if (spacedim > dim)
2047  update_flags |= update_jacobian_grads;
2049  get_data(update_flags, point_quadrature));
2050 
2052 
2053  // dispatch to the various specializations for spacedim=dim,
2054  // spacedim=dim+1, etc
2055  return internal::MappingQGenericImplementation::
2056  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
2057  p,
2058  initial_p_unit,
2059  *mdata);
2060 }
2061 
2062 template <>
2063 Point<1>
2066  const Point<3> &,
2067  const Point<1> &) const
2068 {
2069  Assert(false, ExcNotImplemented());
2070  return {};
2071 }
2072 
2073 
2074 
2075 template <int dim, int spacedim>
2076 Point<dim>
2078  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2079  const Point<spacedim> & p) const
2080 {
2081  // Use an exact formula if one is available. this is only the case
2082  // for Q1 mappings in 1d, and in 2d if dim==spacedim
2083  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
2084  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
2085  {
2086  // The dimension-dependent algorithms are much faster (about 25-45x in
2087  // 2D) but fail most of the time when the given point (p) is not in the
2088  // cell. The dimension-independent Newton algorithm given below is
2089  // slower, but more robust (though it still sometimes fails). Therefore
2090  // this function implements the following strategy based on the
2091  // p's dimension:
2092  //
2093  // * In 1D this mapping is linear, so the mapping is always invertible
2094  // (and the exact formula is known) as long as the cell has non-zero
2095  // length.
2096  // * In 2D the exact (quadratic) formula is called first. If either the
2097  // exact formula does not succeed (negative discriminant in the
2098  // quadratic formula) or succeeds but finds a solution outside of the
2099  // unit cell, then the Newton solver is called. The rationale for the
2100  // second choice is that the exact formula may provide two different
2101  // answers when mapping a point outside of the real cell, but the
2102  // Newton solver (if it converges) will only return one answer.
2103  // Otherwise the exact formula successfully found a point in the unit
2104  // cell and that value is returned.
2105  // * In 3D there is no (known to the authors) exact formula, so the Newton
2106  // algorithm is used.
2107  const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
2108  vertices = this->get_vertices(cell);
2109  try
2110  {
2111  switch (dim)
2112  {
2113  case 1:
2114  {
2115  // formula not subject to any issues in 1d
2116  if (spacedim == 1)
2117  return internal::MappingQ1::transform_real_to_unit_cell(
2118  vertices, p);
2119  else
2120  break;
2121  }
2122 
2123  case 2:
2124  {
2125  const Point<dim> point =
2126  internal::MappingQ1::transform_real_to_unit_cell(vertices,
2127  p);
2128 
2129  // formula not guaranteed to work for points outside of
2130  // the cell. only take the computed point if it lies
2131  // inside the reference cell
2132  const double eps = 1e-15;
2133  if (-eps <= point(1) && point(1) <= 1 + eps &&
2134  -eps <= point(0) && point(0) <= 1 + eps)
2135  {
2136  return point;
2137  }
2138  else
2139  break;
2140  }
2141 
2142  default:
2143  {
2144  // we should get here, based on the if-condition at the top
2145  Assert(false, ExcInternalError());
2146  }
2147  }
2148  }
2149  catch (
2151  {
2152  // simply fall through and continue on to the standard Newton code
2153  }
2154  }
2155  else
2156  {
2157  // we can't use an explicit formula,
2158  }
2159 
2160 
2161  // Find the initial value for the Newton iteration by a normal
2162  // projection to the least square plane determined by the vertices
2163  // of the cell
2164  Point<dim> initial_p_unit;
2165  if (this->preserves_vertex_locations())
2166  {
2167  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
2168  // in 1d with spacedim > 1 the affine approximation is exact
2169  if (dim == 1 && polynomial_degree == 1)
2170  return initial_p_unit;
2171  }
2172  else
2173  {
2174  // else, we simply use the mid point
2175  for (unsigned int d = 0; d < dim; ++d)
2176  initial_p_unit[d] = 0.5;
2177  }
2178 
2179  // in case the function above should have given us something back that lies
2180  // outside the unit cell, then project it back into the reference cell in
2181  // hopes that this gives a better starting point to the following iteration
2182  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2183 
2184  // perform the Newton iteration and return the result. note that this
2185  // statement may throw an exception, which we simply pass up to the caller
2186  const Point<dim> p_unit =
2187  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
2188  if (p_unit[0] == std::numeric_limits<double>::infinity())
2189  AssertThrow(false,
2191  return p_unit;
2192 }
2193 
2194 
2195 
2196 template <int dim, int spacedim>
2197 void
2199  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2200  const ArrayView<const Point<spacedim>> & real_points,
2201  const ArrayView<Point<dim>> & unit_points) const
2202 {
2203  AssertDimension(real_points.size(), unit_points.size());
2204  const std::vector<Point<spacedim>> support_points =
2205  this->compute_mapping_support_points(cell);
2206 
2207  // From the chosen (high-order) support points, now only pick the first
2208  // 2^dim points and construct an affine approximation from those.
2209  const std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
2210  affine_factors = GridTools::affine_cell_approximation<dim>(
2211  ArrayView<const Point<spacedim>>(support_points.data(),
2213  const DerivativeForm<1, spacedim, dim> A_inv =
2214  affine_factors.first.covariant_form().transpose();
2215 
2216  const unsigned int n_points = real_points.size();
2217  const unsigned int n_lanes = VectorizedArray<double>::size();
2218 
2219  // Use the more heavy VectorizedArray code path if there is more than
2220  // one point left to compute
2221  for (unsigned int i = 0; i < n_points; i += n_lanes)
2222  if (n_points - i > 1)
2223  {
2225  for (unsigned int j = 0; j < n_lanes; ++j)
2226  if (i + j < n_points)
2227  for (unsigned int d = 0; d < spacedim; ++d)
2228  p_vec[d][j] = real_points[i + j][d];
2229  else
2230  for (unsigned int d = 0; d < spacedim; ++d)
2231  p_vec[d][j] = real_points[i][d];
2232 
2233  // Compute an initial guess by inverting the affine approximation
2234  // A * x_unit + b = x_real
2236  for (unsigned int d = 0; d < spacedim; ++d)
2237  rhs[d] = affine_factors.second[d];
2238  rhs = p_vec - rhs;
2239 
2240  Point<dim, VectorizedArray<double>> initial_guess;
2241  for (unsigned int d = 0; d < dim; ++d)
2242  {
2243  initial_guess[d] = A_inv[d][0] * rhs[0];
2244  for (unsigned int e = 1; e < spacedim; ++e)
2245  initial_guess[d] += A_inv[d][e] * rhs[e];
2246  }
2248  internal::MappingQGenericImplementation::
2249  do_transform_real_to_unit_cell_internal<dim, spacedim>(
2250  p_vec,
2252  support_points,
2255 
2256  // If the vectorized computation failed, it could be that only some of
2257  // the lanes failed but others would have succeeded if we had let them
2258  // compute alone without interference (like negative Jacobian
2259  // determinants) from other SIMD lanes. Repeat the computation in this
2260  // unlikely case with scalar arguments.
2261  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
2262  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
2263  unit_points[i + j] = internal::MappingQGenericImplementation::
2264  do_transform_real_to_unit_cell_internal<dim, spacedim>(
2265  real_points[i + j],
2268  A_inv, real_points[i + j] - affine_factors.second))),
2269  support_points,
2272  else
2273  for (unsigned int d = 0; d < dim; ++d)
2274  unit_points[i + j][d] = unit_point[d][j];
2275  }
2276  else
2277  unit_points[i] = internal::MappingQGenericImplementation::
2278  do_transform_real_to_unit_cell_internal<dim, spacedim>(
2279  real_points[i],
2281  apply_transformation(A_inv,
2282  real_points[i] - affine_factors.second))),
2283  support_points,
2286 }
2287 
2288 
2289 
2290 template <int dim, int spacedim>
2293  const UpdateFlags in) const
2294 {
2295  // add flags if the respective quantities are necessary to compute
2296  // what we need. note that some flags appear in both the conditions
2297  // and in subsequent set operations. this leads to some circular
2298  // logic. the only way to treat this is to iterate. since there are
2299  // 5 if-clauses in the loop, it will take at most 5 iterations to
2300  // converge. do them:
2301  UpdateFlags out = in;
2302  for (unsigned int i = 0; i < 5; ++i)
2303  {
2304  // The following is a little incorrect:
2305  // If not applied on a face,
2306  // update_boundary_forms does not
2307  // make sense. On the other hand,
2308  // it is necessary on a
2309  // face. Currently,
2310  // update_boundary_forms is simply
2311  // ignored for the interior of a
2312  // cell.
2314  out |= update_boundary_forms;
2315 
2320 
2321  if (out &
2326 
2327  // The contravariant transformation is used in the Piola
2328  // transformation, which requires the determinant of the Jacobi
2329  // matrix of the transformation. Because we have no way of
2330  // knowing here whether the finite element wants to use the
2331  // contravariant or the Piola transforms, we add the JxW values
2332  // to the list of flags to be updated for each cell.
2334  out |= update_volume_elements;
2335 
2336  // the same is true when computing normal vectors: they require
2337  // the determinant of the Jacobian
2338  if (out & update_normal_vectors)
2339  out |= update_volume_elements;
2340  }
2341 
2342  return out;
2343 }
2344 
2345 
2346 
2347 template <int dim, int spacedim>
2348 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2350  const Quadrature<dim> &q) const
2351 {
2352  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2353  std::make_unique<InternalData>(polynomial_degree);
2354  auto &data = dynamic_cast<InternalData &>(*data_ptr);
2355  data.initialize(this->requires_update_flags(update_flags), q, q.size());
2356 
2357  return data_ptr;
2358 }
2359 
2360 
2361 
2362 template <int dim, int spacedim>
2363 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2365  const UpdateFlags update_flags,
2366  const Quadrature<dim - 1> &quadrature) const
2367 {
2368  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2369  std::make_unique<InternalData>(polynomial_degree);
2370  auto &data = dynamic_cast<InternalData &>(*data_ptr);
2371  data.initialize_face(this->requires_update_flags(update_flags),
2373  ReferenceCell::get_hypercube(dim), quadrature),
2374  quadrature.size());
2375 
2376  return data_ptr;
2377 }
2378 
2379 
2380 
2381 template <int dim, int spacedim>
2382 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2384  const UpdateFlags update_flags,
2385  const Quadrature<dim - 1> &quadrature) const
2386 {
2387  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2388  std::make_unique<InternalData>(polynomial_degree);
2389  auto &data = dynamic_cast<InternalData &>(*data_ptr);
2390  data.initialize_face(this->requires_update_flags(update_flags),
2392  ReferenceCell::get_hypercube(dim), quadrature),
2393  quadrature.size());
2394 
2395  return data_ptr;
2396 }
2397 
2398 
2399 
2400 template <int dim, int spacedim>
2403  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2404  const CellSimilarity::Similarity cell_similarity,
2405  const Quadrature<dim> & quadrature,
2406  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
2408  &output_data) const
2409 {
2410  // ensure that the following static_cast is really correct:
2411  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
2412  ExcInternalError());
2413  const InternalData &data = static_cast<const InternalData &>(internal_data);
2414 
2415  const unsigned int n_q_points = quadrature.size();
2416 
2417  // recompute the support points of the transformation of this
2418  // cell. we tried to be clever here in an earlier version of the
2419  // library by checking whether the cell is the same as the one we
2420  // had visited last, but it turns out to be difficult to determine
2421  // that because a cell for the purposes of a mapping is
2422  // characterized not just by its (triangulation, level, index)
2423  // triple, but also by the locations of its vertices, the manifold
2424  // object attached to the cell and all of its bounding faces/edges,
2425  // etc. to reliably test that the "cell" we are on is, therefore,
2426  // not easily done
2428  data.cell_of_current_support_points = cell;
2429 
2430  // if the order of the mapping is greater than 1, then do not reuse any cell
2431  // similarity information. This is necessary because the cell similarity
2432  // value is computed with just cell vertices and does not take into account
2433  // cell curvature.
2434  const CellSimilarity::Similarity computed_cell_similarity =
2435  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
2436 
2437  if (dim > 1 && data.tensor_product_quadrature)
2438  {
2439  internal::MappingQGenericImplementation::
2440  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2441  computed_cell_similarity,
2442  data,
2443  output_data.quadrature_points,
2444  output_data.jacobian_grads);
2445  }
2446  else
2447  {
2448  internal::MappingQGenericImplementation::maybe_compute_q_points<dim,
2449  spacedim>(
2451  data,
2452  output_data.quadrature_points);
2453 
2454  internal::MappingQGenericImplementation::maybe_update_Jacobians<dim,
2455  spacedim>(
2456  computed_cell_similarity,
2458  data);
2459 
2460  internal::MappingQGenericImplementation::maybe_update_jacobian_grads<
2461  dim,
2462  spacedim>(computed_cell_similarity,
2464  data,
2465  output_data.jacobian_grads);
2466  }
2467 
2468  internal::MappingQGenericImplementation::
2469  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2470  computed_cell_similarity,
2472  data,
2473  output_data.jacobian_pushed_forward_grads);
2474 
2475  internal::MappingQGenericImplementation::
2476  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2477  computed_cell_similarity,
2479  data,
2480  output_data.jacobian_2nd_derivatives);
2481 
2482  internal::MappingQGenericImplementation::
2483  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2484  computed_cell_similarity,
2486  data,
2488 
2489  internal::MappingQGenericImplementation::
2490  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2491  computed_cell_similarity,
2493  data,
2494  output_data.jacobian_3rd_derivatives);
2495 
2496  internal::MappingQGenericImplementation::
2497  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2498  computed_cell_similarity,
2500  data,
2502 
2503  const UpdateFlags update_flags = data.update_each;
2504  const std::vector<double> &weights = quadrature.get_weights();
2505 
2506  // Multiply quadrature weights by absolute value of Jacobian determinants or
2507  // the area element g=sqrt(DX^t DX) in case of codim > 0
2508 
2509  if (update_flags & (update_normal_vectors | update_JxW_values))
2510  {
2511  AssertDimension(output_data.JxW_values.size(), n_q_points);
2512 
2513  Assert(!(update_flags & update_normal_vectors) ||
2514  (output_data.normal_vectors.size() == n_q_points),
2515  ExcDimensionMismatch(output_data.normal_vectors.size(),
2516  n_q_points));
2517 
2518 
2519  if (computed_cell_similarity != CellSimilarity::translation)
2520  for (unsigned int point = 0; point < n_q_points; ++point)
2521  {
2522  if (dim == spacedim)
2523  {
2524  const double det = data.contravariant[point].determinant();
2525 
2526  // check for distorted cells.
2527 
2528  // TODO: this allows for anisotropies of up to 1e6 in 3D and
2529  // 1e12 in 2D. might want to find a finer
2530  // (dimension-independent) criterion
2531  Assert(det >
2532  1e-12 * Utilities::fixed_power<dim>(
2533  cell->diameter() / std::sqrt(double(dim))),
2535  cell->center(), det, point)));
2536 
2537  output_data.JxW_values[point] = weights[point] * det;
2538  }
2539  // if dim==spacedim, then there is no cell normal to
2540  // compute. since this is for FEValues (and not FEFaceValues),
2541  // there are also no face normals to compute
2542  else // codim>0 case
2543  {
2544  Tensor<1, spacedim> DX_t[dim];
2545  for (unsigned int i = 0; i < spacedim; ++i)
2546  for (unsigned int j = 0; j < dim; ++j)
2547  DX_t[j][i] = data.contravariant[point][i][j];
2548 
2549  Tensor<2, dim> G; // First fundamental form
2550  for (unsigned int i = 0; i < dim; ++i)
2551  for (unsigned int j = 0; j < dim; ++j)
2552  G[i][j] = DX_t[i] * DX_t[j];
2553 
2554  output_data.JxW_values[point] =
2555  std::sqrt(determinant(G)) * weights[point];
2556 
2557  if (computed_cell_similarity ==
2559  {
2560  // we only need to flip the normal
2561  if (update_flags & update_normal_vectors)
2562  output_data.normal_vectors[point] *= -1.;
2563  }
2564  else
2565  {
2566  if (update_flags & update_normal_vectors)
2567  {
2568  Assert(spacedim == dim + 1,
2569  ExcMessage(
2570  "There is no (unique) cell normal for " +
2572  "-dimensional cells in " +
2573  Utilities::int_to_string(spacedim) +
2574  "-dimensional space. This only works if the "
2575  "space dimension is one greater than the "
2576  "dimensionality of the mesh cells."));
2577 
2578  if (dim == 1)
2579  output_data.normal_vectors[point] =
2580  cross_product_2d(-DX_t[0]);
2581  else // dim == 2
2582  output_data.normal_vectors[point] =
2583  cross_product_3d(DX_t[0], DX_t[1]);
2584 
2585  output_data.normal_vectors[point] /=
2586  output_data.normal_vectors[point].norm();
2587 
2588  if (cell->direction_flag() == false)
2589  output_data.normal_vectors[point] *= -1.;
2590  }
2591  }
2592  } // codim>0 case
2593  }
2594  }
2595 
2596 
2597 
2598  // copy values from InternalData to vector given by reference
2599  if (update_flags & update_jacobians)
2600  {
2601  AssertDimension(output_data.jacobians.size(), n_q_points);
2602  if (computed_cell_similarity != CellSimilarity::translation)
2603  for (unsigned int point = 0; point < n_q_points; ++point)
2604  output_data.jacobians[point] = data.contravariant[point];
2605  }
2606 
2607  // copy values from InternalData to vector given by reference
2608  if (update_flags & update_inverse_jacobians)
2609  {
2610  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
2611  if (computed_cell_similarity != CellSimilarity::translation)
2612  for (unsigned int point = 0; point < n_q_points; ++point)
2613  output_data.inverse_jacobians[point] =
2614  data.covariant[point].transpose();
2615  }
2616 
2617  return computed_cell_similarity;
2618 }
2619 
2620 
2621 
2622 namespace internal
2623 {
2624  namespace MappingQGenericImplementation
2625  {
2626  namespace
2627  {
2637  template <int dim, int spacedim>
2638  void
2639  maybe_compute_face_data(
2640  const ::MappingQGeneric<dim, spacedim> &mapping,
2641  const typename ::Triangulation<dim, spacedim>::cell_iterator
2642  & cell,
2643  const unsigned int face_no,
2644  const unsigned int subface_no,
2645  const unsigned int n_q_points,
2646  const std::vector<double> &weights,
2647  const typename ::MappingQGeneric<dim, spacedim>::InternalData
2648  &data,
2650  &output_data)
2651  {
2652  const UpdateFlags update_flags = data.update_each;
2653 
2654  if (update_flags &
2657  {
2658  if (update_flags & update_boundary_forms)
2659  AssertDimension(output_data.boundary_forms.size(), n_q_points);
2660  if (update_flags & update_normal_vectors)
2661  AssertDimension(output_data.normal_vectors.size(), n_q_points);
2662  if (update_flags & update_JxW_values)
2663  AssertDimension(output_data.JxW_values.size(), n_q_points);
2664 
2665  Assert(data.aux.size() + 1 >= dim, ExcInternalError());
2666 
2667  // first compute some common data that is used for evaluating
2668  // all of the flags below
2669 
2670  // map the unit tangentials to the real cell. checking for d!=dim-1
2671  // eliminates compiler warnings regarding unsigned int expressions <
2672  // 0.
2673  for (unsigned int d = 0; d != dim - 1; ++d)
2674  {
2676  data.unit_tangentials.size(),
2677  ExcInternalError());
2678  Assert(
2679  data.aux[d].size() <=
2680  data
2681  .unit_tangentials[face_no +
2683  .size(),
2684  ExcInternalError());
2685 
2686  mapping.transform(
2688  data
2689  .unit_tangentials[face_no +
2692  data,
2693  make_array_view(data.aux[d]));
2694  }
2695 
2696  if (update_flags & update_boundary_forms)
2697  {
2698  // if dim==spacedim, we can use the unit tangentials to compute
2699  // the boundary form by simply taking the cross product
2700  if (dim == spacedim)
2701  {
2702  for (unsigned int i = 0; i < n_q_points; ++i)
2703  switch (dim)
2704  {
2705  case 1:
2706  // in 1d, we don't have access to any of the
2707  // data.aux fields (because it has only dim-1
2708  // components), but we can still compute the
2709  // boundary form by simply looking at the number of
2710  // the face
2711  output_data.boundary_forms[i][0] =
2712  (face_no == 0 ? -1 : +1);
2713  break;
2714  case 2:
2715  output_data.boundary_forms[i] =
2716  cross_product_2d(data.aux[0][i]);
2717  break;
2718  case 3:
2719  output_data.boundary_forms[i] =
2720  cross_product_3d(data.aux[0][i], data.aux[1][i]);
2721  break;
2722  default:
2723  Assert(false, ExcNotImplemented());
2724  }
2725  }
2726  else //(dim < spacedim)
2727  {
2728  // in the codim-one case, the boundary form results from the
2729  // cross product of all the face tangential vectors and the
2730  // cell normal vector
2731  //
2732  // to compute the cell normal, use the same method used in
2733  // fill_fe_values for cells above
2734  AssertDimension(data.contravariant.size(), n_q_points);
2735 
2736  for (unsigned int point = 0; point < n_q_points; ++point)
2737  {
2738  if (dim == 1)
2739  {
2740  // J is a tangent vector
2741  output_data.boundary_forms[point] =
2742  data.contravariant[point].transpose()[0];
2743  output_data.boundary_forms[point] /=
2744  (face_no == 0 ? -1. : +1.) *
2745  output_data.boundary_forms[point].norm();
2746  }
2747 
2748  if (dim == 2)
2749  {
2751  data.contravariant[point].transpose();
2752 
2753  Tensor<1, spacedim> cell_normal =
2754  cross_product_3d(DX_t[0], DX_t[1]);
2755  cell_normal /= cell_normal.norm();
2756 
2757  // then compute the face normal from the face
2758  // tangent and the cell normal:
2759  output_data.boundary_forms[point] =
2760  cross_product_3d(data.aux[0][point], cell_normal);
2761  }
2762  }
2763  }
2764  }
2765 
2766  if (update_flags & update_JxW_values)
2767  for (unsigned int i = 0; i < output_data.boundary_forms.size();
2768  ++i)
2769  {
2770  output_data.JxW_values[i] =
2771  output_data.boundary_forms[i].norm() * weights[i];
2772 
2773  if (subface_no != numbers::invalid_unsigned_int)
2774  {
2775  const double area_ratio =
2777  cell->subface_case(face_no), subface_no);
2778  output_data.JxW_values[i] *= area_ratio;
2779  }
2780  }
2781 
2782  if (update_flags & update_normal_vectors)
2783  for (unsigned int i = 0; i < output_data.normal_vectors.size();
2784  ++i)
2785  output_data.normal_vectors[i] =
2786  Point<spacedim>(output_data.boundary_forms[i] /
2787  output_data.boundary_forms[i].norm());
2788 
2789  if (update_flags & update_jacobians)
2790  for (unsigned int point = 0; point < n_q_points; ++point)
2791  output_data.jacobians[point] = data.contravariant[point];
2792 
2793  if (update_flags & update_inverse_jacobians)
2794  for (unsigned int point = 0; point < n_q_points; ++point)
2795  output_data.inverse_jacobians[point] =
2796  data.covariant[point].transpose();
2797  }
2798  }
2799 
2800 
2807  template <int dim, int spacedim>
2808  void
2809  do_fill_fe_face_values(
2810  const ::MappingQGeneric<dim, spacedim> &mapping,
2811  const typename ::Triangulation<dim, spacedim>::cell_iterator
2812  & cell,
2813  const unsigned int face_no,
2814  const unsigned int subface_no,
2815  const typename QProjector<dim>::DataSetDescriptor data_set,
2816  const Quadrature<dim - 1> & quadrature,
2817  const typename ::MappingQGeneric<dim, spacedim>::InternalData
2818  &data,
2820  &output_data)
2821  {
2822  if (dim > 1 && data.tensor_product_quadrature)
2823  {
2824  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2826  data,
2827  output_data.quadrature_points,
2828  output_data.jacobian_grads);
2829  }
2830  else
2831  {
2832  maybe_compute_q_points<dim, spacedim>(
2833  data_set, data, output_data.quadrature_points);
2834  maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
2835  data_set,
2836  data);
2837  maybe_update_jacobian_grads<dim, spacedim>(
2838  CellSimilarity::none, data_set, data, output_data.jacobian_grads);
2839  }
2840  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2842  data_set,
2843  data,
2844  output_data.jacobian_pushed_forward_grads);
2845  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2847  data_set,
2848  data,
2849  output_data.jacobian_2nd_derivatives);
2850  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2852  data_set,
2853  data,
2855  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2857  data_set,
2858  data,
2859  output_data.jacobian_3rd_derivatives);
2860  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2862  data_set,
2863  data,
2865 
2866  maybe_compute_face_data(mapping,
2867  cell,
2868  face_no,
2869  subface_no,
2870  quadrature.size(),
2871  quadrature.get_weights(),
2872  data,
2873  output_data);
2874  }
2875  } // namespace
2876  } // namespace MappingQGenericImplementation
2877 } // namespace internal
2878 
2879 
2880 
2881 template <int dim, int spacedim>
2882 void
2884  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2885  const unsigned int face_no,
2886  const Quadrature<dim - 1> & quadrature,
2887  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
2889  &output_data) const
2890 {
2891  // ensure that the following cast is really correct:
2892  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
2893  ExcInternalError());
2894  const InternalData &data = static_cast<const InternalData &>(internal_data);
2895 
2896  // if necessary, recompute the support points of the transformation of this
2897  // cell (note that we need to first check the triangulation pointer, since
2898  // otherwise the second test might trigger an exception if the triangulations
2899  // are not the same)
2900  if ((data.mapping_support_points.size() == 0) ||
2901  (&cell->get_triangulation() !=
2902  &data.cell_of_current_support_points->get_triangulation()) ||
2903  (cell != data.cell_of_current_support_points))
2904  {
2906  data.cell_of_current_support_points = cell;
2907  }
2908 
2909  internal::MappingQGenericImplementation::do_fill_fe_face_values(
2910  *this,
2911  cell,
2912  face_no,
2915  face_no,
2916  cell->face_orientation(face_no),
2917  cell->face_flip(face_no),
2918  cell->face_rotation(face_no),
2919  quadrature.size()),
2920  quadrature,
2921  data,
2922  output_data);
2923 }
2924 
2925 
2926 
2927 template <int dim, int spacedim>
2928 void
2930  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2931  const unsigned int face_no,
2932  const unsigned int subface_no,
2933  const Quadrature<dim - 1> & quadrature,
2934  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
2936  &output_data) const
2937 {
2938  // ensure that the following cast is really correct:
2939  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
2940  ExcInternalError());
2941  const InternalData &data = static_cast<const InternalData &>(internal_data);
2942 
2943  // if necessary, recompute the support points of the transformation of this
2944  // cell (note that we need to first check the triangulation pointer, since
2945  // otherwise the second test might trigger an exception if the triangulations
2946  // are not the same)
2947  if ((data.mapping_support_points.size() == 0) ||
2948  (&cell->get_triangulation() !=
2949  &data.cell_of_current_support_points->get_triangulation()) ||
2950  (cell != data.cell_of_current_support_points))
2951  {
2953  data.cell_of_current_support_points = cell;
2954  }
2955 
2956  internal::MappingQGenericImplementation::do_fill_fe_face_values(
2957  *this,
2958  cell,
2959  face_no,
2960  subface_no,
2962  dim),
2963  face_no,
2964  subface_no,
2965  cell->face_orientation(face_no),
2966  cell->face_flip(face_no),
2967  cell->face_rotation(face_no),
2968  quadrature.size(),
2969  cell->subface_case(face_no)),
2970  quadrature,
2971  data,
2972  output_data);
2973 }
2974 
2975 
2976 
2977 namespace internal
2978 {
2979  namespace MappingQGenericImplementation
2980  {
2981  namespace
2982  {
2983  template <int dim, int spacedim, int rank>
2984  void
2985  transform_fields(
2986  const ArrayView<const Tensor<rank, dim>> & input,
2987  const MappingKind mapping_kind,
2988  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2989  const ArrayView<Tensor<rank, spacedim>> & output)
2990  {
2991  AssertDimension(input.size(), output.size());
2992  Assert((dynamic_cast<const typename ::
2993  MappingQGeneric<dim, spacedim>::InternalData *>(
2994  &mapping_data) != nullptr),
2995  ExcInternalError());
2996  const typename ::MappingQGeneric<dim, spacedim>::InternalData
2997  &data =
2998  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
2999  InternalData &>(mapping_data);
3000 
3001  switch (mapping_kind)
3002  {
3003  case mapping_contravariant:
3004  {
3005  Assert(
3006  data.update_each & update_contravariant_transformation,
3008  "update_contravariant_transformation"));
3009 
3010  for (unsigned int i = 0; i < output.size(); ++i)
3011  output[i] =
3012  apply_transformation(data.contravariant[i], input[i]);
3013 
3014  return;
3015  }
3016 
3017  case mapping_piola:
3018  {
3019  Assert(
3020  data.update_each & update_contravariant_transformation,
3022  "update_contravariant_transformation"));
3023  Assert(
3024  data.update_each & update_volume_elements,
3026  "update_volume_elements"));
3027  Assert(rank == 1, ExcMessage("Only for rank 1"));
3028  if (rank != 1)
3029  return;
3030 
3031  for (unsigned int i = 0; i < output.size(); ++i)
3032  {
3033  output[i] =
3034  apply_transformation(data.contravariant[i], input[i]);
3035  output[i] /= data.volume_elements[i];
3036  }
3037  return;
3038  }
3039  // We still allow this operation as in the
3040  // reference cell Derivatives are Tensor
3041  // rather than DerivativeForm
3042  case mapping_covariant:
3043  {
3044  Assert(
3045  data.update_each & update_contravariant_transformation,
3047  "update_covariant_transformation"));
3048 
3049  for (unsigned int i = 0; i < output.size(); ++i)
3050  output[i] = apply_transformation(data.covariant[i], input[i]);
3051 
3052  return;
3053  }
3054 
3055  default:
3056  Assert(false, ExcNotImplemented());
3057  }
3058  }
3059 
3060 
3061  template <int dim, int spacedim, int rank>
3062  void
3063  transform_gradients(
3064  const ArrayView<const Tensor<rank, dim>> & input,
3065  const MappingKind mapping_kind,
3066  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3067  const ArrayView<Tensor<rank, spacedim>> & output)
3068  {
3069  AssertDimension(input.size(), output.size());
3070  Assert((dynamic_cast<const typename ::
3071  MappingQGeneric<dim, spacedim>::InternalData *>(
3072  &mapping_data) != nullptr),
3073  ExcInternalError());
3074  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3075  &data =
3076  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3077  InternalData &>(mapping_data);
3078 
3079  switch (mapping_kind)
3080  {
3082  {
3083  Assert(
3084  data.update_each & update_covariant_transformation,
3086  "update_covariant_transformation"));
3087  Assert(
3088  data.update_each & update_contravariant_transformation,
3090  "update_contravariant_transformation"));
3091  Assert(rank == 2, ExcMessage("Only for rank 2"));
3092 
3093  for (unsigned int i = 0; i < output.size(); ++i)
3094  {
3096  apply_transformation(data.contravariant[i],
3097  transpose(input[i]));
3098  output[i] =
3099  apply_transformation(data.covariant[i], A.transpose());
3100  }
3101 
3102  return;
3103  }
3104 
3106  {
3107  Assert(
3108  data.update_each & update_covariant_transformation,
3110  "update_covariant_transformation"));
3111  Assert(rank == 2, ExcMessage("Only for rank 2"));
3112 
3113  for (unsigned int i = 0; i < output.size(); ++i)
3114  {
3116  apply_transformation(data.covariant[i],
3117  transpose(input[i]));
3118  output[i] =
3119  apply_transformation(data.covariant[i], A.transpose());
3120  }
3121 
3122  return;
3123  }
3124 
3126  {
3127  Assert(
3128  data.update_each & update_covariant_transformation,
3130  "update_covariant_transformation"));
3131  Assert(
3132  data.update_each & update_contravariant_transformation,
3134  "update_contravariant_transformation"));
3135  Assert(
3136  data.update_each & update_volume_elements,
3138  "update_volume_elements"));
3139  Assert(rank == 2, ExcMessage("Only for rank 2"));
3140 
3141  for (unsigned int i = 0; i < output.size(); ++i)
3142  {
3144  apply_transformation(data.covariant[i], input[i]);
3145  const Tensor<2, spacedim> T =
3146  apply_transformation(data.contravariant[i],
3147  A.transpose());
3148 
3149  output[i] = transpose(T);
3150  output[i] /= data.volume_elements[i];
3151  }
3152 
3153  return;
3154  }
3155 
3156  default:
3157  Assert(false, ExcNotImplemented());
3158  }
3159  }
3160 
3161 
3162 
3163  template <int dim, int spacedim>
3164  void
3165  transform_hessians(
3166  const ArrayView<const Tensor<3, dim>> & input,
3167  const MappingKind mapping_kind,
3168  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3169  const ArrayView<Tensor<3, spacedim>> & output)
3170  {
3171  AssertDimension(input.size(), output.size());
3172  Assert((dynamic_cast<const typename ::
3173  MappingQGeneric<dim, spacedim>::InternalData *>(
3174  &mapping_data) != nullptr),
3175  ExcInternalError());
3176  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3177  &data =
3178  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3179  InternalData &>(mapping_data);
3180 
3181  switch (mapping_kind)
3182  {
3184  {
3185  Assert(
3186  data.update_each & update_covariant_transformation,
3188  "update_covariant_transformation"));
3189  Assert(
3190  data.update_each & update_contravariant_transformation,
3192  "update_contravariant_transformation"));
3193 
3194  for (unsigned int q = 0; q < output.size(); ++q)
3195  for (unsigned int i = 0; i < spacedim; ++i)
3196  {
3197  double tmp1[dim][dim];
3198  for (unsigned int J = 0; J < dim; ++J)
3199  for (unsigned int K = 0; K < dim; ++K)
3200  {
3201  tmp1[J][K] =
3202  data.contravariant[q][i][0] * input[q][0][J][K];
3203  for (unsigned int I = 1; I < dim; ++I)
3204  tmp1[J][K] +=
3205  data.contravariant[q][i][I] * input[q][I][J][K];
3206  }
3207  for (unsigned int j = 0; j < spacedim; ++j)
3208  {
3209  double tmp2[dim];
3210  for (unsigned int K = 0; K < dim; ++K)
3211  {
3212  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3213  for (unsigned int J = 1; J < dim; ++J)
3214  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3215  }
3216  for (unsigned int k = 0; k < spacedim; ++k)
3217  {
3218  output[q][i][j][k] =
3219  data.covariant[q][k][0] * tmp2[0];
3220  for (unsigned int K = 1; K < dim; ++K)
3221  output[q][i][j][k] +=
3222  data.covariant[q][k][K] * tmp2[K];
3223  }
3224  }
3225  }
3226  return;
3227  }
3228 
3230  {
3231  Assert(
3232  data.update_each & update_covariant_transformation,
3234  "update_covariant_transformation"));
3235 
3236  for (unsigned int q = 0; q < output.size(); ++q)
3237  for (unsigned int i = 0; i < spacedim; ++i)
3238  {
3239  double tmp1[dim][dim];
3240  for (unsigned int J = 0; J < dim; ++J)
3241  for (unsigned int K = 0; K < dim; ++K)
3242  {
3243  tmp1[J][K] =
3244  data.covariant[q][i][0] * input[q][0][J][K];
3245  for (unsigned int I = 1; I < dim; ++I)
3246  tmp1[J][K] +=
3247  data.covariant[q][i][I] * input[q][I][J][K];
3248  }
3249  for (unsigned int j = 0; j < spacedim; ++j)
3250  {
3251  double tmp2[dim];
3252  for (unsigned int K = 0; K < dim; ++K)
3253  {
3254  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3255  for (unsigned int J = 1; J < dim; ++J)
3256  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3257  }
3258  for (unsigned int k = 0; k < spacedim; ++k)
3259  {
3260  output[q][i][j][k] =
3261  data.covariant[q][k][0] * tmp2[0];
3262  for (unsigned int K = 1; K < dim; ++K)
3263  output[q][i][j][k] +=
3264  data.covariant[q][k][K] * tmp2[K];
3265  }
3266  }
3267  }
3268 
3269  return;
3270  }
3271 
3272  case mapping_piola_hessian:
3273  {
3274  Assert(
3275  data.update_each & update_covariant_transformation,
3277  "update_covariant_transformation"));
3278  Assert(
3279  data.update_each & update_contravariant_transformation,
3281  "update_contravariant_transformation"));
3282  Assert(
3283  data.update_each & update_volume_elements,
3285  "update_volume_elements"));
3286 
3287  for (unsigned int q = 0; q < output.size(); ++q)
3288  for (unsigned int i = 0; i < spacedim; ++i)
3289  {
3290  double factor[dim];
3291  for (unsigned int I = 0; I < dim; ++I)
3292  factor[I] =
3293  data.contravariant[q][i][I] / data.volume_elements[q];
3294  double tmp1[dim][dim];
3295  for (unsigned int J = 0; J < dim; ++J)
3296  for (unsigned int K = 0; K < dim; ++K)
3297  {
3298  tmp1[J][K] = factor[0] * input[q][0][J][K];
3299  for (unsigned int I = 1; I < dim; ++I)
3300  tmp1[J][K] += factor[I] * input[q][I][J][K];
3301  }
3302  for (unsigned int j = 0; j < spacedim; ++j)
3303  {
3304  double tmp2[dim];
3305  for (unsigned int K = 0; K < dim; ++K)
3306  {
3307  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3308  for (unsigned int J = 1; J < dim; ++J)
3309  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3310  }
3311  for (unsigned int k = 0; k < spacedim; ++k)
3312  {
3313  output[q][i][j][k] =
3314  data.covariant[q][k][0] * tmp2[0];
3315  for (unsigned int K = 1; K < dim; ++K)
3316  output[q][i][j][k] +=
3317  data.covariant[q][k][K] * tmp2[K];
3318  }
3319  }
3320  }
3321 
3322  return;
3323  }
3324 
3325  default:
3326  Assert(false, ExcNotImplemented());
3327  }
3328  }
3329 
3330 
3331 
3332  template <int dim, int spacedim, int rank>
3333  void
3334  transform_differential_forms(
3335  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
3336  const MappingKind mapping_kind,
3337  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3338  const ArrayView<Tensor<rank + 1, spacedim>> & output)
3339  {
3340  AssertDimension(input.size(), output.size());
3341  Assert((dynamic_cast<const typename ::
3342  MappingQGeneric<dim, spacedim>::InternalData *>(
3343  &mapping_data) != nullptr),
3344  ExcInternalError());
3345  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3346  &data =
3347  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3348  InternalData &>(mapping_data);
3349 
3350  switch (mapping_kind)
3351  {
3352  case mapping_covariant:
3353  {
3354  Assert(
3355  data.update_each & update_contravariant_transformation,
3357  "update_covariant_transformation"));
3358 
3359  for (unsigned int i = 0; i < output.size(); ++i)
3360  output[i] = apply_transformation(data.covariant[i], input[i]);
3361 
3362  return;
3363  }
3364  default:
3365  Assert(false, ExcNotImplemented());
3366  }
3367  }
3368  } // namespace
3369  } // namespace MappingQGenericImplementation
3370 } // namespace internal
3371 
3372 
3373 
3374 template <int dim, int spacedim>
3375 void
3377  const ArrayView<const Tensor<1, dim>> & input,
3378  const MappingKind mapping_kind,
3379  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3380  const ArrayView<Tensor<1, spacedim>> & output) const
3381 {
3382  internal::MappingQGenericImplementation::transform_fields(input,
3383  mapping_kind,
3384  mapping_data,
3385  output);
3386 }
3387 
3388 
3389 
3390 template <int dim, int spacedim>
3391 void
3393  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
3394  const MappingKind mapping_kind,
3395  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3396  const ArrayView<Tensor<2, spacedim>> & output) const
3397 {
3398  internal::MappingQGenericImplementation::transform_differential_forms(
3399  input, mapping_kind, mapping_data, output);
3400 }
3401 
3402 
3403 
3404 template <int dim, int spacedim>
3405 void
3407  const ArrayView<const Tensor<2, dim>> & input,
3408  const MappingKind mapping_kind,
3409  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3410  const ArrayView<Tensor<2, spacedim>> & output) const
3411 {
3412  switch (mapping_kind)
3413  {
3414  case mapping_contravariant:
3415  internal::MappingQGenericImplementation::transform_fields(input,
3416  mapping_kind,
3417  mapping_data,
3418  output);
3419  return;
3420 
3424  internal::MappingQGenericImplementation::transform_gradients(
3425  input, mapping_kind, mapping_data, output);
3426  return;
3427  default:
3428  Assert(false, ExcNotImplemented());
3429  }
3430 }
3431 
3432 
3433 
3434 template <int dim, int spacedim>
3435 void
3437  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
3438  const MappingKind mapping_kind,
3439  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3440  const ArrayView<Tensor<3, spacedim>> & output) const
3441 {
3442  AssertDimension(input.size(), output.size());
3443  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
3444  ExcInternalError());
3445  const InternalData &data = static_cast<const InternalData &>(mapping_data);
3446 
3447  switch (mapping_kind)
3448  {
3450  {
3453  "update_covariant_transformation"));
3454 
3455  for (unsigned int q = 0; q < output.size(); ++q)
3456  for (unsigned int i = 0; i < spacedim; ++i)
3457  for (unsigned int j = 0; j < spacedim; ++j)
3458  {
3459  double tmp[dim];
3460  for (unsigned int K = 0; K < dim; ++K)
3461  {
3462  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
3463  for (unsigned int J = 1; J < dim; ++J)
3464  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
3465  }
3466  for (unsigned int k = 0; k < spacedim; ++k)
3467  {
3468  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
3469  for (unsigned int K = 1; K < dim; ++K)
3470  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
3471  }
3472  }
3473  return;
3474  }
3475 
3476  default:
3477  Assert(false, ExcNotImplemented());
3478  }
3479 }
3480 
3481 
3482 
3483 template <int dim, int spacedim>
3484 void
3486  const ArrayView<const Tensor<3, dim>> & input,
3487  const MappingKind mapping_kind,
3488  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3489  const ArrayView<Tensor<3, spacedim>> & output) const
3490 {
3491  switch (mapping_kind)
3492  {
3493  case mapping_piola_hessian:
3496  internal::MappingQGenericImplementation::transform_hessians(
3497  input, mapping_kind, mapping_data, output);
3498  return;
3499  default:
3500  Assert(false, ExcNotImplemented());
3501  }
3502 }
3503 
3504 
3505 
3506 template <int dim, int spacedim>
3507 void
3509  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
3510  std::vector<Point<spacedim>> & a) const
3511 {
3512  // if we only need the midpoint, then ask for it.
3513  if (this->polynomial_degree == 2)
3514  {
3515  for (unsigned int line_no = 0;
3516  line_no < GeometryInfo<dim>::lines_per_cell;
3517  ++line_no)
3518  {
3519  const typename Triangulation<dim, spacedim>::line_iterator line =
3520  (dim == 1 ?
3521  static_cast<
3523  cell->line(line_no));
3524 
3525  const Manifold<dim, spacedim> &manifold =
3526  ((line->manifold_id() == numbers::flat_manifold_id) &&
3527  (dim < spacedim) ?
3528  cell->get_manifold() :
3529  line->get_manifold());
3530  a.push_back(manifold.get_new_point_on_line(line));
3531  }
3532  }
3533  else
3534  // otherwise call the more complicated functions and ask for inner points
3535  // from the manifold description
3536  {
3537  std::vector<Point<spacedim>> tmp_points;
3538  for (unsigned int line_no = 0;
3539  line_no < GeometryInfo<dim>::lines_per_cell;
3540  ++line_no)
3541  {
3542  const typename Triangulation<dim, spacedim>::line_iterator line =
3543  (dim == 1 ?
3544  static_cast<
3546  cell->line(line_no));
3547 
3548  const Manifold<dim, spacedim> &manifold =
3549  ((line->manifold_id() == numbers::flat_manifold_id) &&
3550  (dim < spacedim) ?
3551  cell->get_manifold() :
3552  line->get_manifold());
3553 
3554  const std::array<Point<spacedim>, 2> vertices{
3555  {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
3556  cell->vertex(
3558 
3559  const std::size_t n_rows =
3561  a.resize(a.size() + n_rows);
3562  auto a_view = make_array_view(a.end() - n_rows, a.end());
3563  manifold.get_new_points(
3564  make_array_view(vertices.begin(), vertices.end()),
3566  a_view);
3567  }
3568  }
3569 }
3570 
3571 
3572 
3573 template <>
3574 void
3577  std::vector<Point<3>> & a) const
3578 {
3579  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
3580 
3581  // used if face quad at boundary or entirely in the interior of the domain
3582  std::vector<Point<3>> tmp_points;
3583 
3584  // loop over all faces and collect points on them
3585  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
3586  {
3587  const Triangulation<3>::face_iterator face = cell->face(face_no);
3588 
3589 #ifdef DEBUG
3590  const bool face_orientation = cell->face_orientation(face_no),
3591  face_flip = cell->face_flip(face_no),
3592  face_rotation = cell->face_rotation(face_no);
3593  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
3594  lines_per_face = GeometryInfo<3>::lines_per_face;
3595 
3596  // some sanity checks up front
3597  for (unsigned int i = 0; i < vertices_per_face; ++i)
3598  Assert(face->vertex_index(i) ==
3599  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
3600  face_no, i, face_orientation, face_flip, face_rotation)),
3601  ExcInternalError());
3602 
3603  // indices of the lines that bound a face are given by GeometryInfo<3>::
3604  // face_to_cell_lines
3605  for (unsigned int i = 0; i < lines_per_face; ++i)
3606  Assert(face->line(i) ==
3608  face_no, i, face_orientation, face_flip, face_rotation)),
3609  ExcInternalError());
3610 #endif
3611  // extract the points surrounding a quad from the points
3612  // already computed. First get the 4 vertices and then the points on
3613  // the four lines
3614  boost::container::small_vector<Point<3>, 200> tmp_points(
3617  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
3618  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
3619  if (polynomial_degree > 1)
3620  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
3621  ++line)
3622  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
3623  tmp_points[4 + line * (polynomial_degree - 1) + i] =
3625  (polynomial_degree - 1) *
3626  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
3627  i];
3628 
3629  const std::size_t n_rows =
3631  a.resize(a.size() + n_rows);
3632  auto a_view = make_array_view(a.end() - n_rows, a.end());
3633  face->get_manifold().get_new_points(
3634  make_array_view(tmp_points.begin(), tmp_points.end()),
3636  a_view);
3637  }
3638 }
3639 
3640 
3641 
3642 template <>
3643 void
3646  std::vector<Point<3>> & a) const
3647 {
3648  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
3649  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
3650  vertices[i] = cell->vertex(i);
3651 
3652  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
3654  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
3655  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
3656  {
3657  Point<2> point(line_support_points[q1 + 1][0],
3658  line_support_points[q2 + 1][0]);
3659  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
3660  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
3661  }
3662 
3663  const std::size_t n_rows = weights.size(0);
3664  a.resize(a.size() + n_rows);
3665  auto a_view = make_array_view(a.end() - n_rows, a.end());
3666  cell->get_manifold().get_new_points(
3667  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
3668 }
3669 
3670 
3671 
3672 template <int dim, int spacedim>
3673 void
3676  std::vector<Point<spacedim>> &) const
3677 {
3678  Assert(false, ExcInternalError());
3679 }
3680 
3681 
3682 
3683 template <int dim, int spacedim>
3684 std::vector<Point<spacedim>>
3686  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
3687 {
3688  // get the vertices first
3689  std::vector<Point<spacedim>> a;
3690  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
3691  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
3692  a.push_back(cell->vertex(i));
3693 
3694  if (this->polynomial_degree > 1)
3695  {
3696  // check if all entities have the same manifold id which is when we can
3697  // simply ask the manifold for all points. the transfinite manifold can
3698  // do the interpolation better than this class, so if we detect that we
3699  // do not have to change anything here
3700  Assert(dim <= 3, ExcImpossibleInDim(dim));
3701  bool all_manifold_ids_are_equal = (dim == spacedim);
3702  if (all_manifold_ids_are_equal &&
3704  &cell->get_manifold()) == nullptr)
3705  {
3706  for (auto f : GeometryInfo<dim>::face_indices())
3707  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
3708  all_manifold_ids_are_equal = false;
3709 
3710  if (dim == 3)
3711  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3712  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
3713  all_manifold_ids_are_equal = false;
3714  }
3715 
3716  if (all_manifold_ids_are_equal)
3717  {
3718  const std::size_t n_rows = support_point_weights_cell.size(0);
3719  a.resize(a.size() + n_rows);
3720  auto a_view = make_array_view(a.end() - n_rows, a.end());
3721  cell->get_manifold().get_new_points(make_array_view(a.begin(),
3722  a.end() - n_rows),
3724  a_view);
3725  }
3726  else
3727  switch (dim)
3728  {
3729  case 1:
3730  add_line_support_points(cell, a);
3731  break;
3732  case 2:
3733  // in 2d, add the points on the four bounding lines to the
3734  // exterior (outer) points
3735  add_line_support_points(cell, a);
3736 
3737  // then get the interior support points
3738  if (dim != spacedim)
3739  add_quad_support_points(cell, a);
3740  else
3741  {
3742  const std::size_t n_rows =
3744  a.resize(a.size() + n_rows);
3745  auto a_view = make_array_view(a.end() - n_rows, a.end());
3746  cell->get_manifold().get_new_points(
3747  make_array_view(a.begin(), a.end() - n_rows),
3749  a_view);
3750  }
3751  break;
3752 
3753  case 3:
3754  // in 3d also add the points located on the boundary faces
3755  add_line_support_points(cell, a);
3756  add_quad_support_points(cell, a);
3757 
3758  // then compute the interior points
3759  {
3760  const std::size_t n_rows =
3762  a.resize(a.size() + n_rows);
3763  auto a_view = make_array_view(a.end() - n_rows, a.end());
3764  cell->get_manifold().get_new_points(
3765  make_array_view(a.begin(), a.end() - n_rows),
3767  a_view);
3768  }
3769  break;
3770 
3771  default:
3772  Assert(false, ExcNotImplemented());
3773  break;
3774  }
3775  }
3776 
3777  return a;
3778 }
3779 
3780 
3781 
3782 //--------------------------- Explicit instantiations -----------------------
3783 #include "mapping_q_generic.inst"
3784 
3785 
Transformed quadrature weights.
std::vector< Tensor< 2, dim > > shape_second_derivatives
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 > &)
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
Tensor< 1, spacedim, Number > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number > &grad_F, const Tensor< 1, dim, Number > &d_x)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
const unsigned int polynomial_degree
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1438
Contravariant transformation.
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
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=StaticMappingQ1< dim, spacedim >::mapping, const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:444
Table< 2, double > support_point_weights_cell
const std::vector< Point< dim > > & get_points() const
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
const std::vector< double > & get_weights() const
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Type get_hypercube(const unsigned int dim)
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
Definition: fe_dgq.h:110
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping.cc:28
std::vector< unsigned int > renumber_lexicographic_to_hierarchic
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
constexpr void clear()
const Point< dim > & point(const unsigned int i) const
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Determinant of the Jacobian.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:323
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
MappingQGeneric(const unsigned int polynomial_degree)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
Point< 2 > second
Definition: grid_out.cc:4336
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:564
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
InternalData(const unsigned int polynomial_degree)
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:2441
const unsigned int polynomial_degree
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:752
QGaussLobatto< 1 > line_support_points
T fixed_power(const T t)
Definition: utilities.h:1045
std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:316
std::vector< Polynomials::Polynomial< double > > polynomials_1d
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual void transform_points_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< spacedim >> &real_points, const ArrayView< Point< dim >> &unit_points) const override
static const char T
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
std::vector< std::vector< Tensor< 1, spacedim > > > aux
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1423
UpdateFlags
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const override
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
Abstract base class for mapping classes.
Definition: mapping.h:301
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:665
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
std::vector< Point< spacedim > > mapping_support_points
std::vector< Tensor< 3, dim > > shape_third_derivatives
VectorType::value_type * end(VectorType &V)
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose() const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Expression fabs(const Expression &x)
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< double > volume_elements
Gradient of volume element.
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int n_shape_functions
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
std::vector< Tensor< 1, dim > > shape_derivatives
unsigned int size() const
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4335
Point< dim > transform_real_to_unit_cell_internal(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit) const
std::vector< Point< spacedim > > quadrature_points
static const char A
unsigned int get_degree() const
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
Definition: manifold.cc:123
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: cuda.h:32
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
size_type size(const unsigned int i) const
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
std::vector< Point< 1 > > line_support_points
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
virtual bool preserves_vertex_locations() const override
VectorType::value_type * begin(VectorType &V)
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
T min(const T &t, const MPI_Comm &mpi_communicator)
Normal vectors.
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
virtual std::size_t memory_consumption() const override
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
static ::ExceptionBase & ExcNotImplemented()
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
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< unsigned int > lexicographic_numbering
Definition: shape_info.h:357
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2416
std::vector< double > shape_values
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
EvaluationFlags
The EvaluationFlags enum.
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
bool is_tensor_product() const
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:9481
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701
Covariant transformation.
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
std::vector< Tensor< 1, spacedim > > normal_vectors
internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< double > > shape_info
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
UpdateFlags update_each
Definition: mapping.h:642
static ::ExceptionBase & ExcInternalError()