Reference documentation for deal.II version Git 50c3491829 2021-08-01 13:40:40 +0200
\(\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.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
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_q.h>
30 #include <deal.II/fe/mapping_q1.h>
32 
34 #include <deal.II/grid/tria.h>
36 
38 #include <boost/container/small_vector.hpp>
40 
41 #include <algorithm>
42 #include <array>
43 #include <cmath>
44 #include <memory>
45 #include <numeric>
46 
47 
49 
50 
51 template <int dim, int spacedim>
53  const unsigned int polynomial_degree)
54  : polynomial_degree(polynomial_degree)
55  , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
56  , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
57  , tensor_product_quadrature(false)
58 {}
59 
60 
61 
62 template <int dim, int spacedim>
63 std::size_t
65 {
66  return (
79 }
80 
81 
82 
83 template <int dim, int spacedim>
84 void
86  const UpdateFlags update_flags,
87  const Quadrature<dim> &q,
88  const unsigned int n_original_q_points)
89 {
90  // store the flags in the internal data object so we can access them
91  // in fill_fe_*_values()
92  this->update_each = update_flags;
93 
94  const unsigned int n_q_points = q.size();
95 
96  const bool needs_higher_order_terms =
97  this->update_each &
102 
104  covariant.resize(n_original_q_points);
105 
107  contravariant.resize(n_original_q_points);
108 
110  volume_elements.resize(n_original_q_points);
111 
113 
114  // use of MatrixFree only for higher order elements and with more than one
115  // point where tensor products do not make sense
116  if (polynomial_degree < 2 || n_q_points == 1)
118 
119  if (dim > 1)
120  {
121  // find out if the one-dimensional formula is the same
122  // in all directions
124  {
125  const std::array<Quadrature<1>, dim> quad_array =
126  q.get_tensor_basis();
127  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
128  {
129  if (quad_array[i - 1].size() != quad_array[i].size())
130  {
131  tensor_product_quadrature = false;
132  break;
133  }
134  else
135  {
136  const std::vector<Point<1>> &points_1 =
137  quad_array[i - 1].get_points();
138  const std::vector<Point<1>> &points_2 =
139  quad_array[i].get_points();
140  const std::vector<double> &weights_1 =
141  quad_array[i - 1].get_weights();
142  const std::vector<double> &weights_2 =
143  quad_array[i].get_weights();
144  for (unsigned int j = 0; j < quad_array[i].size(); ++j)
145  {
146  if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
147  std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
148  {
149  tensor_product_quadrature = false;
150  break;
151  }
152  }
153  }
154  }
155 
156  if (tensor_product_quadrature)
157  {
158  // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
159  // numbering manually (building an FE_Q<dim> is relatively
160  // expensive due to constraints)
161  const FE_DGQ<1> fe(polynomial_degree);
162  shape_info.reinit(q.get_tensor_basis()[0], fe);
164  FETools::lexicographic_to_hierarchic_numbering<dim>(
166  shape_info.n_q_points = q.size();
169  }
170  }
171  }
172 
173  // Only fill the big arrays on demand in case we cannot use the tensor
174  // product quadrature code path
175  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
176  {
177  // see if we need the (transformation) shape function values
178  // and/or gradients and resize the necessary arrays
180  shape_values.resize(n_shape_functions * n_q_points);
181 
182  if (this->update_each &
183  (update_covariant_transformation |
184  update_contravariant_transformation | update_JxW_values |
193 
194  if (this->update_each &
197 
201 
205 
206  // now also fill the various fields with their correct values
208  }
209 }
210 
211 
212 
213 template <int dim, int spacedim>
214 void
216  const UpdateFlags update_flags,
217  const Quadrature<dim> &q,
218  const unsigned int n_original_q_points)
219 {
220  initialize(update_flags, q, n_original_q_points);
221 
222  if (dim > 1 && tensor_product_quadrature)
223  {
224  constexpr unsigned int facedim = dim - 1;
225  const FE_DGQ<1> fe(polynomial_degree);
226  shape_info.reinit(q.get_tensor_basis()[0], fe);
228  FETools::lexicographic_to_hierarchic_numbering<facedim>(
230  shape_info.n_q_points = n_original_q_points;
233  }
234 
235  if (dim > 1)
236  {
237  if (this->update_each &
240  {
241  aux.resize(dim - 1,
242  AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
243 
244  // Compute tangentials to the unit cell.
245  for (const unsigned int i : GeometryInfo<dim>::face_indices())
246  {
247  unit_tangentials[i].resize(n_original_q_points);
248  std::fill(unit_tangentials[i].begin(),
249  unit_tangentials[i].end(),
251  if (dim > 2)
252  {
254  .resize(n_original_q_points);
255  std::fill(
257  .begin(),
259  .end(),
261  }
262  }
263  }
264  }
265 }
266 
267 
268 
269 template <int dim, int spacedim>
270 void
272  const std::vector<Point<dim>> &unit_points)
273 {
274  const unsigned int n_points = unit_points.size();
275 
276  // Construct the tensor product polynomials used as shape functions for
277  // the Qp mapping of cells at the boundary.
278  const TensorProductPolynomials<dim> tensor_pols(
281  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
282 
283  // then also construct the mapping from lexicographic to the Qp shape
284  // function numbering
285  const std::vector<unsigned int> renumber =
286  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
287 
288  std::vector<double> values;
289  std::vector<Tensor<1, dim>> grads;
290  if (shape_values.size() != 0)
291  {
292  Assert(shape_values.size() == n_shape_functions * n_points,
293  ExcInternalError());
294  values.resize(n_shape_functions);
295  }
296  if (shape_derivatives.size() != 0)
297  {
299  ExcInternalError());
300  grads.resize(n_shape_functions);
301  }
302 
303  std::vector<Tensor<2, dim>> grad2;
304  if (shape_second_derivatives.size() != 0)
305  {
307  ExcInternalError());
308  grad2.resize(n_shape_functions);
309  }
310 
311  std::vector<Tensor<3, dim>> grad3;
312  if (shape_third_derivatives.size() != 0)
313  {
315  ExcInternalError());
316  grad3.resize(n_shape_functions);
317  }
318 
319  std::vector<Tensor<4, dim>> grad4;
320  if (shape_fourth_derivatives.size() != 0)
321  {
323  ExcInternalError());
324  grad4.resize(n_shape_functions);
325  }
326 
327 
328  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
330  shape_third_derivatives.size() != 0 ||
332  for (unsigned int point = 0; point < n_points; ++point)
333  {
334  tensor_pols.evaluate(
335  unit_points[point], values, grads, grad2, grad3, grad4);
336 
337  if (shape_values.size() != 0)
338  for (unsigned int i = 0; i < n_shape_functions; ++i)
339  shape(point, i) = values[renumber[i]];
340 
341  if (shape_derivatives.size() != 0)
342  for (unsigned int i = 0; i < n_shape_functions; ++i)
343  derivative(point, i) = grads[renumber[i]];
344 
345  if (shape_second_derivatives.size() != 0)
346  for (unsigned int i = 0; i < n_shape_functions; ++i)
347  second_derivative(point, i) = grad2[renumber[i]];
348 
349  if (shape_third_derivatives.size() != 0)
350  for (unsigned int i = 0; i < n_shape_functions; ++i)
351  third_derivative(point, i) = grad3[renumber[i]];
352 
353  if (shape_fourth_derivatives.size() != 0)
354  for (unsigned int i = 0; i < n_shape_functions; ++i)
355  fourth_derivative(point, i) = grad4[renumber[i]];
356  }
357 }
358 
359 
360 
361 template <int dim, int spacedim>
363  : polynomial_degree(p)
365  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
366  , polynomials_1d(
371  internal::MappingQImplementation::unit_support_points<dim>(
375  internal::MappingQImplementation::
377  this->polynomial_degree,
378  dim))
380  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
381  this->polynomial_degree))
382 {
383  Assert(p >= 1,
384  ExcMessage("It only makes sense to create polynomial mappings "
385  "with a polynomial degree greater or equal to one."));
386 }
387 
388 
389 
390 template <int dim, int spacedim>
391 MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
392  : polynomial_degree(p)
394  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
395  , polynomials_1d(
400  internal::MappingQImplementation::unit_support_points<dim>(
404  internal::MappingQImplementation::
406  this->polynomial_degree,
407  dim))
409  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
410  this->polynomial_degree))
411 {
412  Assert(p >= 1,
413  ExcMessage("It only makes sense to create polynomial mappings "
414  "with a polynomial degree greater or equal to one."));
415 }
416 
417 
418 
419 template <int dim, int spacedim>
423  , polynomials_1d(mapping.polynomials_1d)
429 {}
430 
431 
432 
433 template <int dim, int spacedim>
434 std::unique_ptr<Mapping<dim, spacedim>>
436 {
437  return std::make_unique<MappingQ<dim, spacedim>>(*this);
438 }
439 
440 
441 
442 template <int dim, int spacedim>
443 unsigned int
445 {
446  return polynomial_degree;
447 }
448 
449 
450 
451 template <int dim, int spacedim>
454  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
455  const Point<dim> & p) const
456 {
459  this->compute_mapping_support_points(cell),
460  p,
461  polynomials_1d.size() == 2,
463  .first);
464 }
465 
466 
467 // In the code below, GCC tries to instantiate MappingQ<3,4> when
468 // seeing which of the overloaded versions of
469 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
470 // error messages and, generally, nothing very good. Avoid this by ensuring
471 // that this class exists, but does not have an inner InternalData
472 // type, thereby ruling out the codim-1 version of the function
473 // below when doing overload resolution.
474 template <>
475 class MappingQ<3, 4>
476 {};
477 
478 
479 
480 // visual studio freaks out when trying to determine if
481 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
482 // candidate. So instead of letting the compiler pick the correct overload, we
483 // use template specialization to make sure we pick up the right function to
484 // call:
485 
486 template <int dim, int spacedim>
490  const Point<spacedim> &,
491  const Point<dim> &) const
492 {
493  // default implementation (should never be called)
494  Assert(false, ExcInternalError());
495  return {};
496 }
497 
498 
499 
500 template <>
501 Point<1>
504  const Point<1> & p,
505  const Point<1> & initial_p_unit) const
506 {
507  // dispatch to the various specializations for spacedim=dim,
508  // spacedim=dim+1, etc
509  return internal::MappingQImplementation::
510  do_transform_real_to_unit_cell_internal<1>(
511  p,
512  initial_p_unit,
513  this->compute_mapping_support_points(cell),
516 }
517 
518 
519 
520 template <>
521 Point<2>
524  const Point<2> & p,
525  const Point<2> & initial_p_unit) const
526 {
527  return internal::MappingQImplementation::
528  do_transform_real_to_unit_cell_internal<2>(
529  p,
530  initial_p_unit,
531  this->compute_mapping_support_points(cell),
534 }
535 
536 
537 
538 template <>
539 Point<3>
542  const Point<3> & p,
543  const Point<3> & initial_p_unit) const
544 {
545  return internal::MappingQImplementation::
546  do_transform_real_to_unit_cell_internal<3>(
547  p,
548  initial_p_unit,
549  this->compute_mapping_support_points(cell),
552 }
553 
554 
555 
556 template <>
557 Point<1>
560  const Point<2> & p,
561  const Point<1> & initial_p_unit) const
562 {
563  const int dim = 1;
564  const int spacedim = 2;
565 
566  const Quadrature<dim> point_quadrature(initial_p_unit);
567 
569  if (spacedim > dim)
570  update_flags |= update_jacobian_grads;
572  get_data(update_flags, point_quadrature));
573 
575 
576  // dispatch to the various specializations for spacedim=dim,
577  // spacedim=dim+1, etc
578  return internal::MappingQImplementation::
579  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
580  p,
581  initial_p_unit,
582  *mdata);
583 }
584 
585 
586 
587 template <>
588 Point<2>
591  const Point<3> & p,
592  const Point<2> & initial_p_unit) const
593 {
594  const int dim = 2;
595  const int spacedim = 3;
596 
597  const Quadrature<dim> point_quadrature(initial_p_unit);
598 
600  if (spacedim > dim)
601  update_flags |= update_jacobian_grads;
603  get_data(update_flags, point_quadrature));
604 
606 
607  // dispatch to the various specializations for spacedim=dim,
608  // spacedim=dim+1, etc
609  return internal::MappingQImplementation::
610  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
611  p,
612  initial_p_unit,
613  *mdata);
614 }
615 
616 template <>
617 Point<1>
620  const Point<3> &,
621  const Point<1> &) const
622 {
623  Assert(false, ExcNotImplemented());
624  return {};
625 }
626 
627 
628 
629 template <int dim, int spacedim>
632  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
633  const Point<spacedim> & p) const
634 {
635  // Use an exact formula if one is available. this is only the case
636  // for Q1 mappings in 1d, and in 2d if dim==spacedim
637  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
638  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
639  {
640  // The dimension-dependent algorithms are much faster (about 25-45x in
641  // 2D) but fail most of the time when the given point (p) is not in the
642  // cell. The dimension-independent Newton algorithm given below is
643  // slower, but more robust (though it still sometimes fails). Therefore
644  // this function implements the following strategy based on the
645  // p's dimension:
646  //
647  // * In 1D this mapping is linear, so the mapping is always invertible
648  // (and the exact formula is known) as long as the cell has non-zero
649  // length.
650  // * In 2D the exact (quadratic) formula is called first. If either the
651  // exact formula does not succeed (negative discriminant in the
652  // quadratic formula) or succeeds but finds a solution outside of the
653  // unit cell, then the Newton solver is called. The rationale for the
654  // second choice is that the exact formula may provide two different
655  // answers when mapping a point outside of the real cell, but the
656  // Newton solver (if it converges) will only return one answer.
657  // Otherwise the exact formula successfully found a point in the unit
658  // cell and that value is returned.
659  // * In 3D there is no (known to the authors) exact formula, so the Newton
660  // algorithm is used.
661  const auto vertices_ = this->get_vertices(cell);
662 
663  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
664  vertices;
665  for (unsigned int i = 0; i < vertices.size(); ++i)
666  vertices[i] = vertices_[i];
667 
668  try
669  {
670  switch (dim)
671  {
672  case 1:
673  {
674  // formula not subject to any issues in 1d
675  if (spacedim == 1)
677  vertices, p);
678  else
679  break;
680  }
681 
682  case 2:
683  {
684  const Point<dim> point =
686  p);
687 
688  // formula not guaranteed to work for points outside of
689  // the cell. only take the computed point if it lies
690  // inside the reference cell
691  const double eps = 1e-15;
692  if (-eps <= point(1) && point(1) <= 1 + eps &&
693  -eps <= point(0) && point(0) <= 1 + eps)
694  {
695  return point;
696  }
697  else
698  break;
699  }
700 
701  default:
702  {
703  // we should get here, based on the if-condition at the top
704  Assert(false, ExcInternalError());
705  }
706  }
707  }
708  catch (
710  {
711  // simply fall through and continue on to the standard Newton code
712  }
713  }
714  else
715  {
716  // we can't use an explicit formula,
717  }
718 
719 
720  // Find the initial value for the Newton iteration by a normal
721  // projection to the least square plane determined by the vertices
722  // of the cell
723  Point<dim> initial_p_unit;
724  if (this->preserves_vertex_locations())
725  {
726  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
727  // in 1d with spacedim > 1 the affine approximation is exact
728  if (dim == 1 && polynomial_degree == 1)
729  return initial_p_unit;
730  }
731  else
732  {
733  // else, we simply use the mid point
734  for (unsigned int d = 0; d < dim; ++d)
735  initial_p_unit[d] = 0.5;
736  }
737 
738  // perform the Newton iteration and return the result. note that this
739  // statement may throw an exception, which we simply pass up to the caller
740  const Point<dim> p_unit =
741  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
742  if (p_unit[0] == std::numeric_limits<double>::infinity())
743  AssertThrow(false,
745  return p_unit;
746 }
747 
748 
749 
750 template <int dim, int spacedim>
751 void
753  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
754  const ArrayView<const Point<spacedim>> & real_points,
755  const ArrayView<Point<dim>> & unit_points) const
756 {
757  // Go to base class functions for dim < spacedim because it is not yet
758  // implemented with optimized code.
759  if (dim < spacedim)
760  {
762  real_points,
763  unit_points);
764  return;
765  }
766 
767  AssertDimension(real_points.size(), unit_points.size());
768  const std::vector<Point<spacedim>> support_points =
769  this->compute_mapping_support_points(cell);
770 
771  // From the given (high-order) support points, now only pick the first
772  // 2^dim points and construct an affine approximation from those.
774  inverse_approximation(support_points, unit_cell_support_points);
775 
776  const unsigned int n_points = real_points.size();
777  const unsigned int n_lanes = VectorizedArray<double>::size();
778 
779  // Use the more heavy VectorizedArray code path if there is more than
780  // one point left to compute
781  for (unsigned int i = 0; i < n_points; i += n_lanes)
782  if (n_points - i > 1)
783  {
785  for (unsigned int j = 0; j < n_lanes; ++j)
786  if (i + j < n_points)
787  for (unsigned int d = 0; d < spacedim; ++d)
788  p_vec[d][j] = real_points[i + j][d];
789  else
790  for (unsigned int d = 0; d < spacedim; ++d)
791  p_vec[d][j] = real_points[i][d];
792 
794  internal::MappingQImplementation::
795  do_transform_real_to_unit_cell_internal<dim, spacedim>(
796  p_vec,
797  inverse_approximation.compute(p_vec),
798  support_points,
801 
802  // If the vectorized computation failed, it could be that only some of
803  // the lanes failed but others would have succeeded if we had let them
804  // compute alone without interference (like negative Jacobian
805  // determinants) from other SIMD lanes. Repeat the computation in this
806  // unlikely case with scalar arguments.
807  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
808  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
809  unit_points[i + j] = internal::MappingQImplementation::
810  do_transform_real_to_unit_cell_internal<dim, spacedim>(
811  real_points[i + j],
812  inverse_approximation.compute(real_points[i + j]),
813  support_points,
816  else
817  for (unsigned int d = 0; d < dim; ++d)
818  unit_points[i + j][d] = unit_point[d][j];
819  }
820  else
821  unit_points[i] = internal::MappingQImplementation::
822  do_transform_real_to_unit_cell_internal<dim, spacedim>(
823  real_points[i],
824  inverse_approximation.compute(real_points[i]),
825  support_points,
828 }
829 
830 
831 
832 template <int dim, int spacedim>
835 {
836  // add flags if the respective quantities are necessary to compute
837  // what we need. note that some flags appear in both the conditions
838  // and in subsequent set operations. this leads to some circular
839  // logic. the only way to treat this is to iterate. since there are
840  // 5 if-clauses in the loop, it will take at most 5 iterations to
841  // converge. do them:
842  UpdateFlags out = in;
843  for (unsigned int i = 0; i < 5; ++i)
844  {
845  // The following is a little incorrect:
846  // If not applied on a face,
847  // update_boundary_forms does not
848  // make sense. On the other hand,
849  // it is necessary on a
850  // face. Currently,
851  // update_boundary_forms is simply
852  // ignored for the interior of a
853  // cell.
855  out |= update_boundary_forms;
856 
861 
862  if (out &
867 
868  // The contravariant transformation is used in the Piola
869  // transformation, which requires the determinant of the Jacobi
870  // matrix of the transformation. Because we have no way of
871  // knowing here whether the finite element wants to use the
872  // contravariant or the Piola transforms, we add the JxW values
873  // to the list of flags to be updated for each cell.
875  out |= update_volume_elements;
876 
877  // the same is true when computing normal vectors: they require
878  // the determinant of the Jacobian
879  if (out & update_normal_vectors)
880  out |= update_volume_elements;
881  }
882 
883  return out;
884 }
885 
886 
887 
888 template <int dim, int spacedim>
889 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
891  const Quadrature<dim> &q) const
892 {
893  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
894  std::make_unique<InternalData>(polynomial_degree);
895  auto &data = dynamic_cast<InternalData &>(*data_ptr);
896  data.initialize(this->requires_update_flags(update_flags), q, q.size());
897 
898  return data_ptr;
899 }
900 
901 
902 
903 template <int dim, int spacedim>
904 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
906  const UpdateFlags update_flags,
907  const hp::QCollection<dim - 1> &quadrature) const
908 {
909  AssertDimension(quadrature.size(), 1);
910 
911  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
912  std::make_unique<InternalData>(polynomial_degree);
913  auto &data = dynamic_cast<InternalData &>(*data_ptr);
914  data.initialize_face(this->requires_update_flags(update_flags),
916  ReferenceCells::get_hypercube<dim>(), quadrature[0]),
917  quadrature[0].size());
918 
919  return data_ptr;
920 }
921 
922 
923 
924 template <int dim, int spacedim>
925 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
927  const UpdateFlags update_flags,
928  const Quadrature<dim - 1> &quadrature) const
929 {
930  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
931  std::make_unique<InternalData>(polynomial_degree);
932  auto &data = dynamic_cast<InternalData &>(*data_ptr);
933  data.initialize_face(this->requires_update_flags(update_flags),
935  ReferenceCells::get_hypercube<dim>(), quadrature),
936  quadrature.size());
937 
938  return data_ptr;
939 }
940 
941 
942 
943 template <int dim, int spacedim>
946  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
947  const CellSimilarity::Similarity cell_similarity,
948  const Quadrature<dim> & quadrature,
949  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
951  &output_data) const
952 {
953  // ensure that the following static_cast is really correct:
954  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
955  ExcInternalError());
956  const InternalData &data = static_cast<const InternalData &>(internal_data);
957 
958  const unsigned int n_q_points = quadrature.size();
959 
960  // recompute the support points of the transformation of this
961  // cell. we tried to be clever here in an earlier version of the
962  // library by checking whether the cell is the same as the one we
963  // had visited last, but it turns out to be difficult to determine
964  // that because a cell for the purposes of a mapping is
965  // characterized not just by its (triangulation, level, index)
966  // triple, but also by the locations of its vertices, the manifold
967  // object attached to the cell and all of its bounding faces/edges,
968  // etc. to reliably test that the "cell" we are on is, therefore,
969  // not easily done
971  data.cell_of_current_support_points = cell;
972 
973  // if the order of the mapping is greater than 1, then do not reuse any cell
974  // similarity information. This is necessary because the cell similarity
975  // value is computed with just cell vertices and does not take into account
976  // cell curvature.
977  const CellSimilarity::Similarity computed_cell_similarity =
978  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
979 
980  if (dim > 1 && data.tensor_product_quadrature)
981  {
982  internal::MappingQImplementation::
983  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
984  computed_cell_similarity,
985  data,
986  output_data.quadrature_points,
987  output_data.jacobian_grads);
988  }
989  else
990  {
991  internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
993  data,
994  output_data.quadrature_points);
995 
996  internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
997  computed_cell_similarity,
999  data);
1000 
1002  spacedim>(
1003  computed_cell_similarity,
1005  data,
1006  output_data.jacobian_grads);
1007  }
1008 
1010  dim,
1011  spacedim>(computed_cell_similarity,
1013  data,
1014  output_data.jacobian_pushed_forward_grads);
1015 
1017  dim,
1018  spacedim>(computed_cell_similarity,
1020  data,
1021  output_data.jacobian_2nd_derivatives);
1022 
1023  internal::MappingQImplementation::
1024  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1025  computed_cell_similarity,
1027  data,
1029 
1031  dim,
1032  spacedim>(computed_cell_similarity,
1034  data,
1035  output_data.jacobian_3rd_derivatives);
1036 
1037  internal::MappingQImplementation::
1038  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1039  computed_cell_similarity,
1041  data,
1043 
1044  const UpdateFlags update_flags = data.update_each;
1045  const std::vector<double> &weights = quadrature.get_weights();
1046 
1047  // Multiply quadrature weights by absolute value of Jacobian determinants or
1048  // the area element g=sqrt(DX^t DX) in case of codim > 0
1049 
1050  if (update_flags & (update_normal_vectors | update_JxW_values))
1051  {
1052  AssertDimension(output_data.JxW_values.size(), n_q_points);
1053 
1054  Assert(!(update_flags & update_normal_vectors) ||
1055  (output_data.normal_vectors.size() == n_q_points),
1056  ExcDimensionMismatch(output_data.normal_vectors.size(),
1057  n_q_points));
1058 
1059 
1060  if (computed_cell_similarity != CellSimilarity::translation)
1061  for (unsigned int point = 0; point < n_q_points; ++point)
1062  {
1063  if (dim == spacedim)
1064  {
1065  const double det = data.contravariant[point].determinant();
1066 
1067  // check for distorted cells.
1068 
1069  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1070  // 1e12 in 2D. might want to find a finer
1071  // (dimension-independent) criterion
1072  Assert(det >
1073  1e-12 * Utilities::fixed_power<dim>(
1074  cell->diameter() / std::sqrt(double(dim))),
1076  cell->center(), det, point)));
1077 
1078  output_data.JxW_values[point] = weights[point] * det;
1079  }
1080  // if dim==spacedim, then there is no cell normal to
1081  // compute. since this is for FEValues (and not FEFaceValues),
1082  // there are also no face normals to compute
1083  else // codim>0 case
1084  {
1085  Tensor<1, spacedim> DX_t[dim];
1086  for (unsigned int i = 0; i < spacedim; ++i)
1087  for (unsigned int j = 0; j < dim; ++j)
1088  DX_t[j][i] = data.contravariant[point][i][j];
1089 
1090  Tensor<2, dim> G; // First fundamental form
1091  for (unsigned int i = 0; i < dim; ++i)
1092  for (unsigned int j = 0; j < dim; ++j)
1093  G[i][j] = DX_t[i] * DX_t[j];
1094 
1095  output_data.JxW_values[point] =
1096  std::sqrt(determinant(G)) * weights[point];
1097 
1098  if (computed_cell_similarity ==
1100  {
1101  // we only need to flip the normal
1102  if (update_flags & update_normal_vectors)
1103  output_data.normal_vectors[point] *= -1.;
1104  }
1105  else
1106  {
1107  if (update_flags & update_normal_vectors)
1108  {
1109  Assert(spacedim == dim + 1,
1110  ExcMessage(
1111  "There is no (unique) cell normal for " +
1113  "-dimensional cells in " +
1114  Utilities::int_to_string(spacedim) +
1115  "-dimensional space. This only works if the "
1116  "space dimension is one greater than the "
1117  "dimensionality of the mesh cells."));
1118 
1119  if (dim == 1)
1120  output_data.normal_vectors[point] =
1121  cross_product_2d(-DX_t[0]);
1122  else // dim == 2
1123  output_data.normal_vectors[point] =
1124  cross_product_3d(DX_t[0], DX_t[1]);
1125 
1126  output_data.normal_vectors[point] /=
1127  output_data.normal_vectors[point].norm();
1128 
1129  if (cell->direction_flag() == false)
1130  output_data.normal_vectors[point] *= -1.;
1131  }
1132  }
1133  } // codim>0 case
1134  }
1135  }
1136 
1137 
1138 
1139  // copy values from InternalData to vector given by reference
1140  if (update_flags & update_jacobians)
1141  {
1142  AssertDimension(output_data.jacobians.size(), n_q_points);
1143  if (computed_cell_similarity != CellSimilarity::translation)
1144  for (unsigned int point = 0; point < n_q_points; ++point)
1145  output_data.jacobians[point] = data.contravariant[point];
1146  }
1147 
1148  // copy values from InternalData to vector given by reference
1149  if (update_flags & update_inverse_jacobians)
1150  {
1151  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1152  if (computed_cell_similarity != CellSimilarity::translation)
1153  for (unsigned int point = 0; point < n_q_points; ++point)
1154  output_data.inverse_jacobians[point] =
1155  data.covariant[point].transpose();
1156  }
1157 
1158  return computed_cell_similarity;
1159 }
1160 
1161 
1162 
1163 template <int dim, int spacedim>
1164 void
1166  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1167  const unsigned int face_no,
1168  const hp::QCollection<dim - 1> & quadrature,
1169  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1171  &output_data) const
1172 {
1173  AssertDimension(quadrature.size(), 1);
1174 
1175  // ensure that the following cast is really correct:
1176  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1177  ExcInternalError());
1178  const InternalData &data = static_cast<const InternalData &>(internal_data);
1179 
1180  // if necessary, recompute the support points of the transformation of this
1181  // cell (note that we need to first check the triangulation pointer, since
1182  // otherwise the second test might trigger an exception if the triangulations
1183  // are not the same)
1184  if ((data.mapping_support_points.size() == 0) ||
1185  (&cell->get_triangulation() !=
1186  &data.cell_of_current_support_points->get_triangulation()) ||
1187  (cell != data.cell_of_current_support_points))
1188  {
1190  data.cell_of_current_support_points = cell;
1191  }
1192 
1194  *this,
1195  cell,
1196  face_no,
1199  ReferenceCells::get_hypercube<dim>(),
1200  face_no,
1201  cell->face_orientation(face_no),
1202  cell->face_flip(face_no),
1203  cell->face_rotation(face_no),
1204  quadrature[0].size()),
1205  quadrature[0],
1206  data,
1207  output_data);
1208 }
1209 
1210 
1211 
1212 template <int dim, int spacedim>
1213 void
1215  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1216  const unsigned int face_no,
1217  const unsigned int subface_no,
1218  const Quadrature<dim - 1> & quadrature,
1219  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1221  &output_data) const
1222 {
1223  // ensure that the following cast is really correct:
1224  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1225  ExcInternalError());
1226  const InternalData &data = static_cast<const InternalData &>(internal_data);
1227 
1228  // if necessary, recompute the support points of the transformation of this
1229  // cell (note that we need to first check the triangulation pointer, since
1230  // otherwise the second test might trigger an exception if the triangulations
1231  // are not the same)
1232  if ((data.mapping_support_points.size() == 0) ||
1233  (&cell->get_triangulation() !=
1234  &data.cell_of_current_support_points->get_triangulation()) ||
1235  (cell != data.cell_of_current_support_points))
1236  {
1238  data.cell_of_current_support_points = cell;
1239  }
1240 
1242  *this,
1243  cell,
1244  face_no,
1245  subface_no,
1247  ReferenceCells::get_hypercube<dim>(),
1248  face_no,
1249  subface_no,
1250  cell->face_orientation(face_no),
1251  cell->face_flip(face_no),
1252  cell->face_rotation(face_no),
1253  quadrature.size(),
1254  cell->subface_case(face_no)),
1255  quadrature,
1256  data,
1257  output_data);
1258 }
1259 
1260 
1261 
1262 template <int dim, int spacedim>
1263 void
1265  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1266  const ArrayView<const Point<dim>> & unit_points,
1267  const UpdateFlags update_flags,
1269  &output_data) const
1270 {
1271  if (update_flags == update_default)
1272  return;
1273 
1274  Assert(update_flags & update_inverse_jacobians ||
1275  update_flags & update_jacobians ||
1276  update_flags & update_quadrature_points,
1277  ExcNotImplemented());
1278 
1279  output_data.initialize(unit_points.size(), update_flags);
1280  const std::vector<Point<spacedim>> support_points =
1281  this->compute_mapping_support_points(cell);
1282 
1283  const unsigned int n_points = unit_points.size();
1284  const unsigned int n_lanes = VectorizedArray<double>::size();
1285 
1286  // Use the more heavy VectorizedArray code path if there is more than
1287  // one point left to compute
1288  for (unsigned int i = 0; i < n_points; i += n_lanes)
1289  if (n_points - i > 1)
1290  {
1292  for (unsigned int j = 0; j < n_lanes; ++j)
1293  if (i + j < n_points)
1294  for (unsigned int d = 0; d < dim; ++d)
1295  p_vec[d][j] = unit_points[i + j][d];
1296  else
1297  for (unsigned int d = 0; d < dim; ++d)
1298  p_vec[d][j] = unit_points[i][d];
1299 
1300  const auto result =
1303  support_points,
1304  p_vec,
1305  polynomial_degree == 1,
1307 
1308  if (update_flags & update_quadrature_points)
1309  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1310  for (unsigned int d = 0; d < spacedim; ++d)
1311  output_data.quadrature_points[i + j][d] = result.first[d][j];
1312 
1313  if (update_flags & update_jacobians)
1314  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1315  for (unsigned int d = 0; d < spacedim; ++d)
1316  for (unsigned int e = 0; e < dim; ++e)
1317  output_data.jacobians[i + j][d][e] = result.second[e][d][j];
1318 
1319  if (update_flags & update_inverse_jacobians)
1320  {
1322  result.second);
1324  inv_jac = jac.covariant_form();
1325  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1326  for (unsigned int d = 0; d < dim; ++d)
1327  for (unsigned int e = 0; e < spacedim; ++e)
1328  output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j];
1329  }
1330  }
1331  else
1332  {
1333  const auto result =
1336  support_points,
1337  unit_points[i],
1338  polynomial_degree == 1,
1340 
1341  if (update_flags & update_quadrature_points)
1342  output_data.quadrature_points[i] = result.first;
1343 
1344  if (update_flags & update_jacobians)
1345  {
1346  DerivativeForm<1, spacedim, dim> jac = result.second;
1347  output_data.jacobians[i] = jac.transpose();
1348  }
1349 
1350  if (update_flags & update_inverse_jacobians)
1351  {
1352  DerivativeForm<1, spacedim, dim> jac(result.second);
1354  for (unsigned int d = 0; d < dim; ++d)
1355  for (unsigned int e = 0; e < spacedim; ++e)
1356  output_data.inverse_jacobians[i][d][e] = inv_jac[d][e];
1357  }
1358  }
1359 }
1360 
1361 
1362 
1363 template <int dim, int spacedim>
1364 void
1366  const ArrayView<const Tensor<1, dim>> & input,
1367  const MappingKind mapping_kind,
1368  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1369  const ArrayView<Tensor<1, spacedim>> & output) const
1370 {
1372  mapping_kind,
1373  mapping_data,
1374  output);
1375 }
1376 
1377 
1378 
1379 template <int dim, int spacedim>
1380 void
1382  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1383  const MappingKind mapping_kind,
1384  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1385  const ArrayView<Tensor<2, spacedim>> & output) const
1386 {
1388  mapping_kind,
1389  mapping_data,
1390  output);
1391 }
1392 
1393 
1394 
1395 template <int dim, int spacedim>
1396 void
1398  const ArrayView<const Tensor<2, dim>> & input,
1399  const MappingKind mapping_kind,
1400  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1401  const ArrayView<Tensor<2, spacedim>> & output) const
1402 {
1403  switch (mapping_kind)
1404  {
1405  case mapping_contravariant:
1407  mapping_kind,
1408  mapping_data,
1409  output);
1410  return;
1411 
1416  mapping_kind,
1417  mapping_data,
1418  output);
1419  return;
1420  default:
1421  Assert(false, ExcNotImplemented());
1422  }
1423 }
1424 
1425 
1426 
1427 template <int dim, int spacedim>
1428 void
1430  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1431  const MappingKind mapping_kind,
1432  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1433  const ArrayView<Tensor<3, spacedim>> & output) const
1434 {
1435  AssertDimension(input.size(), output.size());
1436  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1437  ExcInternalError());
1438  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1439 
1440  switch (mapping_kind)
1441  {
1443  {
1446  "update_covariant_transformation"));
1447 
1448  for (unsigned int q = 0; q < output.size(); ++q)
1449  for (unsigned int i = 0; i < spacedim; ++i)
1450  for (unsigned int j = 0; j < spacedim; ++j)
1451  {
1452  double tmp[dim];
1453  for (unsigned int K = 0; K < dim; ++K)
1454  {
1455  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1456  for (unsigned int J = 1; J < dim; ++J)
1457  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1458  }
1459  for (unsigned int k = 0; k < spacedim; ++k)
1460  {
1461  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1462  for (unsigned int K = 1; K < dim; ++K)
1463  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1464  }
1465  }
1466  return;
1467  }
1468 
1469  default:
1470  Assert(false, ExcNotImplemented());
1471  }
1472 }
1473 
1474 
1475 
1476 template <int dim, int spacedim>
1477 void
1479  const ArrayView<const Tensor<3, dim>> & input,
1480  const MappingKind mapping_kind,
1481  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1482  const ArrayView<Tensor<3, spacedim>> & output) const
1483 {
1484  switch (mapping_kind)
1485  {
1486  case mapping_piola_hessian:
1490  mapping_kind,
1491  mapping_data,
1492  output);
1493  return;
1494  default:
1495  Assert(false, ExcNotImplemented());
1496  }
1497 }
1498 
1499 
1500 
1501 template <int dim, int spacedim>
1502 void
1504  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1505  std::vector<Point<spacedim>> & a) const
1506 {
1507  // if we only need the midpoint, then ask for it.
1508  if (this->polynomial_degree == 2)
1509  {
1510  for (unsigned int line_no = 0;
1511  line_no < GeometryInfo<dim>::lines_per_cell;
1512  ++line_no)
1513  {
1514  const typename Triangulation<dim, spacedim>::line_iterator line =
1515  (dim == 1 ?
1516  static_cast<
1518  cell->line(line_no));
1519 
1520  const Manifold<dim, spacedim> &manifold =
1521  ((line->manifold_id() == numbers::flat_manifold_id) &&
1522  (dim < spacedim) ?
1523  cell->get_manifold() :
1524  line->get_manifold());
1525  a.push_back(manifold.get_new_point_on_line(line));
1526  }
1527  }
1528  else
1529  // otherwise call the more complicated functions and ask for inner points
1530  // from the manifold description
1531  {
1532  std::vector<Point<spacedim>> tmp_points;
1533  for (unsigned int line_no = 0;
1534  line_no < GeometryInfo<dim>::lines_per_cell;
1535  ++line_no)
1536  {
1537  const typename Triangulation<dim, spacedim>::line_iterator line =
1538  (dim == 1 ?
1539  static_cast<
1541  cell->line(line_no));
1542 
1543  const Manifold<dim, spacedim> &manifold =
1544  ((line->manifold_id() == numbers::flat_manifold_id) &&
1545  (dim < spacedim) ?
1546  cell->get_manifold() :
1547  line->get_manifold());
1548 
1549  const std::array<Point<spacedim>, 2> vertices{
1550  {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
1551  cell->vertex(
1553 
1554  const std::size_t n_rows =
1556  a.resize(a.size() + n_rows);
1557  auto a_view = make_array_view(a.end() - n_rows, a.end());
1558  manifold.get_new_points(
1559  make_array_view(vertices.begin(), vertices.end()),
1561  a_view);
1562  }
1563  }
1564 }
1565 
1566 
1567 
1568 template <>
1569 void
1572  std::vector<Point<3>> & a) const
1573 {
1574  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1575 
1576  // used if face quad at boundary or entirely in the interior of the domain
1577  std::vector<Point<3>> tmp_points;
1578 
1579  // loop over all faces and collect points on them
1580  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1581  {
1582  const Triangulation<3>::face_iterator face = cell->face(face_no);
1583 
1584 #ifdef DEBUG
1585  const bool face_orientation = cell->face_orientation(face_no),
1586  face_flip = cell->face_flip(face_no),
1587  face_rotation = cell->face_rotation(face_no);
1588  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1589  lines_per_face = GeometryInfo<3>::lines_per_face;
1590 
1591  // some sanity checks up front
1592  for (unsigned int i = 0; i < vertices_per_face; ++i)
1593  Assert(face->vertex_index(i) ==
1594  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1595  face_no, i, face_orientation, face_flip, face_rotation)),
1596  ExcInternalError());
1597 
1598  // indices of the lines that bound a face are given by GeometryInfo<3>::
1599  // face_to_cell_lines
1600  for (unsigned int i = 0; i < lines_per_face; ++i)
1601  Assert(face->line(i) ==
1603  face_no, i, face_orientation, face_flip, face_rotation)),
1604  ExcInternalError());
1605 #endif
1606  // extract the points surrounding a quad from the points
1607  // already computed. First get the 4 vertices and then the points on
1608  // the four lines
1609  boost::container::small_vector<Point<3>, 200> tmp_points(
1612  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1613  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1614  if (polynomial_degree > 1)
1615  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1616  ++line)
1617  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1618  tmp_points[4 + line * (polynomial_degree - 1) + i] =
1620  (polynomial_degree - 1) *
1621  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
1622  i];
1623 
1624  const std::size_t n_rows =
1626  a.resize(a.size() + n_rows);
1627  auto a_view = make_array_view(a.end() - n_rows, a.end());
1628  face->get_manifold().get_new_points(
1629  make_array_view(tmp_points.begin(), tmp_points.end()),
1631  a_view);
1632  }
1633 }
1634 
1635 
1636 
1637 template <>
1638 void
1641  std::vector<Point<3>> & a) const
1642 {
1643  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1644  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1645  vertices[i] = cell->vertex(i);
1646 
1647  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1649  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1650  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1651  {
1652  Point<2> point(line_support_points[q1 + 1][0],
1653  line_support_points[q2 + 1][0]);
1654  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1655  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1656  }
1657 
1658  const std::size_t n_rows = weights.size(0);
1659  a.resize(a.size() + n_rows);
1660  auto a_view = make_array_view(a.end() - n_rows, a.end());
1661  cell->get_manifold().get_new_points(
1662  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1663 }
1664 
1665 
1666 
1667 template <int dim, int spacedim>
1668 void
1671  std::vector<Point<spacedim>> &) const
1672 {
1673  Assert(false, ExcInternalError());
1674 }
1675 
1676 
1677 
1678 template <int dim, int spacedim>
1679 std::vector<Point<spacedim>>
1681  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1682 {
1683  // get the vertices first
1684  std::vector<Point<spacedim>> a;
1685  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1686  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1687  a.push_back(cell->vertex(i));
1688 
1689  if (this->polynomial_degree > 1)
1690  {
1691  // check if all entities have the same manifold id which is when we can
1692  // simply ask the manifold for all points. the transfinite manifold can
1693  // do the interpolation better than this class, so if we detect that we
1694  // do not have to change anything here
1695  Assert(dim <= 3, ExcImpossibleInDim(dim));
1696  bool all_manifold_ids_are_equal = (dim == spacedim);
1697  if (all_manifold_ids_are_equal &&
1699  &cell->get_manifold()) == nullptr)
1700  {
1701  for (auto f : GeometryInfo<dim>::face_indices())
1702  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1703  all_manifold_ids_are_equal = false;
1704 
1705  if (dim == 3)
1706  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1707  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1708  all_manifold_ids_are_equal = false;
1709  }
1710 
1711  if (all_manifold_ids_are_equal)
1712  {
1713  const std::size_t n_rows = support_point_weights_cell.size(0);
1714  a.resize(a.size() + n_rows);
1715  auto a_view = make_array_view(a.end() - n_rows, a.end());
1716  cell->get_manifold().get_new_points(make_array_view(a.begin(),
1717  a.end() - n_rows),
1719  a_view);
1720  }
1721  else
1722  switch (dim)
1723  {
1724  case 1:
1725  add_line_support_points(cell, a);
1726  break;
1727  case 2:
1728  // in 2d, add the points on the four bounding lines to the
1729  // exterior (outer) points
1730  add_line_support_points(cell, a);
1731 
1732  // then get the interior support points
1733  if (dim != spacedim)
1734  add_quad_support_points(cell, a);
1735  else
1736  {
1737  const std::size_t n_rows =
1739  a.resize(a.size() + n_rows);
1740  auto a_view = make_array_view(a.end() - n_rows, a.end());
1741  cell->get_manifold().get_new_points(
1742  make_array_view(a.begin(), a.end() - n_rows),
1744  a_view);
1745  }
1746  break;
1747 
1748  case 3:
1749  // in 3d also add the points located on the boundary faces
1750  add_line_support_points(cell, a);
1751  add_quad_support_points(cell, a);
1752 
1753  // then compute the interior points
1754  {
1755  const std::size_t n_rows =
1757  a.resize(a.size() + n_rows);
1758  auto a_view = make_array_view(a.end() - n_rows, a.end());
1759  cell->get_manifold().get_new_points(
1760  make_array_view(a.begin(), a.end() - n_rows),
1762  a_view);
1763  }
1764  break;
1765 
1766  default:
1767  Assert(false, ExcNotImplemented());
1768  break;
1769  }
1770  }
1771 
1772  return a;
1773 }
1774 
1775 
1776 
1777 template <int dim, int spacedim>
1780  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1781 {
1783 }
1784 
1785 
1786 
1787 template <int dim, int spacedim>
1788 bool
1790  const ReferenceCell &reference_cell) const
1791 {
1792  Assert(dim == reference_cell.get_dimension(),
1793  ExcMessage("The dimension of your mapping (" +
1794  Utilities::to_string(dim) +
1795  ") and the reference cell cell_type (" +
1796  Utilities::to_string(reference_cell.get_dimension()) +
1797  " ) do not agree."));
1798 
1799  return reference_cell.is_hyper_cube();
1800 }
1801 
1802 
1803 
1804 //--------------------------- Explicit instantiations -----------------------
1805 #include "mapping_q.inst"
1806 
1807 
const unsigned int polynomial_degree
Definition: mapping_q.h:429
Transformed quadrature weights.
unsigned int get_degree() const
Definition: mapping_q.cc:444
const std::vector< Point< dim > > unit_cell_support_points
Definition: mapping_q.h:664
void resize(const size_type new_size)
AlignedVector< Tensor< 4, dim > > shape_fourth_derivatives
Definition: mapping_q.h:406
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
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
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:64
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
Definition: mapping_q.cc:752
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1659
bool is_hyper_cube() const
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
Definition: mapping_q.cc:271
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
AlignedVector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition: mapping_q.h:513
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1442
Contravariant transformation.
MappingQ(const unsigned int polynomial_degree)
Definition: mapping_q.cc:362
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)
InternalData(const unsigned int polynomial_degree)
Definition: mapping_q.cc:52
unsigned int size() const
Definition: collection.h:110
const std::vector< Point< dim > > & get_points() const
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
const std::vector< double > & get_weights() const
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
Definition: mapping.cc:87
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:890
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
Definition: mapping_q.h:423
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
Volume element.
Point< dim, Number > compute(const Point< spacedim, Number > &p) const
Definition: fe_dgq.h:110
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
AlignedVector< double > volume_elements
Definition: mapping_q.h:536
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
Definition: mapping_q.cc:1669
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
Determinant of the Jacobian.
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
Definition: mapping_q.cc:1365
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:322
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1577
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim >> &output)
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)
MappingKind
Definition: mapping.h:64
static DataSetDescriptor cell()
Definition: qprojector.h:563
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
void transform_gradients(const ArrayView< const Tensor< rank, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim >> &output)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition: mapping_q.h:652
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:2569
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:414
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:215
std::vector< AlignedVector< Tensor< 1, spacedim > > > aux
Definition: mapping_q.h:518
T fixed_power(const T t)
Definition: utilities.h:1081
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
Definition: mapping_q.cc:905
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:318
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:85
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
const unsigned int n_shape_functions
Definition: mapping_q.h:440
AlignedVector< double > shape_values
Definition: mapping_q.h:375
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
No update.
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition: mapping_q.cc:1165
#define Assert(cond, exc)
Definition: exceptions.h:1467
UpdateFlags
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
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)
const Table< 2, double > support_point_weights_cell
Definition: mapping_q.h:700
const std::vector< Point< 1 > > line_support_points
Definition: mapping_q.h:638
Abstract base class for mapping classes.
Definition: mapping.h:303
internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< double > > shape_info
Definition: mapping_q.h:458
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:698
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
Definition: fe_values.cc:2871
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
const unsigned int polynomial_degree
Definition: mapping_q.h:628
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
AlignedVector< Tensor< 2, dim > > shape_second_derivatives
Definition: mapping_q.h:390
VectorType::value_type * end(VectorType &V)
Point< 3 > vertices[4]
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_q.cc:1779
DerivativeForm< 1, spacedim, dim, Number > transpose() const
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition: mapping_q.h:645
QGaussLobatto< 1 > line_support_points
Definition: mapping_q.h:450
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
inline ::Table< 2, double > compute_support_point_weights_cell(const unsigned int polynomial_degree)
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition: mapping_q.h:530
std::vector< Point< spacedim > > mapping_support_points
Definition: mapping_q.h:524
Gradient of volume element.
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
Definition: mapping_q.cc:1214
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int size() const
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Point< 2 > first
Definition: grid_out.cc:4586
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition: mapping_q.cc:926
unsigned int get_dimension() const
std::vector< Point< spacedim > > quadrature_points
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition: mapping_q.cc:1789
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void fill_mapping_data_for_generic_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim >> &unit_points, const UpdateFlags update_flags, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
Definition: mapping_q.cc:1264
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:125
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_q.cc:1680
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
Definition: mapping_q.cc:453
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:1407
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
AlignedVector< Tensor< 3, dim > > shape_third_derivatives
Definition: mapping_q.h:398
virtual bool preserves_vertex_locations() const override
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:452
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
AlignedVector< Tensor< 1, dim > > shape_derivatives
Definition: mapping_q.h:382
Normal vectors.
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition: mapping_q.h:686
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
size_type size() const
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping.cc:33
static ::ExceptionBase & ExcNotImplemented()
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
Definition: mapping_q.cc:488
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:386
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:435
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
Definition: mapping_q.cc:1503
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void transform_hessians(const ArrayView< const Tensor< 3, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim >> &output)
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2544
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_q.cc:834
void transform_fields(const ArrayView< const Tensor< rank, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim >> &output)
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
Definition: mapping_q.cc:945
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
Definition: mapping_q.cc:631
bool is_tensor_product() const
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:11441
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701
Covariant transformation.
std::vector< Tensor< 1, spacedim > > normal_vectors
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
UpdateFlags update_each
Definition: mapping.h:652
static ::ExceptionBase & ExcInternalError()
AlignedVector< DerivativeForm< 1, dim, spacedim > > covariant
Definition: mapping_q.h:504