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