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