Reference documentation for deal.II version Git 5b8b897cb2 2021-04-22 22:24:19 -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>
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 |
192  shape_derivatives.resize(n_shape_functions * n_q_points);
193 
194  if (this->update_each &
196  shape_second_derivatives.resize(n_shape_functions * n_q_points);
197 
200  shape_third_derivatives.resize(n_shape_functions * n_q_points);
201 
204  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
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  std::vector<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  {
298  Assert(shape_derivatives.size() == n_shape_functions * n_points,
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  {
314  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
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 ||
329  shape_second_derivatives.size() != 0 ||
330  shape_third_derivatives.size() != 0 ||
331  shape_fourth_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::MappingQGenericImplementation::unit_support_points<dim>(
375  internal::MappingQGenericImplementation::
377  this->polynomial_degree,
378  dim))
380  internal::MappingQGenericImplementation::
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>
392  const MappingQGeneric<dim, spacedim> &mapping)
395  , polynomials_1d(mapping.polynomials_1d)
401 {}
402 
403 
404 
405 template <int dim, int spacedim>
406 std::unique_ptr<Mapping<dim, spacedim>>
408 {
409  return std::make_unique<MappingQGeneric<dim, spacedim>>(*this);
410 }
411 
412 
413 
414 template <int dim, int spacedim>
415 unsigned int
417 {
418  return polynomial_degree;
419 }
420 
421 
422 
423 template <int dim, int spacedim>
426  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
427  const Point<dim> & p) const
428 {
431  this->compute_mapping_support_points(cell),
432  p,
433  polynomials_1d.size() == 2,
435  .first);
436 }
437 
438 
439 // In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
440 // seeing which of the overloaded versions of
441 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
442 // error messages and, generally, nothing very good. Avoid this by ensuring
443 // that this class exists, but does not have an inner InternalData
444 // type, thereby ruling out the codim-1 version of the function
445 // below when doing overload resolution.
446 template <>
447 class MappingQGeneric<3, 4>
448 {};
449 
450 
451 
452 // visual studio freaks out when trying to determine if
453 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
454 // candidate. So instead of letting the compiler pick the correct overload, we
455 // use template specialization to make sure we pick up the right function to
456 // call:
457 
458 template <int dim, int spacedim>
462  const Point<spacedim> &,
463  const Point<dim> &) const
464 {
465  // default implementation (should never be called)
466  Assert(false, ExcInternalError());
467  return {};
468 }
469 
470 
471 
472 template <>
473 Point<1>
476  const Point<1> & p,
477  const Point<1> & initial_p_unit) const
478 {
479  // dispatch to the various specializations for spacedim=dim,
480  // spacedim=dim+1, etc
481  return internal::MappingQGenericImplementation::
482  do_transform_real_to_unit_cell_internal<1>(
483  p,
484  initial_p_unit,
485  this->compute_mapping_support_points(cell),
488 }
489 
490 
491 
492 template <>
493 Point<2>
496  const Point<2> & p,
497  const Point<2> & initial_p_unit) const
498 {
499  return internal::MappingQGenericImplementation::
500  do_transform_real_to_unit_cell_internal<2>(
501  p,
502  initial_p_unit,
503  this->compute_mapping_support_points(cell),
506 }
507 
508 
509 
510 template <>
511 Point<3>
514  const Point<3> & p,
515  const Point<3> & initial_p_unit) const
516 {
517  return internal::MappingQGenericImplementation::
518  do_transform_real_to_unit_cell_internal<3>(
519  p,
520  initial_p_unit,
521  this->compute_mapping_support_points(cell),
524 }
525 
526 
527 
528 template <>
529 Point<1>
532  const Point<2> & p,
533  const Point<1> & initial_p_unit) const
534 {
535  const int dim = 1;
536  const int spacedim = 2;
537 
538  const Quadrature<dim> point_quadrature(initial_p_unit);
539 
541  if (spacedim > dim)
542  update_flags |= update_jacobian_grads;
544  get_data(update_flags, point_quadrature));
545 
547 
548  // dispatch to the various specializations for spacedim=dim,
549  // spacedim=dim+1, etc
550  return internal::MappingQGenericImplementation::
551  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
552  p,
553  initial_p_unit,
554  *mdata);
555 }
556 
557 
558 
559 template <>
560 Point<2>
563  const Point<3> & p,
564  const Point<2> & initial_p_unit) const
565 {
566  const int dim = 2;
567  const int spacedim = 3;
568 
569  const Quadrature<dim> point_quadrature(initial_p_unit);
570 
572  if (spacedim > dim)
573  update_flags |= update_jacobian_grads;
575  get_data(update_flags, point_quadrature));
576 
578 
579  // dispatch to the various specializations for spacedim=dim,
580  // spacedim=dim+1, etc
581  return internal::MappingQGenericImplementation::
582  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
583  p,
584  initial_p_unit,
585  *mdata);
586 }
587 
588 template <>
589 Point<1>
592  const Point<3> &,
593  const Point<1> &) const
594 {
595  Assert(false, ExcNotImplemented());
596  return {};
597 }
598 
599 
600 
601 template <int dim, int spacedim>
604  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
605  const Point<spacedim> & p) const
606 {
607  // Use an exact formula if one is available. this is only the case
608  // for Q1 mappings in 1d, and in 2d if dim==spacedim
609  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
610  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
611  {
612  // The dimension-dependent algorithms are much faster (about 25-45x in
613  // 2D) but fail most of the time when the given point (p) is not in the
614  // cell. The dimension-independent Newton algorithm given below is
615  // slower, but more robust (though it still sometimes fails). Therefore
616  // this function implements the following strategy based on the
617  // p's dimension:
618  //
619  // * In 1D this mapping is linear, so the mapping is always invertible
620  // (and the exact formula is known) as long as the cell has non-zero
621  // length.
622  // * In 2D the exact (quadratic) formula is called first. If either the
623  // exact formula does not succeed (negative discriminant in the
624  // quadratic formula) or succeeds but finds a solution outside of the
625  // unit cell, then the Newton solver is called. The rationale for the
626  // second choice is that the exact formula may provide two different
627  // answers when mapping a point outside of the real cell, but the
628  // Newton solver (if it converges) will only return one answer.
629  // Otherwise the exact formula successfully found a point in the unit
630  // cell and that value is returned.
631  // * In 3D there is no (known to the authors) exact formula, so the Newton
632  // algorithm is used.
633  const auto vertices_ = this->get_vertices(cell);
634 
635  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
636  vertices;
637  for (unsigned int i = 0; i < vertices.size(); ++i)
638  vertices[i] = vertices_[i];
639 
640  try
641  {
642  switch (dim)
643  {
644  case 1:
645  {
646  // formula not subject to any issues in 1d
647  if (spacedim == 1)
649  vertices, p);
650  else
651  break;
652  }
653 
654  case 2:
655  {
656  const Point<dim> point =
658  p);
659 
660  // formula not guaranteed to work for points outside of
661  // the cell. only take the computed point if it lies
662  // inside the reference cell
663  const double eps = 1e-15;
664  if (-eps <= point(1) && point(1) <= 1 + eps &&
665  -eps <= point(0) && point(0) <= 1 + eps)
666  {
667  return point;
668  }
669  else
670  break;
671  }
672 
673  default:
674  {
675  // we should get here, based on the if-condition at the top
676  Assert(false, ExcInternalError());
677  }
678  }
679  }
680  catch (
682  {
683  // simply fall through and continue on to the standard Newton code
684  }
685  }
686  else
687  {
688  // we can't use an explicit formula,
689  }
690 
691 
692  // Find the initial value for the Newton iteration by a normal
693  // projection to the least square plane determined by the vertices
694  // of the cell
695  Point<dim> initial_p_unit;
696  if (this->preserves_vertex_locations())
697  {
698  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
699  // in 1d with spacedim > 1 the affine approximation is exact
700  if (dim == 1 && polynomial_degree == 1)
701  return initial_p_unit;
702  }
703  else
704  {
705  // else, we simply use the mid point
706  for (unsigned int d = 0; d < dim; ++d)
707  initial_p_unit[d] = 0.5;
708  }
709 
710  // perform the Newton iteration and return the result. note that this
711  // statement may throw an exception, which we simply pass up to the caller
712  const Point<dim> p_unit =
713  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
714  if (p_unit[0] == std::numeric_limits<double>::infinity())
715  AssertThrow(false,
717  return p_unit;
718 }
719 
720 
721 
722 template <int dim, int spacedim>
723 void
725  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
726  const ArrayView<const Point<spacedim>> & real_points,
727  const ArrayView<Point<dim>> & unit_points) const
728 {
729  // Go to base class functions for dim < spacedim because it is not yet
730  // implemented with optimized code.
731  if (dim < spacedim)
732  {
734  real_points,
735  unit_points);
736  return;
737  }
738 
739  AssertDimension(real_points.size(), unit_points.size());
740  const std::vector<Point<spacedim>> support_points =
741  this->compute_mapping_support_points(cell);
742 
743  // From the given (high-order) support points, now only pick the first
744  // 2^dim points and construct an affine approximation from those.
745  internal::MappingQGenericImplementation::
746  InverseQuadraticApproximation<dim, spacedim>
747  inverse_approximation(support_points, unit_cell_support_points);
748 
749  const unsigned int n_points = real_points.size();
750  const unsigned int n_lanes = VectorizedArray<double>::size();
751 
752  // Use the more heavy VectorizedArray code path if there is more than
753  // one point left to compute
754  for (unsigned int i = 0; i < n_points; i += n_lanes)
755  if (n_points - i > 1)
756  {
758  for (unsigned int j = 0; j < n_lanes; ++j)
759  if (i + j < n_points)
760  for (unsigned int d = 0; d < spacedim; ++d)
761  p_vec[d][j] = real_points[i + j][d];
762  else
763  for (unsigned int d = 0; d < spacedim; ++d)
764  p_vec[d][j] = real_points[i][d];
765 
767  internal::MappingQGenericImplementation::
768  do_transform_real_to_unit_cell_internal<dim, spacedim>(
769  p_vec,
770  inverse_approximation.compute(p_vec),
771  support_points,
774 
775  // If the vectorized computation failed, it could be that only some of
776  // the lanes failed but others would have succeeded if we had let them
777  // compute alone without interference (like negative Jacobian
778  // determinants) from other SIMD lanes. Repeat the computation in this
779  // unlikely case with scalar arguments.
780  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
781  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
782  unit_points[i + j] = internal::MappingQGenericImplementation::
783  do_transform_real_to_unit_cell_internal<dim, spacedim>(
784  real_points[i + j],
785  inverse_approximation.compute(real_points[i + j]),
786  support_points,
789  else
790  for (unsigned int d = 0; d < dim; ++d)
791  unit_points[i + j][d] = unit_point[d][j];
792  }
793  else
794  unit_points[i] = internal::MappingQGenericImplementation::
795  do_transform_real_to_unit_cell_internal<dim, spacedim>(
796  real_points[i],
797  inverse_approximation.compute(real_points[i]),
798  support_points,
801 }
802 
803 
804 
805 template <int dim, int spacedim>
808  const UpdateFlags in) const
809 {
810  // add flags if the respective quantities are necessary to compute
811  // what we need. note that some flags appear in both the conditions
812  // and in subsequent set operations. this leads to some circular
813  // logic. the only way to treat this is to iterate. since there are
814  // 5 if-clauses in the loop, it will take at most 5 iterations to
815  // converge. do them:
816  UpdateFlags out = in;
817  for (unsigned int i = 0; i < 5; ++i)
818  {
819  // The following is a little incorrect:
820  // If not applied on a face,
821  // update_boundary_forms does not
822  // make sense. On the other hand,
823  // it is necessary on a
824  // face. Currently,
825  // update_boundary_forms is simply
826  // ignored for the interior of a
827  // cell.
829  out |= update_boundary_forms;
830 
835 
836  if (out &
841 
842  // The contravariant transformation is used in the Piola
843  // transformation, which requires the determinant of the Jacobi
844  // matrix of the transformation. Because we have no way of
845  // knowing here whether the finite element wants to use the
846  // contravariant or the Piola transforms, we add the JxW values
847  // to the list of flags to be updated for each cell.
849  out |= update_volume_elements;
850 
851  // the same is true when computing normal vectors: they require
852  // the determinant of the Jacobian
853  if (out & update_normal_vectors)
854  out |= update_volume_elements;
855  }
856 
857  return out;
858 }
859 
860 
861 
862 template <int dim, int spacedim>
863 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
865  const Quadrature<dim> &q) const
866 {
867  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
868  std::make_unique<InternalData>(polynomial_degree);
869  auto &data = dynamic_cast<InternalData &>(*data_ptr);
870  data.initialize(this->requires_update_flags(update_flags), q, q.size());
871 
872  return data_ptr;
873 }
874 
875 
876 
877 template <int dim, int spacedim>
878 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
880  const UpdateFlags update_flags,
881  const hp::QCollection<dim - 1> &quadrature) const
882 {
883  AssertDimension(quadrature.size(), 1);
884 
885  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
886  std::make_unique<InternalData>(polynomial_degree);
887  auto &data = dynamic_cast<InternalData &>(*data_ptr);
888  data.initialize_face(this->requires_update_flags(update_flags),
890  ReferenceCells::get_hypercube<dim>(), quadrature[0]),
891  quadrature[0].size());
892 
893  return data_ptr;
894 }
895 
896 
897 
898 template <int dim, int spacedim>
899 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
901  const UpdateFlags update_flags,
902  const Quadrature<dim - 1> &quadrature) const
903 {
904  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
905  std::make_unique<InternalData>(polynomial_degree);
906  auto &data = dynamic_cast<InternalData &>(*data_ptr);
907  data.initialize_face(this->requires_update_flags(update_flags),
909  ReferenceCells::get_hypercube<dim>(), quadrature),
910  quadrature.size());
911 
912  return data_ptr;
913 }
914 
915 
916 
917 template <int dim, int spacedim>
920  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
921  const CellSimilarity::Similarity cell_similarity,
922  const Quadrature<dim> & quadrature,
923  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
925  &output_data) const
926 {
927  // ensure that the following static_cast is really correct:
928  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
929  ExcInternalError());
930  const InternalData &data = static_cast<const InternalData &>(internal_data);
931 
932  const unsigned int n_q_points = quadrature.size();
933 
934  // recompute the support points of the transformation of this
935  // cell. we tried to be clever here in an earlier version of the
936  // library by checking whether the cell is the same as the one we
937  // had visited last, but it turns out to be difficult to determine
938  // that because a cell for the purposes of a mapping is
939  // characterized not just by its (triangulation, level, index)
940  // triple, but also by the locations of its vertices, the manifold
941  // object attached to the cell and all of its bounding faces/edges,
942  // etc. to reliably test that the "cell" we are on is, therefore,
943  // not easily done
945  data.cell_of_current_support_points = cell;
946 
947  // if the order of the mapping is greater than 1, then do not reuse any cell
948  // similarity information. This is necessary because the cell similarity
949  // value is computed with just cell vertices and does not take into account
950  // cell curvature.
951  const CellSimilarity::Similarity computed_cell_similarity =
952  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
953 
954  if (dim > 1 && data.tensor_product_quadrature)
955  {
956  internal::MappingQGenericImplementation::
957  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
958  computed_cell_similarity,
959  data,
960  output_data.quadrature_points,
961  output_data.jacobian_grads);
962  }
963  else
964  {
966  spacedim>(
968  data,
969  output_data.quadrature_points);
970 
972  spacedim>(
973  computed_cell_similarity,
975  data);
976 
978  dim,
979  spacedim>(computed_cell_similarity,
981  data,
982  output_data.jacobian_grads);
983  }
984 
985  internal::MappingQGenericImplementation::
986  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
987  computed_cell_similarity,
989  data,
990  output_data.jacobian_pushed_forward_grads);
991 
992  internal::MappingQGenericImplementation::
993  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
994  computed_cell_similarity,
996  data,
997  output_data.jacobian_2nd_derivatives);
998 
999  internal::MappingQGenericImplementation::
1000  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1001  computed_cell_similarity,
1003  data,
1005 
1006  internal::MappingQGenericImplementation::
1007  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1008  computed_cell_similarity,
1010  data,
1011  output_data.jacobian_3rd_derivatives);
1012 
1013  internal::MappingQGenericImplementation::
1014  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1015  computed_cell_similarity,
1017  data,
1019 
1020  const UpdateFlags update_flags = data.update_each;
1021  const std::vector<double> &weights = quadrature.get_weights();
1022 
1023  // Multiply quadrature weights by absolute value of Jacobian determinants or
1024  // the area element g=sqrt(DX^t DX) in case of codim > 0
1025 
1026  if (update_flags & (update_normal_vectors | update_JxW_values))
1027  {
1028  AssertDimension(output_data.JxW_values.size(), n_q_points);
1029 
1030  Assert(!(update_flags & update_normal_vectors) ||
1031  (output_data.normal_vectors.size() == n_q_points),
1032  ExcDimensionMismatch(output_data.normal_vectors.size(),
1033  n_q_points));
1034 
1035 
1036  if (computed_cell_similarity != CellSimilarity::translation)
1037  for (unsigned int point = 0; point < n_q_points; ++point)
1038  {
1039  if (dim == spacedim)
1040  {
1041  const double det = data.contravariant[point].determinant();
1042 
1043  // check for distorted cells.
1044 
1045  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1046  // 1e12 in 2D. might want to find a finer
1047  // (dimension-independent) criterion
1048  Assert(det >
1049  1e-12 * Utilities::fixed_power<dim>(
1050  cell->diameter() / std::sqrt(double(dim))),
1052  cell->center(), det, point)));
1053 
1054  output_data.JxW_values[point] = weights[point] * det;
1055  }
1056  // if dim==spacedim, then there is no cell normal to
1057  // compute. since this is for FEValues (and not FEFaceValues),
1058  // there are also no face normals to compute
1059  else // codim>0 case
1060  {
1061  Tensor<1, spacedim> DX_t[dim];
1062  for (unsigned int i = 0; i < spacedim; ++i)
1063  for (unsigned int j = 0; j < dim; ++j)
1064  DX_t[j][i] = data.contravariant[point][i][j];
1065 
1066  Tensor<2, dim> G; // First fundamental form
1067  for (unsigned int i = 0; i < dim; ++i)
1068  for (unsigned int j = 0; j < dim; ++j)
1069  G[i][j] = DX_t[i] * DX_t[j];
1070 
1071  output_data.JxW_values[point] =
1072  std::sqrt(determinant(G)) * weights[point];
1073 
1074  if (computed_cell_similarity ==
1076  {
1077  // we only need to flip the normal
1078  if (update_flags & update_normal_vectors)
1079  output_data.normal_vectors[point] *= -1.;
1080  }
1081  else
1082  {
1083  if (update_flags & update_normal_vectors)
1084  {
1085  Assert(spacedim == dim + 1,
1086  ExcMessage(
1087  "There is no (unique) cell normal for " +
1089  "-dimensional cells in " +
1090  Utilities::int_to_string(spacedim) +
1091  "-dimensional space. This only works if the "
1092  "space dimension is one greater than the "
1093  "dimensionality of the mesh cells."));
1094 
1095  if (dim == 1)
1096  output_data.normal_vectors[point] =
1097  cross_product_2d(-DX_t[0]);
1098  else // dim == 2
1099  output_data.normal_vectors[point] =
1100  cross_product_3d(DX_t[0], DX_t[1]);
1101 
1102  output_data.normal_vectors[point] /=
1103  output_data.normal_vectors[point].norm();
1104 
1105  if (cell->direction_flag() == false)
1106  output_data.normal_vectors[point] *= -1.;
1107  }
1108  }
1109  } // codim>0 case
1110  }
1111  }
1112 
1113 
1114 
1115  // copy values from InternalData to vector given by reference
1116  if (update_flags & update_jacobians)
1117  {
1118  AssertDimension(output_data.jacobians.size(), n_q_points);
1119  if (computed_cell_similarity != CellSimilarity::translation)
1120  for (unsigned int point = 0; point < n_q_points; ++point)
1121  output_data.jacobians[point] = data.contravariant[point];
1122  }
1123 
1124  // copy values from InternalData to vector given by reference
1125  if (update_flags & update_inverse_jacobians)
1126  {
1127  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1128  if (computed_cell_similarity != CellSimilarity::translation)
1129  for (unsigned int point = 0; point < n_q_points; ++point)
1130  output_data.inverse_jacobians[point] =
1131  data.covariant[point].transpose();
1132  }
1133 
1134  return computed_cell_similarity;
1135 }
1136 
1137 
1138 
1139 template <int dim, int spacedim>
1140 void
1142  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1143  const unsigned int face_no,
1144  const hp::QCollection<dim - 1> & quadrature,
1145  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1147  &output_data) const
1148 {
1149  AssertDimension(quadrature.size(), 1);
1150 
1151  // ensure that the following cast is really correct:
1152  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1153  ExcInternalError());
1154  const InternalData &data = static_cast<const InternalData &>(internal_data);
1155 
1156  // if necessary, recompute the support points of the transformation of this
1157  // cell (note that we need to first check the triangulation pointer, since
1158  // otherwise the second test might trigger an exception if the triangulations
1159  // are not the same)
1160  if ((data.mapping_support_points.size() == 0) ||
1161  (&cell->get_triangulation() !=
1162  &data.cell_of_current_support_points->get_triangulation()) ||
1163  (cell != data.cell_of_current_support_points))
1164  {
1166  data.cell_of_current_support_points = cell;
1167  }
1168 
1170  *this,
1171  cell,
1172  face_no,
1175  ReferenceCells::get_hypercube<dim>(),
1176  face_no,
1177  cell->face_orientation(face_no),
1178  cell->face_flip(face_no),
1179  cell->face_rotation(face_no),
1180  quadrature[0].size()),
1181  quadrature[0],
1182  data,
1183  output_data);
1184 }
1185 
1186 
1187 
1188 template <int dim, int spacedim>
1189 void
1191  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1192  const unsigned int face_no,
1193  const unsigned int subface_no,
1194  const Quadrature<dim - 1> & quadrature,
1195  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1197  &output_data) const
1198 {
1199  // ensure that the following cast is really correct:
1200  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1201  ExcInternalError());
1202  const InternalData &data = static_cast<const InternalData &>(internal_data);
1203 
1204  // if necessary, recompute the support points of the transformation of this
1205  // cell (note that we need to first check the triangulation pointer, since
1206  // otherwise the second test might trigger an exception if the triangulations
1207  // are not the same)
1208  if ((data.mapping_support_points.size() == 0) ||
1209  (&cell->get_triangulation() !=
1210  &data.cell_of_current_support_points->get_triangulation()) ||
1211  (cell != data.cell_of_current_support_points))
1212  {
1214  data.cell_of_current_support_points = cell;
1215  }
1216 
1218  *this,
1219  cell,
1220  face_no,
1221  subface_no,
1223  ReferenceCells::get_hypercube<dim>(),
1224  face_no,
1225  subface_no,
1226  cell->face_orientation(face_no),
1227  cell->face_flip(face_no),
1228  cell->face_rotation(face_no),
1229  quadrature.size(),
1230  cell->subface_case(face_no)),
1231  quadrature,
1232  data,
1233  output_data);
1234 }
1235 
1236 
1237 
1238 template <int dim, int spacedim>
1239 void
1241  const ArrayView<const Tensor<1, dim>> & input,
1242  const MappingKind mapping_kind,
1243  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1244  const ArrayView<Tensor<1, spacedim>> & output) const
1245 {
1247  mapping_kind,
1248  mapping_data,
1249  output);
1250 }
1251 
1252 
1253 
1254 template <int dim, int spacedim>
1255 void
1257  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1258  const MappingKind mapping_kind,
1259  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1260  const ArrayView<Tensor<2, spacedim>> & output) const
1261 {
1263  input, mapping_kind, mapping_data, output);
1264 }
1265 
1266 
1267 
1268 template <int dim, int spacedim>
1269 void
1271  const ArrayView<const Tensor<2, dim>> & input,
1272  const MappingKind mapping_kind,
1273  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1274  const ArrayView<Tensor<2, spacedim>> & output) const
1275 {
1276  switch (mapping_kind)
1277  {
1278  case mapping_contravariant:
1280  mapping_kind,
1281  mapping_data,
1282  output);
1283  return;
1284 
1289  input, mapping_kind, mapping_data, output);
1290  return;
1291  default:
1292  Assert(false, ExcNotImplemented());
1293  }
1294 }
1295 
1296 
1297 
1298 template <int dim, int spacedim>
1299 void
1301  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1302  const MappingKind mapping_kind,
1303  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1304  const ArrayView<Tensor<3, spacedim>> & output) const
1305 {
1306  AssertDimension(input.size(), output.size());
1307  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1308  ExcInternalError());
1309  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1310 
1311  switch (mapping_kind)
1312  {
1314  {
1317  "update_covariant_transformation"));
1318 
1319  for (unsigned int q = 0; q < output.size(); ++q)
1320  for (unsigned int i = 0; i < spacedim; ++i)
1321  for (unsigned int j = 0; j < spacedim; ++j)
1322  {
1323  double tmp[dim];
1324  for (unsigned int K = 0; K < dim; ++K)
1325  {
1326  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1327  for (unsigned int J = 1; J < dim; ++J)
1328  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1329  }
1330  for (unsigned int k = 0; k < spacedim; ++k)
1331  {
1332  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1333  for (unsigned int K = 1; K < dim; ++K)
1334  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1335  }
1336  }
1337  return;
1338  }
1339 
1340  default:
1341  Assert(false, ExcNotImplemented());
1342  }
1343 }
1344 
1345 
1346 
1347 template <int dim, int spacedim>
1348 void
1350  const ArrayView<const Tensor<3, dim>> & input,
1351  const MappingKind mapping_kind,
1352  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1353  const ArrayView<Tensor<3, spacedim>> & output) const
1354 {
1355  switch (mapping_kind)
1356  {
1357  case mapping_piola_hessian:
1361  input, mapping_kind, mapping_data, output);
1362  return;
1363  default:
1364  Assert(false, ExcNotImplemented());
1365  }
1366 }
1367 
1368 
1369 
1370 template <int dim, int spacedim>
1371 void
1373  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1374  std::vector<Point<spacedim>> & a) const
1375 {
1376  // if we only need the midpoint, then ask for it.
1377  if (this->polynomial_degree == 2)
1378  {
1379  for (unsigned int line_no = 0;
1380  line_no < GeometryInfo<dim>::lines_per_cell;
1381  ++line_no)
1382  {
1383  const typename Triangulation<dim, spacedim>::line_iterator line =
1384  (dim == 1 ?
1385  static_cast<
1387  cell->line(line_no));
1388 
1389  const Manifold<dim, spacedim> &manifold =
1390  ((line->manifold_id() == numbers::flat_manifold_id) &&
1391  (dim < spacedim) ?
1392  cell->get_manifold() :
1393  line->get_manifold());
1394  a.push_back(manifold.get_new_point_on_line(line));
1395  }
1396  }
1397  else
1398  // otherwise call the more complicated functions and ask for inner points
1399  // from the manifold description
1400  {
1401  std::vector<Point<spacedim>> tmp_points;
1402  for (unsigned int line_no = 0;
1403  line_no < GeometryInfo<dim>::lines_per_cell;
1404  ++line_no)
1405  {
1406  const typename Triangulation<dim, spacedim>::line_iterator line =
1407  (dim == 1 ?
1408  static_cast<
1410  cell->line(line_no));
1411 
1412  const Manifold<dim, spacedim> &manifold =
1413  ((line->manifold_id() == numbers::flat_manifold_id) &&
1414  (dim < spacedim) ?
1415  cell->get_manifold() :
1416  line->get_manifold());
1417 
1418  const std::array<Point<spacedim>, 2> vertices{
1419  {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
1420  cell->vertex(
1422 
1423  const std::size_t n_rows =
1425  a.resize(a.size() + n_rows);
1426  auto a_view = make_array_view(a.end() - n_rows, a.end());
1427  manifold.get_new_points(
1428  make_array_view(vertices.begin(), vertices.end()),
1430  a_view);
1431  }
1432  }
1433 }
1434 
1435 
1436 
1437 template <>
1438 void
1441  std::vector<Point<3>> & a) const
1442 {
1443  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1444 
1445  // used if face quad at boundary or entirely in the interior of the domain
1446  std::vector<Point<3>> tmp_points;
1447 
1448  // loop over all faces and collect points on them
1449  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1450  {
1451  const Triangulation<3>::face_iterator face = cell->face(face_no);
1452 
1453 #ifdef DEBUG
1454  const bool face_orientation = cell->face_orientation(face_no),
1455  face_flip = cell->face_flip(face_no),
1456  face_rotation = cell->face_rotation(face_no);
1457  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1458  lines_per_face = GeometryInfo<3>::lines_per_face;
1459 
1460  // some sanity checks up front
1461  for (unsigned int i = 0; i < vertices_per_face; ++i)
1462  Assert(face->vertex_index(i) ==
1463  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1464  face_no, i, face_orientation, face_flip, face_rotation)),
1465  ExcInternalError());
1466 
1467  // indices of the lines that bound a face are given by GeometryInfo<3>::
1468  // face_to_cell_lines
1469  for (unsigned int i = 0; i < lines_per_face; ++i)
1470  Assert(face->line(i) ==
1472  face_no, i, face_orientation, face_flip, face_rotation)),
1473  ExcInternalError());
1474 #endif
1475  // extract the points surrounding a quad from the points
1476  // already computed. First get the 4 vertices and then the points on
1477  // the four lines
1478  boost::container::small_vector<Point<3>, 200> tmp_points(
1481  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1482  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1483  if (polynomial_degree > 1)
1484  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1485  ++line)
1486  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1487  tmp_points[4 + line * (polynomial_degree - 1) + i] =
1489  (polynomial_degree - 1) *
1490  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
1491  i];
1492 
1493  const std::size_t n_rows =
1495  a.resize(a.size() + n_rows);
1496  auto a_view = make_array_view(a.end() - n_rows, a.end());
1497  face->get_manifold().get_new_points(
1498  make_array_view(tmp_points.begin(), tmp_points.end()),
1500  a_view);
1501  }
1502 }
1503 
1504 
1505 
1506 template <>
1507 void
1510  std::vector<Point<3>> & a) const
1511 {
1512  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1513  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1514  vertices[i] = cell->vertex(i);
1515 
1516  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1518  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1519  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1520  {
1521  Point<2> point(line_support_points[q1 + 1][0],
1522  line_support_points[q2 + 1][0]);
1523  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1524  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1525  }
1526 
1527  const std::size_t n_rows = weights.size(0);
1528  a.resize(a.size() + n_rows);
1529  auto a_view = make_array_view(a.end() - n_rows, a.end());
1530  cell->get_manifold().get_new_points(
1531  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1532 }
1533 
1534 
1535 
1536 template <int dim, int spacedim>
1537 void
1540  std::vector<Point<spacedim>> &) const
1541 {
1542  Assert(false, ExcInternalError());
1543 }
1544 
1545 
1546 
1547 template <int dim, int spacedim>
1548 std::vector<Point<spacedim>>
1550  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1551 {
1552  // get the vertices first
1553  std::vector<Point<spacedim>> a;
1554  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1555  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1556  a.push_back(cell->vertex(i));
1557 
1558  if (this->polynomial_degree > 1)
1559  {
1560  // check if all entities have the same manifold id which is when we can
1561  // simply ask the manifold for all points. the transfinite manifold can
1562  // do the interpolation better than this class, so if we detect that we
1563  // do not have to change anything here
1564  Assert(dim <= 3, ExcImpossibleInDim(dim));
1565  bool all_manifold_ids_are_equal = (dim == spacedim);
1566  if (all_manifold_ids_are_equal &&
1568  &cell->get_manifold()) == nullptr)
1569  {
1570  for (auto f : GeometryInfo<dim>::face_indices())
1571  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1572  all_manifold_ids_are_equal = false;
1573 
1574  if (dim == 3)
1575  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1576  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1577  all_manifold_ids_are_equal = false;
1578  }
1579 
1580  if (all_manifold_ids_are_equal)
1581  {
1582  const std::size_t n_rows = support_point_weights_cell.size(0);
1583  a.resize(a.size() + n_rows);
1584  auto a_view = make_array_view(a.end() - n_rows, a.end());
1585  cell->get_manifold().get_new_points(make_array_view(a.begin(),
1586  a.end() - n_rows),
1588  a_view);
1589  }
1590  else
1591  switch (dim)
1592  {
1593  case 1:
1594  add_line_support_points(cell, a);
1595  break;
1596  case 2:
1597  // in 2d, add the points on the four bounding lines to the
1598  // exterior (outer) points
1599  add_line_support_points(cell, a);
1600 
1601  // then get the interior support points
1602  if (dim != spacedim)
1603  add_quad_support_points(cell, a);
1604  else
1605  {
1606  const std::size_t n_rows =
1608  a.resize(a.size() + n_rows);
1609  auto a_view = make_array_view(a.end() - n_rows, a.end());
1610  cell->get_manifold().get_new_points(
1611  make_array_view(a.begin(), a.end() - n_rows),
1613  a_view);
1614  }
1615  break;
1616 
1617  case 3:
1618  // in 3d also add the points located on the boundary faces
1619  add_line_support_points(cell, a);
1620  add_quad_support_points(cell, a);
1621 
1622  // then compute the interior points
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  cell->get_manifold().get_new_points(
1629  make_array_view(a.begin(), a.end() - n_rows),
1631  a_view);
1632  }
1633  break;
1634 
1635  default:
1636  Assert(false, ExcNotImplemented());
1637  break;
1638  }
1639  }
1640 
1641  return a;
1642 }
1643 
1644 
1645 
1646 template <int dim, int spacedim>
1649  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1650 {
1652 }
1653 
1654 
1655 
1656 template <int dim, int spacedim>
1657 bool
1659  const ReferenceCell &reference_cell) const
1660 {
1661  Assert(dim == reference_cell.get_dimension(),
1662  ExcMessage("The dimension of your mapping (" +
1663  Utilities::to_string(dim) +
1664  ") and the reference cell cell_type (" +
1665  Utilities::to_string(reference_cell.get_dimension()) +
1666  " ) do not agree."));
1667 
1668  return reference_cell.is_hyper_cube();
1669 }
1670 
1671 
1672 
1673 //--------------------------- Explicit instantiations -----------------------
1674 #include "mapping_q_generic.inst"
1675 
1676 
Transformed quadrature weights.
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
std::vector< Tensor< 2, dim > > shape_second_derivatives
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
void transform_fields(const ArrayView< const Tensor< rank, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim >> &output)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1630
bool is_hyper_cube() const
const unsigned int polynomial_degree
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1448
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)
unsigned int size() const
Definition: collection.h:109
const std::vector< Point< dim > > & get_points() const
inline ::Table< 2, double > compute_support_point_weights_cell(const unsigned int polynomial_degree)
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
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
void maybe_compute_q_points(const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< Point< spacedim >> &quadrature_points)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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
Volume element.
Definition: fe_dgq.h:110
void do_fill_fe_face_values(const ::MappingQGeneric< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Determinant of the Jacobian.
const Table< 2, double > support_point_weights_cell
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:1583
MappingKind
Definition: mapping.h:64
static DataSetDescriptor cell()
Definition: qprojector.h:563
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
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:2442
const unsigned int polynomial_degree
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:409
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
QGaussLobatto< 1 > line_support_points
T fixed_power(const T t)
Definition: utilities.h:1081
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
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
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
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
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< std::vector< Tensor< 1, spacedim > > > aux
#define Assert(cond, exc)
Definition: exceptions.h:1473
const std::vector< Point< dim > > unit_cell_support_points
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)
void reference_cell(const ReferenceCell &reference_cell, Triangulation< dim, spacedim > &tria)
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:303
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:697
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
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]
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
void transform_hessians(const ArrayView< const Tensor< 3, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim >> &output)
void 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:473
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
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
unsigned int size() const
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Point< dim, Number > compute(const Point< spacedim, Number > &p) const
Point< 2 > first
Definition: grid_out.cc:4587
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
unsigned int get_dimension() 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
std::vector< Point< spacedim > > quadrature_points
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:125
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
size_type size(const unsigned int i) const
const std::vector< Point< 1 > > line_support_points
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:1403
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={})
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:446
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
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
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
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
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
static ::ExceptionBase & ExcNotImplemented()
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:380
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
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:2417
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
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)
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:10330
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)
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()
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)