Reference documentation for deal.II version Git fe553f7db3 2020-12-03 08:11:48 +0100
\(\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 
35 #include <deal.II/grid/tria.h>
37 
39 
40 #include <boost/container/small_vector.hpp>
41 
42 #include <algorithm>
43 #include <array>
44 #include <cmath>
45 #include <memory>
46 #include <numeric>
47 
48 
50 
51 
52 template <int dim, int spacedim>
54  const unsigned int polynomial_degree)
55  : polynomial_degree(polynomial_degree)
56  , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
57  , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
58  , tensor_product_quadrature(false)
59 {}
60 
61 
62 
63 template <int dim, int spacedim>
64 std::size_t
66 {
67  return (
80 }
81 
82 
83 
84 template <int dim, int spacedim>
85 void
87  const UpdateFlags update_flags,
88  const Quadrature<dim> &q,
89  const unsigned int n_original_q_points)
90 {
91  // store the flags in the internal data object so we can access them
92  // in fill_fe_*_values()
93  this->update_each = update_flags;
94 
95  const unsigned int n_q_points = q.size();
96 
97  const bool needs_higher_order_terms =
98  this->update_each &
103 
105  covariant.resize(n_original_q_points);
106 
108  contravariant.resize(n_original_q_points);
109 
111  volume_elements.resize(n_original_q_points);
112 
114 
115  // use of MatrixFree only for higher order elements and with more than one
116  // point where tensor products do not make sense
117  if (polynomial_degree < 2 || n_q_points == 1)
119 
120  if (dim > 1)
121  {
122  // find out if the one-dimensional formula is the same
123  // in all directions
125  {
126  const std::array<Quadrature<1>, dim> quad_array =
127  q.get_tensor_basis();
128  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
129  {
130  if (quad_array[i - 1].size() != quad_array[i].size())
131  {
132  tensor_product_quadrature = false;
133  break;
134  }
135  else
136  {
137  const std::vector<Point<1>> &points_1 =
138  quad_array[i - 1].get_points();
139  const std::vector<Point<1>> &points_2 =
140  quad_array[i].get_points();
141  const std::vector<double> &weights_1 =
142  quad_array[i - 1].get_weights();
143  const std::vector<double> &weights_2 =
144  quad_array[i].get_weights();
145  for (unsigned int j = 0; j < quad_array[i].size(); ++j)
146  {
147  if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
148  std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
149  {
150  tensor_product_quadrature = false;
151  break;
152  }
153  }
154  }
155  }
156 
157  if (tensor_product_quadrature)
158  {
159  // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
160  // numbering manually (building an FE_Q<dim> is relatively
161  // expensive due to constraints)
162  const FE_DGQ<1> fe(polynomial_degree);
163  shape_info.reinit(q.get_tensor_basis()[0], fe);
165  FETools::lexicographic_to_hierarchic_numbering<dim>(
167  shape_info.n_q_points = q.size();
170  }
171  }
172  }
173 
174  // Only fill the big arrays on demand in case we cannot use the tensor
175  // product quadrature code path
176  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
177  {
178  // see if we need the (transformation) shape function values
179  // and/or gradients and resize the necessary arrays
181  shape_values.resize(n_shape_functions * n_q_points);
182 
183  if (this->update_each &
184  (update_covariant_transformation |
185  update_contravariant_transformation | update_JxW_values |
193  shape_derivatives.resize(n_shape_functions * n_q_points);
194 
195  if (this->update_each &
197  shape_second_derivatives.resize(n_shape_functions * n_q_points);
198 
201  shape_third_derivatives.resize(n_shape_functions * n_q_points);
202 
205  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
206 
207  // now also fill the various fields with their correct values
209  }
210 }
211 
212 
213 
214 template <int dim, int spacedim>
215 void
217  const UpdateFlags update_flags,
218  const Quadrature<dim> &q,
219  const unsigned int n_original_q_points)
220 {
221  initialize(update_flags, q, n_original_q_points);
222 
223  if (dim > 1 && tensor_product_quadrature)
224  {
225  constexpr unsigned int facedim = dim - 1;
226  const FE_DGQ<1> fe(polynomial_degree);
227  shape_info.reinit(q.get_tensor_basis()[0], fe);
229  FETools::lexicographic_to_hierarchic_numbering<facedim>(
231  shape_info.n_q_points = n_original_q_points;
234  }
235 
236  if (dim > 1)
237  {
238  if (this->update_each &
241  {
242  aux.resize(dim - 1,
243  std::vector<Tensor<1, spacedim>>(n_original_q_points));
244 
245  // Compute tangentials to the unit cell.
246  for (const unsigned int i : GeometryInfo<dim>::face_indices())
247  {
248  unit_tangentials[i].resize(n_original_q_points);
249  std::fill(unit_tangentials[i].begin(),
250  unit_tangentials[i].end(),
252  if (dim > 2)
253  {
255  .resize(n_original_q_points);
256  std::fill(
258  .begin(),
260  .end(),
262  }
263  }
264  }
265  }
266 }
267 
268 
269 
270 template <int dim, int spacedim>
271 void
273  const std::vector<Point<dim>> &unit_points)
274 {
275  const unsigned int n_points = unit_points.size();
276 
277  // Construct the tensor product polynomials used as shape functions for
278  // the Qp mapping of cells at the boundary.
279  const TensorProductPolynomials<dim> tensor_pols(
282  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
283 
284  // then also construct the mapping from lexicographic to the Qp shape
285  // function numbering
286  const std::vector<unsigned int> renumber =
287  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
288 
289  std::vector<double> values;
290  std::vector<Tensor<1, dim>> grads;
291  if (shape_values.size() != 0)
292  {
293  Assert(shape_values.size() == n_shape_functions * n_points,
294  ExcInternalError());
295  values.resize(n_shape_functions);
296  }
297  if (shape_derivatives.size() != 0)
298  {
299  Assert(shape_derivatives.size() == n_shape_functions * n_points,
300  ExcInternalError());
301  grads.resize(n_shape_functions);
302  }
303 
304  std::vector<Tensor<2, dim>> grad2;
305  if (shape_second_derivatives.size() != 0)
306  {
308  ExcInternalError());
309  grad2.resize(n_shape_functions);
310  }
311 
312  std::vector<Tensor<3, dim>> grad3;
313  if (shape_third_derivatives.size() != 0)
314  {
315  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
316  ExcInternalError());
317  grad3.resize(n_shape_functions);
318  }
319 
320  std::vector<Tensor<4, dim>> grad4;
321  if (shape_fourth_derivatives.size() != 0)
322  {
324  ExcInternalError());
325  grad4.resize(n_shape_functions);
326  }
327 
328 
329  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
330  shape_second_derivatives.size() != 0 ||
331  shape_third_derivatives.size() != 0 ||
332  shape_fourth_derivatives.size() != 0)
333  for (unsigned int point = 0; point < n_points; ++point)
334  {
335  tensor_pols.evaluate(
336  unit_points[point], values, grads, grad2, grad3, grad4);
337 
338  if (shape_values.size() != 0)
339  for (unsigned int i = 0; i < n_shape_functions; ++i)
340  shape(point, i) = values[renumber[i]];
341 
342  if (shape_derivatives.size() != 0)
343  for (unsigned int i = 0; i < n_shape_functions; ++i)
344  derivative(point, i) = grads[renumber[i]];
345 
346  if (shape_second_derivatives.size() != 0)
347  for (unsigned int i = 0; i < n_shape_functions; ++i)
348  second_derivative(point, i) = grad2[renumber[i]];
349 
350  if (shape_third_derivatives.size() != 0)
351  for (unsigned int i = 0; i < n_shape_functions; ++i)
352  third_derivative(point, i) = grad3[renumber[i]];
353 
354  if (shape_fourth_derivatives.size() != 0)
355  for (unsigned int i = 0; i < n_shape_functions; ++i)
356  fourth_derivative(point, i) = grad4[renumber[i]];
357  }
358 }
359 
360 
361 
362 template <int dim, int spacedim>
364  : polynomial_degree(p)
366  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
367  , polynomials_1d(
372  internal::MappingQGenericImplementation::unit_support_points<dim>(
376  internal::MappingQGenericImplementation::
378  this->polynomial_degree,
379  dim))
381  internal::MappingQGenericImplementation::
383 {
384  Assert(p >= 1,
385  ExcMessage("It only makes sense to create polynomial mappings "
386  "with a polynomial degree greater or equal to one."));
387 }
388 
389 
390 
391 template <int dim, int spacedim>
393  const MappingQGeneric<dim, spacedim> &mapping)
396  , polynomials_1d(mapping.polynomials_1d)
402 {}
403 
404 
405 
406 template <int dim, int spacedim>
407 std::unique_ptr<Mapping<dim, spacedim>>
409 {
410  return std::make_unique<MappingQGeneric<dim, spacedim>>(*this);
411 }
412 
413 
414 
415 template <int dim, int spacedim>
416 unsigned int
418 {
419  return polynomial_degree;
420 }
421 
422 
423 
424 template <int dim, int spacedim>
427  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
428  const Point<dim> & p) const
429 {
432  this->compute_mapping_support_points(cell),
433  p,
434  polynomials_1d.size() == 2,
436  .first);
437 }
438 
439 
440 // In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
441 // seeing which of the overloaded versions of
442 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
443 // error messages and, generally, nothing very good. Avoid this by ensuring
444 // that this class exists, but does not have an inner InternalData
445 // type, thereby ruling out the codim-1 version of the function
446 // below when doing overload resolution.
447 template <>
448 class MappingQGeneric<3, 4>
449 {};
450 
451 
452 
453 // visual studio freaks out when trying to determine if
454 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
455 // candidate. So instead of letting the compiler pick the correct overload, we
456 // use template specialization to make sure we pick up the right function to
457 // call:
458 
459 template <int dim, int spacedim>
463  const Point<spacedim> &,
464  const Point<dim> &) const
465 {
466  // default implementation (should never be called)
467  Assert(false, ExcInternalError());
468  return Point<dim>();
469 }
470 
471 
472 
473 template <>
474 Point<1>
477  const Point<1> & p,
478  const Point<1> & initial_p_unit) const
479 {
480  // dispatch to the various specializations for spacedim=dim,
481  // spacedim=dim+1, etc
482  return internal::MappingQGenericImplementation::
483  do_transform_real_to_unit_cell_internal<1>(
484  p,
485  initial_p_unit,
486  this->compute_mapping_support_points(cell),
489 }
490 
491 
492 
493 template <>
494 Point<2>
497  const Point<2> & p,
498  const Point<2> & initial_p_unit) const
499 {
500  return internal::MappingQGenericImplementation::
501  do_transform_real_to_unit_cell_internal<2>(
502  p,
503  initial_p_unit,
504  this->compute_mapping_support_points(cell),
507 }
508 
509 
510 
511 template <>
512 Point<3>
515  const Point<3> & p,
516  const Point<3> & initial_p_unit) const
517 {
518  return internal::MappingQGenericImplementation::
519  do_transform_real_to_unit_cell_internal<3>(
520  p,
521  initial_p_unit,
522  this->compute_mapping_support_points(cell),
525 }
526 
527 
528 
529 template <>
530 Point<1>
533  const Point<2> & p,
534  const Point<1> & initial_p_unit) const
535 {
536  const int dim = 1;
537  const int spacedim = 2;
538 
539  const Quadrature<dim> point_quadrature(initial_p_unit);
540 
542  if (spacedim > dim)
543  update_flags |= update_jacobian_grads;
545  get_data(update_flags, point_quadrature));
546 
548 
549  // dispatch to the various specializations for spacedim=dim,
550  // spacedim=dim+1, etc
551  return internal::MappingQGenericImplementation::
552  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
553  p,
554  initial_p_unit,
555  *mdata);
556 }
557 
558 
559 
560 template <>
561 Point<2>
564  const Point<3> & p,
565  const Point<2> & initial_p_unit) const
566 {
567  const int dim = 2;
568  const int spacedim = 3;
569 
570  const Quadrature<dim> point_quadrature(initial_p_unit);
571 
573  if (spacedim > dim)
574  update_flags |= update_jacobian_grads;
576  get_data(update_flags, point_quadrature));
577 
579 
580  // dispatch to the various specializations for spacedim=dim,
581  // spacedim=dim+1, etc
582  return internal::MappingQGenericImplementation::
583  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
584  p,
585  initial_p_unit,
586  *mdata);
587 }
588 
589 template <>
590 Point<1>
593  const Point<3> &,
594  const Point<1> &) const
595 {
596  Assert(false, ExcNotImplemented());
597  return {};
598 }
599 
600 
601 
602 template <int dim, int spacedim>
605  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
606  const Point<spacedim> & p) const
607 {
608  // Use an exact formula if one is available. this is only the case
609  // for Q1 mappings in 1d, and in 2d if dim==spacedim
610  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
611  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
612  {
613  // The dimension-dependent algorithms are much faster (about 25-45x in
614  // 2D) but fail most of the time when the given point (p) is not in the
615  // cell. The dimension-independent Newton algorithm given below is
616  // slower, but more robust (though it still sometimes fails). Therefore
617  // this function implements the following strategy based on the
618  // p's dimension:
619  //
620  // * In 1D this mapping is linear, so the mapping is always invertible
621  // (and the exact formula is known) as long as the cell has non-zero
622  // length.
623  // * In 2D the exact (quadratic) formula is called first. If either the
624  // exact formula does not succeed (negative discriminant in the
625  // quadratic formula) or succeeds but finds a solution outside of the
626  // unit cell, then the Newton solver is called. The rationale for the
627  // second choice is that the exact formula may provide two different
628  // answers when mapping a point outside of the real cell, but the
629  // Newton solver (if it converges) will only return one answer.
630  // Otherwise the exact formula successfully found a point in the unit
631  // cell and that value is returned.
632  // * In 3D there is no (known to the authors) exact formula, so the Newton
633  // algorithm is used.
634  const auto vertices_ = this->get_vertices(cell);
635 
636  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
637  vertices;
638  for (unsigned int i = 0; i < vertices.size(); ++i)
639  vertices[i] = vertices_[i];
640 
641  try
642  {
643  switch (dim)
644  {
645  case 1:
646  {
647  // formula not subject to any issues in 1d
648  if (spacedim == 1)
650  vertices, p);
651  else
652  break;
653  }
654 
655  case 2:
656  {
657  const Point<dim> point =
659  p);
660 
661  // formula not guaranteed to work for points outside of
662  // the cell. only take the computed point if it lies
663  // inside the reference cell
664  const double eps = 1e-15;
665  if (-eps <= point(1) && point(1) <= 1 + eps &&
666  -eps <= point(0) && point(0) <= 1 + eps)
667  {
668  return point;
669  }
670  else
671  break;
672  }
673 
674  default:
675  {
676  // we should get here, based on the if-condition at the top
677  Assert(false, ExcInternalError());
678  }
679  }
680  }
681  catch (
683  {
684  // simply fall through and continue on to the standard Newton code
685  }
686  }
687  else
688  {
689  // we can't use an explicit formula,
690  }
691 
692 
693  // Find the initial value for the Newton iteration by a normal
694  // projection to the least square plane determined by the vertices
695  // of the cell
696  Point<dim> initial_p_unit;
697  if (this->preserves_vertex_locations())
698  {
699  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
700  // in 1d with spacedim > 1 the affine approximation is exact
701  if (dim == 1 && polynomial_degree == 1)
702  return initial_p_unit;
703  }
704  else
705  {
706  // else, we simply use the mid point
707  for (unsigned int d = 0; d < dim; ++d)
708  initial_p_unit[d] = 0.5;
709  }
710 
711  // perform the Newton iteration and return the result. note that this
712  // statement may throw an exception, which we simply pass up to the caller
713  const Point<dim> p_unit =
714  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
715  if (p_unit[0] == std::numeric_limits<double>::infinity())
716  AssertThrow(false,
718  return p_unit;
719 }
720 
721 
722 
723 template <int dim, int spacedim>
724 void
726  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
727  const ArrayView<const Point<spacedim>> & real_points,
728  const ArrayView<Point<dim>> & unit_points) const
729 {
730  // Go to base class functions for dim < spacedim because it is not yet
731  // implemented with optimized code.
732  if (dim < spacedim)
733  {
735  real_points,
736  unit_points);
737  return;
738  }
739 
740  AssertDimension(real_points.size(), unit_points.size());
741  const std::vector<Point<spacedim>> support_points =
742  this->compute_mapping_support_points(cell);
743 
744  // From the given (high-order) support points, now only pick the first
745  // 2^dim points and construct an affine approximation from those.
746  internal::MappingQGenericImplementation::
747  InverseQuadraticApproximation<dim, spacedim>
748  inverse_approximation(support_points, unit_cell_support_points);
749 
750  const unsigned int n_points = real_points.size();
751  const unsigned int n_lanes = VectorizedArray<double>::size();
752 
753  // Use the more heavy VectorizedArray code path if there is more than
754  // one point left to compute
755  for (unsigned int i = 0; i < n_points; i += n_lanes)
756  if (n_points - i > 1)
757  {
759  for (unsigned int j = 0; j < n_lanes; ++j)
760  if (i + j < n_points)
761  for (unsigned int d = 0; d < spacedim; ++d)
762  p_vec[d][j] = real_points[i + j][d];
763  else
764  for (unsigned int d = 0; d < spacedim; ++d)
765  p_vec[d][j] = real_points[i][d];
766 
768  internal::MappingQGenericImplementation::
769  do_transform_real_to_unit_cell_internal<dim, spacedim>(
770  p_vec,
771  inverse_approximation.compute(p_vec),
772  support_points,
775 
776  // If the vectorized computation failed, it could be that only some of
777  // the lanes failed but others would have succeeded if we had let them
778  // compute alone without interference (like negative Jacobian
779  // determinants) from other SIMD lanes. Repeat the computation in this
780  // unlikely case with scalar arguments.
781  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
782  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
783  unit_points[i + j] = internal::MappingQGenericImplementation::
784  do_transform_real_to_unit_cell_internal<dim, spacedim>(
785  real_points[i + j],
786  inverse_approximation.compute(real_points[i + j]),
787  support_points,
790  else
791  for (unsigned int d = 0; d < dim; ++d)
792  unit_points[i + j][d] = unit_point[d][j];
793  }
794  else
795  unit_points[i] = internal::MappingQGenericImplementation::
796  do_transform_real_to_unit_cell_internal<dim, spacedim>(
797  real_points[i],
798  inverse_approximation.compute(real_points[i]),
799  support_points,
802 }
803 
804 
805 
806 template <int dim, int spacedim>
809  const UpdateFlags in) const
810 {
811  // add flags if the respective quantities are necessary to compute
812  // what we need. note that some flags appear in both the conditions
813  // and in subsequent set operations. this leads to some circular
814  // logic. the only way to treat this is to iterate. since there are
815  // 5 if-clauses in the loop, it will take at most 5 iterations to
816  // converge. do them:
817  UpdateFlags out = in;
818  for (unsigned int i = 0; i < 5; ++i)
819  {
820  // The following is a little incorrect:
821  // If not applied on a face,
822  // update_boundary_forms does not
823  // make sense. On the other hand,
824  // it is necessary on a
825  // face. Currently,
826  // update_boundary_forms is simply
827  // ignored for the interior of a
828  // cell.
830  out |= update_boundary_forms;
831 
836 
837  if (out &
842 
843  // The contravariant transformation is used in the Piola
844  // transformation, which requires the determinant of the Jacobi
845  // matrix of the transformation. Because we have no way of
846  // knowing here whether the finite element wants to use the
847  // contravariant or the Piola transforms, we add the JxW values
848  // to the list of flags to be updated for each cell.
850  out |= update_volume_elements;
851 
852  // the same is true when computing normal vectors: they require
853  // the determinant of the Jacobian
854  if (out & update_normal_vectors)
855  out |= update_volume_elements;
856  }
857 
858  return out;
859 }
860 
861 
862 
863 template <int dim, int spacedim>
864 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
866  const Quadrature<dim> &q) const
867 {
868  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
869  std::make_unique<InternalData>(polynomial_degree);
870  auto &data = dynamic_cast<InternalData &>(*data_ptr);
871  data.initialize(this->requires_update_flags(update_flags), q, q.size());
872 
873  return data_ptr;
874 }
875 
876 
877 
878 template <int dim, int spacedim>
879 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
881  const UpdateFlags update_flags,
882  const Quadrature<dim - 1> &quadrature) const
883 {
884  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
885  std::make_unique<InternalData>(polynomial_degree);
886  auto &data = dynamic_cast<InternalData &>(*data_ptr);
887  data.initialize_face(this->requires_update_flags(update_flags),
889  ReferenceCell::get_hypercube(dim), quadrature),
890  quadrature.size());
891 
892  return data_ptr;
893 }
894 
895 
896 
897 template <int dim, int spacedim>
898 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
900  const UpdateFlags update_flags,
901  const Quadrature<dim - 1> &quadrature) const
902 {
903  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
904  std::make_unique<InternalData>(polynomial_degree);
905  auto &data = dynamic_cast<InternalData &>(*data_ptr);
906  data.initialize_face(this->requires_update_flags(update_flags),
908  ReferenceCell::get_hypercube(dim), quadrature),
909  quadrature.size());
910 
911  return data_ptr;
912 }
913 
914 
915 
916 template <int dim, int spacedim>
919  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
920  const CellSimilarity::Similarity cell_similarity,
921  const Quadrature<dim> & quadrature,
922  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
924  &output_data) const
925 {
926  // ensure that the following static_cast is really correct:
927  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
928  ExcInternalError());
929  const InternalData &data = static_cast<const InternalData &>(internal_data);
930 
931  const unsigned int n_q_points = quadrature.size();
932 
933  // recompute the support points of the transformation of this
934  // cell. we tried to be clever here in an earlier version of the
935  // library by checking whether the cell is the same as the one we
936  // had visited last, but it turns out to be difficult to determine
937  // that because a cell for the purposes of a mapping is
938  // characterized not just by its (triangulation, level, index)
939  // triple, but also by the locations of its vertices, the manifold
940  // object attached to the cell and all of its bounding faces/edges,
941  // etc. to reliably test that the "cell" we are on is, therefore,
942  // not easily done
944  data.cell_of_current_support_points = cell;
945 
946  // if the order of the mapping is greater than 1, then do not reuse any cell
947  // similarity information. This is necessary because the cell similarity
948  // value is computed with just cell vertices and does not take into account
949  // cell curvature.
950  const CellSimilarity::Similarity computed_cell_similarity =
951  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
952 
953  if (dim > 1 && data.tensor_product_quadrature)
954  {
955  internal::MappingQGenericImplementation::
956  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
957  computed_cell_similarity,
958  data,
959  output_data.quadrature_points,
960  output_data.jacobian_grads);
961  }
962  else
963  {
965  spacedim>(
967  data,
968  output_data.quadrature_points);
969 
971  spacedim>(
972  computed_cell_similarity,
974  data);
975 
977  dim,
978  spacedim>(computed_cell_similarity,
980  data,
981  output_data.jacobian_grads);
982  }
983 
984  internal::MappingQGenericImplementation::
985  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
986  computed_cell_similarity,
988  data,
989  output_data.jacobian_pushed_forward_grads);
990 
991  internal::MappingQGenericImplementation::
992  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
993  computed_cell_similarity,
995  data,
996  output_data.jacobian_2nd_derivatives);
997 
998  internal::MappingQGenericImplementation::
999  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1000  computed_cell_similarity,
1002  data,
1004 
1005  internal::MappingQGenericImplementation::
1006  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1007  computed_cell_similarity,
1009  data,
1010  output_data.jacobian_3rd_derivatives);
1011 
1012  internal::MappingQGenericImplementation::
1013  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1014  computed_cell_similarity,
1016  data,
1018 
1019  const UpdateFlags update_flags = data.update_each;
1020  const std::vector<double> &weights = quadrature.get_weights();
1021 
1022  // Multiply quadrature weights by absolute value of Jacobian determinants or
1023  // the area element g=sqrt(DX^t DX) in case of codim > 0
1024 
1025  if (update_flags & (update_normal_vectors | update_JxW_values))
1026  {
1027  AssertDimension(output_data.JxW_values.size(), n_q_points);
1028 
1029  Assert(!(update_flags & update_normal_vectors) ||
1030  (output_data.normal_vectors.size() == n_q_points),
1031  ExcDimensionMismatch(output_data.normal_vectors.size(),
1032  n_q_points));
1033 
1034 
1035  if (computed_cell_similarity != CellSimilarity::translation)
1036  for (unsigned int point = 0; point < n_q_points; ++point)
1037  {
1038  if (dim == spacedim)
1039  {
1040  const double det = data.contravariant[point].determinant();
1041 
1042  // check for distorted cells.
1043 
1044  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1045  // 1e12 in 2D. might want to find a finer
1046  // (dimension-independent) criterion
1047  Assert(det >
1048  1e-12 * Utilities::fixed_power<dim>(
1049  cell->diameter() / std::sqrt(double(dim))),
1051  cell->center(), det, point)));
1052 
1053  output_data.JxW_values[point] = weights[point] * det;
1054  }
1055  // if dim==spacedim, then there is no cell normal to
1056  // compute. since this is for FEValues (and not FEFaceValues),
1057  // there are also no face normals to compute
1058  else // codim>0 case
1059  {
1060  Tensor<1, spacedim> DX_t[dim];
1061  for (unsigned int i = 0; i < spacedim; ++i)
1062  for (unsigned int j = 0; j < dim; ++j)
1063  DX_t[j][i] = data.contravariant[point][i][j];
1064 
1065  Tensor<2, dim> G; // First fundamental form
1066  for (unsigned int i = 0; i < dim; ++i)
1067  for (unsigned int j = 0; j < dim; ++j)
1068  G[i][j] = DX_t[i] * DX_t[j];
1069 
1070  output_data.JxW_values[point] =
1071  std::sqrt(determinant(G)) * weights[point];
1072 
1073  if (computed_cell_similarity ==
1075  {
1076  // we only need to flip the normal
1077  if (update_flags & update_normal_vectors)
1078  output_data.normal_vectors[point] *= -1.;
1079  }
1080  else
1081  {
1082  if (update_flags & update_normal_vectors)
1083  {
1084  Assert(spacedim == dim + 1,
1085  ExcMessage(
1086  "There is no (unique) cell normal for " +
1088  "-dimensional cells in " +
1089  Utilities::int_to_string(spacedim) +
1090  "-dimensional space. This only works if the "
1091  "space dimension is one greater than the "
1092  "dimensionality of the mesh cells."));
1093 
1094  if (dim == 1)
1095  output_data.normal_vectors[point] =
1096  cross_product_2d(-DX_t[0]);
1097  else // dim == 2
1098  output_data.normal_vectors[point] =
1099  cross_product_3d(DX_t[0], DX_t[1]);
1100 
1101  output_data.normal_vectors[point] /=
1102  output_data.normal_vectors[point].norm();
1103 
1104  if (cell->direction_flag() == false)
1105  output_data.normal_vectors[point] *= -1.;
1106  }
1107  }
1108  } // codim>0 case
1109  }
1110  }
1111 
1112 
1113 
1114  // copy values from InternalData to vector given by reference
1115  if (update_flags & update_jacobians)
1116  {
1117  AssertDimension(output_data.jacobians.size(), n_q_points);
1118  if (computed_cell_similarity != CellSimilarity::translation)
1119  for (unsigned int point = 0; point < n_q_points; ++point)
1120  output_data.jacobians[point] = data.contravariant[point];
1121  }
1122 
1123  // copy values from InternalData to vector given by reference
1124  if (update_flags & update_inverse_jacobians)
1125  {
1126  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1127  if (computed_cell_similarity != CellSimilarity::translation)
1128  for (unsigned int point = 0; point < n_q_points; ++point)
1129  output_data.inverse_jacobians[point] =
1130  data.covariant[point].transpose();
1131  }
1132 
1133  return computed_cell_similarity;
1134 }
1135 
1136 
1137 
1138 template <int dim, int spacedim>
1139 void
1141  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1142  const unsigned int face_no,
1143  const Quadrature<dim - 1> & quadrature,
1144  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1146  &output_data) const
1147 {
1148  // ensure that the following cast is really correct:
1149  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1150  ExcInternalError());
1151  const InternalData &data = static_cast<const InternalData &>(internal_data);
1152 
1153  // if necessary, recompute the support points of the transformation of this
1154  // cell (note that we need to first check the triangulation pointer, since
1155  // otherwise the second test might trigger an exception if the triangulations
1156  // are not the same)
1157  if ((data.mapping_support_points.size() == 0) ||
1158  (&cell->get_triangulation() !=
1159  &data.cell_of_current_support_points->get_triangulation()) ||
1160  (cell != data.cell_of_current_support_points))
1161  {
1163  data.cell_of_current_support_points = cell;
1164  }
1165 
1167  *this,
1168  cell,
1169  face_no,
1172  face_no,
1173  cell->face_orientation(face_no),
1174  cell->face_flip(face_no),
1175  cell->face_rotation(face_no),
1176  quadrature.size()),
1177  quadrature,
1178  data,
1179  output_data);
1180 }
1181 
1182 
1183 
1184 template <int dim, int spacedim>
1185 void
1187  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1188  const unsigned int face_no,
1189  const unsigned int subface_no,
1190  const Quadrature<dim - 1> & quadrature,
1191  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1193  &output_data) const
1194 {
1195  // ensure that the following cast is really correct:
1196  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1197  ExcInternalError());
1198  const InternalData &data = static_cast<const InternalData &>(internal_data);
1199 
1200  // if necessary, recompute the support points of the transformation of this
1201  // cell (note that we need to first check the triangulation pointer, since
1202  // otherwise the second test might trigger an exception if the triangulations
1203  // are not the same)
1204  if ((data.mapping_support_points.size() == 0) ||
1205  (&cell->get_triangulation() !=
1206  &data.cell_of_current_support_points->get_triangulation()) ||
1207  (cell != data.cell_of_current_support_points))
1208  {
1210  data.cell_of_current_support_points = cell;
1211  }
1212 
1214  *this,
1215  cell,
1216  face_no,
1217  subface_no,
1219  dim),
1220  face_no,
1221  subface_no,
1222  cell->face_orientation(face_no),
1223  cell->face_flip(face_no),
1224  cell->face_rotation(face_no),
1225  quadrature.size(),
1226  cell->subface_case(face_no)),
1227  quadrature,
1228  data,
1229  output_data);
1230 }
1231 
1232 
1233 
1234 template <int dim, int spacedim>
1235 void
1237  const ArrayView<const Tensor<1, dim>> & input,
1238  const MappingKind mapping_kind,
1239  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1240  const ArrayView<Tensor<1, spacedim>> & output) const
1241 {
1243  mapping_kind,
1244  mapping_data,
1245  output);
1246 }
1247 
1248 
1249 
1250 template <int dim, int spacedim>
1251 void
1253  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1254  const MappingKind mapping_kind,
1255  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1256  const ArrayView<Tensor<2, spacedim>> & output) const
1257 {
1259  input, mapping_kind, mapping_data, output);
1260 }
1261 
1262 
1263 
1264 template <int dim, int spacedim>
1265 void
1267  const ArrayView<const Tensor<2, dim>> & input,
1268  const MappingKind mapping_kind,
1269  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1270  const ArrayView<Tensor<2, spacedim>> & output) const
1271 {
1272  switch (mapping_kind)
1273  {
1274  case mapping_contravariant:
1276  mapping_kind,
1277  mapping_data,
1278  output);
1279  return;
1280 
1285  input, mapping_kind, mapping_data, output);
1286  return;
1287  default:
1288  Assert(false, ExcNotImplemented());
1289  }
1290 }
1291 
1292 
1293 
1294 template <int dim, int spacedim>
1295 void
1297  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1298  const MappingKind mapping_kind,
1299  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1300  const ArrayView<Tensor<3, spacedim>> & output) const
1301 {
1302  AssertDimension(input.size(), output.size());
1303  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1304  ExcInternalError());
1305  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1306 
1307  switch (mapping_kind)
1308  {
1310  {
1313  "update_covariant_transformation"));
1314 
1315  for (unsigned int q = 0; q < output.size(); ++q)
1316  for (unsigned int i = 0; i < spacedim; ++i)
1317  for (unsigned int j = 0; j < spacedim; ++j)
1318  {
1319  double tmp[dim];
1320  for (unsigned int K = 0; K < dim; ++K)
1321  {
1322  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1323  for (unsigned int J = 1; J < dim; ++J)
1324  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1325  }
1326  for (unsigned int k = 0; k < spacedim; ++k)
1327  {
1328  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1329  for (unsigned int K = 1; K < dim; ++K)
1330  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1331  }
1332  }
1333  return;
1334  }
1335 
1336  default:
1337  Assert(false, ExcNotImplemented());
1338  }
1339 }
1340 
1341 
1342 
1343 template <int dim, int spacedim>
1344 void
1346  const ArrayView<const Tensor<3, dim>> & input,
1347  const MappingKind mapping_kind,
1348  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1349  const ArrayView<Tensor<3, spacedim>> & output) const
1350 {
1351  switch (mapping_kind)
1352  {
1353  case mapping_piola_hessian:
1357  input, mapping_kind, mapping_data, output);
1358  return;
1359  default:
1360  Assert(false, ExcNotImplemented());
1361  }
1362 }
1363 
1364 
1365 
1366 template <int dim, int spacedim>
1367 void
1369  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1370  std::vector<Point<spacedim>> & a) const
1371 {
1372  // if we only need the midpoint, then ask for it.
1373  if (this->polynomial_degree == 2)
1374  {
1375  for (unsigned int line_no = 0;
1376  line_no < GeometryInfo<dim>::lines_per_cell;
1377  ++line_no)
1378  {
1379  const typename Triangulation<dim, spacedim>::line_iterator line =
1380  (dim == 1 ?
1381  static_cast<
1383  cell->line(line_no));
1384 
1385  const Manifold<dim, spacedim> &manifold =
1386  ((line->manifold_id() == numbers::flat_manifold_id) &&
1387  (dim < spacedim) ?
1388  cell->get_manifold() :
1389  line->get_manifold());
1390  a.push_back(manifold.get_new_point_on_line(line));
1391  }
1392  }
1393  else
1394  // otherwise call the more complicated functions and ask for inner points
1395  // from the manifold description
1396  {
1397  std::vector<Point<spacedim>> tmp_points;
1398  for (unsigned int line_no = 0;
1399  line_no < GeometryInfo<dim>::lines_per_cell;
1400  ++line_no)
1401  {
1402  const typename Triangulation<dim, spacedim>::line_iterator line =
1403  (dim == 1 ?
1404  static_cast<
1406  cell->line(line_no));
1407 
1408  const Manifold<dim, spacedim> &manifold =
1409  ((line->manifold_id() == numbers::flat_manifold_id) &&
1410  (dim < spacedim) ?
1411  cell->get_manifold() :
1412  line->get_manifold());
1413 
1414  const std::array<Point<spacedim>, 2> vertices{
1415  {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
1416  cell->vertex(
1418 
1419  const std::size_t n_rows =
1421  a.resize(a.size() + n_rows);
1422  auto a_view = make_array_view(a.end() - n_rows, a.end());
1423  manifold.get_new_points(
1424  make_array_view(vertices.begin(), vertices.end()),
1426  a_view);
1427  }
1428  }
1429 }
1430 
1431 
1432 
1433 template <>
1434 void
1437  std::vector<Point<3>> & a) const
1438 {
1439  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1440 
1441  // used if face quad at boundary or entirely in the interior of the domain
1442  std::vector<Point<3>> tmp_points;
1443 
1444  // loop over all faces and collect points on them
1445  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1446  {
1447  const Triangulation<3>::face_iterator face = cell->face(face_no);
1448 
1449 #ifdef DEBUG
1450  const bool face_orientation = cell->face_orientation(face_no),
1451  face_flip = cell->face_flip(face_no),
1452  face_rotation = cell->face_rotation(face_no);
1453  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1454  lines_per_face = GeometryInfo<3>::lines_per_face;
1455 
1456  // some sanity checks up front
1457  for (unsigned int i = 0; i < vertices_per_face; ++i)
1458  Assert(face->vertex_index(i) ==
1459  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1460  face_no, i, face_orientation, face_flip, face_rotation)),
1461  ExcInternalError());
1462 
1463  // indices of the lines that bound a face are given by GeometryInfo<3>::
1464  // face_to_cell_lines
1465  for (unsigned int i = 0; i < lines_per_face; ++i)
1466  Assert(face->line(i) ==
1468  face_no, i, face_orientation, face_flip, face_rotation)),
1469  ExcInternalError());
1470 #endif
1471  // extract the points surrounding a quad from the points
1472  // already computed. First get the 4 vertices and then the points on
1473  // the four lines
1474  boost::container::small_vector<Point<3>, 200> tmp_points(
1477  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1478  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1479  if (polynomial_degree > 1)
1480  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1481  ++line)
1482  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1483  tmp_points[4 + line * (polynomial_degree - 1) + i] =
1485  (polynomial_degree - 1) *
1486  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
1487  i];
1488 
1489  const std::size_t n_rows =
1491  a.resize(a.size() + n_rows);
1492  auto a_view = make_array_view(a.end() - n_rows, a.end());
1493  face->get_manifold().get_new_points(
1494  make_array_view(tmp_points.begin(), tmp_points.end()),
1496  a_view);
1497  }
1498 }
1499 
1500 
1501 
1502 template <>
1503 void
1506  std::vector<Point<3>> & a) const
1507 {
1508  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1509  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1510  vertices[i] = cell->vertex(i);
1511 
1512  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1514  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1515  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1516  {
1517  Point<2> point(line_support_points[q1 + 1][0],
1518  line_support_points[q2 + 1][0]);
1519  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1520  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1521  }
1522 
1523  const std::size_t n_rows = weights.size(0);
1524  a.resize(a.size() + n_rows);
1525  auto a_view = make_array_view(a.end() - n_rows, a.end());
1526  cell->get_manifold().get_new_points(
1527  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1528 }
1529 
1530 
1531 
1532 template <int dim, int spacedim>
1533 void
1536  std::vector<Point<spacedim>> &) const
1537 {
1538  Assert(false, ExcInternalError());
1539 }
1540 
1541 
1542 
1543 template <int dim, int spacedim>
1544 std::vector<Point<spacedim>>
1546  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1547 {
1548  // get the vertices first
1549  std::vector<Point<spacedim>> a;
1550  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1551  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1552  a.push_back(cell->vertex(i));
1553 
1554  if (this->polynomial_degree > 1)
1555  {
1556  // check if all entities have the same manifold id which is when we can
1557  // simply ask the manifold for all points. the transfinite manifold can
1558  // do the interpolation better than this class, so if we detect that we
1559  // do not have to change anything here
1560  Assert(dim <= 3, ExcImpossibleInDim(dim));
1561  bool all_manifold_ids_are_equal = (dim == spacedim);
1562  if (all_manifold_ids_are_equal &&
1564  &cell->get_manifold()) == nullptr)
1565  {
1566  for (auto f : GeometryInfo<dim>::face_indices())
1567  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1568  all_manifold_ids_are_equal = false;
1569 
1570  if (dim == 3)
1571  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1572  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1573  all_manifold_ids_are_equal = false;
1574  }
1575 
1576  if (all_manifold_ids_are_equal)
1577  {
1578  const std::size_t n_rows = support_point_weights_cell.size(0);
1579  a.resize(a.size() + n_rows);
1580  auto a_view = make_array_view(a.end() - n_rows, a.end());
1581  cell->get_manifold().get_new_points(make_array_view(a.begin(),
1582  a.end() - n_rows),
1584  a_view);
1585  }
1586  else
1587  switch (dim)
1588  {
1589  case 1:
1590  add_line_support_points(cell, a);
1591  break;
1592  case 2:
1593  // in 2d, add the points on the four bounding lines to the
1594  // exterior (outer) points
1595  add_line_support_points(cell, a);
1596 
1597  // then get the interior support points
1598  if (dim != spacedim)
1599  add_quad_support_points(cell, a);
1600  else
1601  {
1602  const std::size_t n_rows =
1604  a.resize(a.size() + n_rows);
1605  auto a_view = make_array_view(a.end() - n_rows, a.end());
1606  cell->get_manifold().get_new_points(
1607  make_array_view(a.begin(), a.end() - n_rows),
1609  a_view);
1610  }
1611  break;
1612 
1613  case 3:
1614  // in 3d also add the points located on the boundary faces
1615  add_line_support_points(cell, a);
1616  add_quad_support_points(cell, a);
1617 
1618  // then compute the interior points
1619  {
1620  const std::size_t n_rows =
1622  a.resize(a.size() + n_rows);
1623  auto a_view = make_array_view(a.end() - n_rows, a.end());
1624  cell->get_manifold().get_new_points(
1625  make_array_view(a.begin(), a.end() - n_rows),
1627  a_view);
1628  }
1629  break;
1630 
1631  default:
1632  Assert(false, ExcNotImplemented());
1633  break;
1634  }
1635  }
1636 
1637  return a;
1638 }
1639 
1640 
1641 
1642 //--------------------------- Explicit instantiations -----------------------
1643 #include "mapping_q_generic.inst"
1644 
1645 
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:1623
const unsigned int polynomial_degree
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1444
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)
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:83
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)
Type get_hypercube(const unsigned int dim)
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
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:1576
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:564
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
InternalData(const unsigned int polynomial_degree)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2441
const unsigned int polynomial_degree
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:752
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQGeneric< dim, spacedim >::InternalData &data)
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
QGaussLobatto< 1 > line_support_points
T fixed_power(const T t)
Definition: utilities.h:1045
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:316
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
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:1466
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)
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
Abstract base class for mapping classes.
Definition: mapping.h:301
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:665
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
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:474
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int n_shape_functions
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
std::vector< Tensor< 1, dim > > shape_derivatives
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:4340
Point< dim > transform_real_to_unit_cell_internal(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit) const
std::vector< Point< spacedim > > quadrature_points
unsigned int get_degree() const
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
Definition: manifold.cc:123
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: cuda.h:32
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
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
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_NAMESPACE_OPEN
Definition: config.h:371
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:29
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
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2416
std::vector< double > shape_values
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
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:9685
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701
Covariant transformation.
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
std::vector< Tensor< 1, spacedim > > normal_vectors
internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< double > > shape_info
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
UpdateFlags update_each
Definition: mapping.h:643
static ::ExceptionBase & ExcInternalError()
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)