Reference documentation for deal.II version GIT 9d46a3bd8c 2023-10-03 12:10:02+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 - 2023 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;
174  const FE_DGQ<1> fe(polynomial_degree);
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())
224  , polynomials_1d(
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>
257  : polynomial_degree(mapping.polynomial_degree)
258  , line_support_points(mapping.line_support_points)
259  , polynomials_1d(mapping.polynomials_1d)
260  , renumber_lexicographic_to_hierarchic(
261  mapping.renumber_lexicographic_to_hierarchic)
262  , unit_cell_support_points(mapping.unit_cell_support_points)
263  , support_point_weights_perimeter_to_interior(
264  mapping.support_point_weights_perimeter_to_interior)
265  , support_point_weights_cell(mapping.support_point_weights_cell)
266 {}
267 
268 
269 
270 template <int dim, int spacedim>
271 std::unique_ptr<Mapping<dim, spacedim>>
273 {
274  return std::make_unique<MappingQ<dim, spacedim>>(*this);
275 }
276 
277 
278 
279 template <int dim, int spacedim>
280 unsigned int
282 {
283  return polynomial_degree;
284 }
285 
286 
287 
288 template <int dim, int spacedim>
291  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
292  const Point<dim> &p) const
293 {
295  polynomials_1d,
296  this->compute_mapping_support_points(cell),
297  p,
298  polynomials_1d.size() == 2,
299  renumber_lexicographic_to_hierarchic));
300 }
301 
302 
303 // In the code below, GCC tries to instantiate MappingQ<3,4> when
304 // seeing which of the overloaded versions of
305 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
306 // error messages and, generally, nothing very good. Avoid this by ensuring
307 // that this class exists, but does not have an inner InternalData
308 // type, thereby ruling out the codim-1 version of the function
309 // below when doing overload resolution.
310 template <>
311 class MappingQ<3, 4>
312 {};
313 
314 
315 
316 // visual studio freaks out when trying to determine if
317 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
318 // candidate. So instead of letting the compiler pick the correct overload, we
319 // use template specialization to make sure we pick up the right function to
320 // call:
321 
322 template <int dim, int spacedim>
326  const Point<spacedim> &,
327  const Point<dim> &) const
328 {
329  // default implementation (should never be called)
330  Assert(false, ExcInternalError());
331  return {};
332 }
333 
334 
335 
336 template <>
337 Point<1>
340  const Point<1> &p,
341  const Point<1> &initial_p_unit) const
342 {
343  // dispatch to the various specializations for spacedim=dim,
344  // spacedim=dim+1, etc
345  return internal::MappingQImplementation::
346  do_transform_real_to_unit_cell_internal<1>(
347  p,
348  initial_p_unit,
349  this->compute_mapping_support_points(cell),
350  polynomials_1d,
351  renumber_lexicographic_to_hierarchic);
352 }
353 
354 
355 
356 template <>
357 Point<2>
360  const Point<2> &p,
361  const Point<2> &initial_p_unit) const
362 {
363  return internal::MappingQImplementation::
364  do_transform_real_to_unit_cell_internal<2>(
365  p,
366  initial_p_unit,
367  this->compute_mapping_support_points(cell),
368  polynomials_1d,
369  renumber_lexicographic_to_hierarchic);
370 }
371 
372 
373 
374 template <>
375 Point<3>
378  const Point<3> &p,
379  const Point<3> &initial_p_unit) const
380 {
381  return internal::MappingQImplementation::
382  do_transform_real_to_unit_cell_internal<3>(
383  p,
384  initial_p_unit,
385  this->compute_mapping_support_points(cell),
386  polynomials_1d,
387  renumber_lexicographic_to_hierarchic);
388 }
389 
390 
391 
392 template <>
393 Point<1>
396  const Point<2> &p,
397  const Point<1> &initial_p_unit) const
398 {
399  const int dim = 1;
400  const int spacedim = 2;
401 
402  const Quadrature<dim> point_quadrature(initial_p_unit);
403 
405  if (spacedim > dim)
406  update_flags |= update_jacobian_grads;
407  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
408  get_data(update_flags, point_quadrature));
409 
410  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
411 
412  // dispatch to the various specializations for spacedim=dim,
413  // spacedim=dim+1, etc
414  return internal::MappingQImplementation::
415  do_transform_real_to_unit_cell_internal_codim1<1>(
416  p,
417  initial_p_unit,
418  mdata->mapping_support_points,
419  polynomials_1d,
420  renumber_lexicographic_to_hierarchic);
421 }
422 
423 
424 
425 template <>
426 Point<2>
429  const Point<3> &p,
430  const Point<2> &initial_p_unit) const
431 {
432  const int dim = 2;
433  const int spacedim = 3;
434 
435  const Quadrature<dim> point_quadrature(initial_p_unit);
436 
438  if (spacedim > dim)
439  update_flags |= update_jacobian_grads;
440  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
441  get_data(update_flags, point_quadrature));
442 
443  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
444 
445  // dispatch to the various specializations for spacedim=dim,
446  // spacedim=dim+1, etc
447  return internal::MappingQImplementation::
448  do_transform_real_to_unit_cell_internal_codim1<2>(
449  p,
450  initial_p_unit,
451  mdata->mapping_support_points,
452  polynomials_1d,
453  renumber_lexicographic_to_hierarchic);
454 }
455 
456 
457 
458 template <>
459 Point<1>
462  const Point<3> &,
463  const Point<1> &) const
464 {
466  return {};
467 }
468 
469 
470 
471 template <int dim, int spacedim>
474  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
475  const Point<spacedim> &p) const
476 {
477  // Use an exact formula if one is available. this is only the case
478  // for Q1 mappings in 1d, and in 2d if dim==spacedim
479  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
480  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
481  {
482  // The dimension-dependent algorithms are much faster (about 25-45x in
483  // 2d) but fail most of the time when the given point (p) is not in the
484  // cell. The dimension-independent Newton algorithm given below is
485  // slower, but more robust (though it still sometimes fails). Therefore
486  // this function implements the following strategy based on the
487  // p's dimension:
488  //
489  // * In 1d this mapping is linear, so the mapping is always invertible
490  // (and the exact formula is known) as long as the cell has non-zero
491  // length.
492  // * In 2d the exact (quadratic) formula is called first. If either the
493  // exact formula does not succeed (negative discriminant in the
494  // quadratic formula) or succeeds but finds a solution outside of the
495  // unit cell, then the Newton solver is called. The rationale for the
496  // second choice is that the exact formula may provide two different
497  // answers when mapping a point outside of the real cell, but the
498  // Newton solver (if it converges) will only return one answer.
499  // Otherwise the exact formula successfully found a point in the unit
500  // cell and that value is returned.
501  // * In 3d there is no (known to the authors) exact formula, so the Newton
502  // algorithm is used.
503  const auto vertices_ = this->get_vertices(cell);
504 
505  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
506  vertices;
507  for (unsigned int i = 0; i < vertices.size(); ++i)
508  vertices[i] = vertices_[i];
509 
510  try
511  {
512  switch (dim)
513  {
514  case 1:
515  {
516  // formula not subject to any issues in 1d
517  if (spacedim == 1)
519  vertices, p);
520  else
521  break;
522  }
523 
524  case 2:
525  {
526  const Point<dim> point =
528  p);
529 
530  // formula not guaranteed to work for points outside of
531  // the cell. only take the computed point if it lies
532  // inside the reference cell
533  const double eps = 1e-15;
534  if (-eps <= point(1) && point(1) <= 1 + eps &&
535  -eps <= point(0) && point(0) <= 1 + eps)
536  {
537  return point;
538  }
539  else
540  break;
541  }
542 
543  default:
544  {
545  // we should get here, based on the if-condition at the top
546  Assert(false, ExcInternalError());
547  }
548  }
549  }
550  catch (
552  {
553  // simply fall through and continue on to the standard Newton code
554  }
555  }
556  else
557  {
558  // we can't use an explicit formula,
559  }
560 
561 
562  // Find the initial value for the Newton iteration by a normal
563  // projection to the least square plane determined by the vertices
564  // of the cell
565  Point<dim> initial_p_unit;
566  if (this->preserves_vertex_locations())
567  {
568  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
569  // in 1d with spacedim > 1 the affine approximation is exact
570  if (dim == 1 && polynomial_degree == 1)
571  return initial_p_unit;
572  }
573  else
574  {
575  // else, we simply use the mid point
576  for (unsigned int d = 0; d < dim; ++d)
577  initial_p_unit[d] = 0.5;
578  }
579 
580  // perform the Newton iteration and return the result. note that this
581  // statement may throw an exception, which we simply pass up to the caller
582  const Point<dim> p_unit =
583  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
584  AssertThrow(numbers::is_finite(p_unit[0]),
586  return p_unit;
587 }
588 
589 
590 
591 template <int dim, int spacedim>
592 void
594  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
595  const ArrayView<const Point<spacedim>> &real_points,
596  const ArrayView<Point<dim>> &unit_points) const
597 {
598  // Go to base class functions for dim < spacedim because it is not yet
599  // implemented with optimized code.
600  if (dim < spacedim)
601  {
603  real_points,
604  unit_points);
605  return;
606  }
607 
608  AssertDimension(real_points.size(), unit_points.size());
609  const std::vector<Point<spacedim>> support_points =
610  this->compute_mapping_support_points(cell);
611 
612  // From the given (high-order) support points, now only pick the first
613  // 2^dim points and construct an affine approximation from those.
615  inverse_approximation(support_points, unit_cell_support_points);
616 
617  const unsigned int n_points = real_points.size();
618  const unsigned int n_lanes = VectorizedArray<double>::size();
619 
620  // Use the more heavy VectorizedArray code path if there is more than
621  // one point left to compute
622  for (unsigned int i = 0; i < n_points; i += n_lanes)
623  if (n_points - i > 1)
624  {
626  for (unsigned int j = 0; j < n_lanes; ++j)
627  if (i + j < n_points)
628  for (unsigned int d = 0; d < spacedim; ++d)
629  p_vec[d][j] = real_points[i + j][d];
630  else
631  for (unsigned int d = 0; d < spacedim; ++d)
632  p_vec[d][j] = real_points[i][d];
633 
635  internal::MappingQImplementation::
636  do_transform_real_to_unit_cell_internal<dim, spacedim>(
637  p_vec,
638  inverse_approximation.compute(p_vec),
639  support_points,
640  polynomials_1d,
641  renumber_lexicographic_to_hierarchic);
642 
643  // If the vectorized computation failed, it could be that only some of
644  // the lanes failed but others would have succeeded if we had let them
645  // compute alone without interference (like negative Jacobian
646  // determinants) from other SIMD lanes. Repeat the computation in this
647  // unlikely case with scalar arguments.
648  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
649  if (numbers::is_finite(unit_point[0][j]))
650  for (unsigned int d = 0; d < dim; ++d)
651  unit_points[i + j][d] = unit_point[d][j];
652  else
653  unit_points[i + j] = internal::MappingQImplementation::
654  do_transform_real_to_unit_cell_internal<dim, spacedim>(
655  real_points[i + j],
656  inverse_approximation.compute(real_points[i + j]),
657  support_points,
658  polynomials_1d,
659  renumber_lexicographic_to_hierarchic);
660  }
661  else
662  unit_points[i] = internal::MappingQImplementation::
663  do_transform_real_to_unit_cell_internal<dim, spacedim>(
664  real_points[i],
665  inverse_approximation.compute(real_points[i]),
666  support_points,
667  polynomials_1d,
668  renumber_lexicographic_to_hierarchic);
669 }
670 
671 
672 
673 template <int dim, int spacedim>
676 {
677  // add flags if the respective quantities are necessary to compute
678  // what we need. note that some flags appear in both the conditions
679  // and in subsequent set operations. this leads to some circular
680  // logic. the only way to treat this is to iterate. since there are
681  // 5 if-clauses in the loop, it will take at most 5 iterations to
682  // converge. do them:
683  UpdateFlags out = in;
684  for (unsigned int i = 0; i < 5; ++i)
685  {
686  // The following is a little incorrect:
687  // If not applied on a face,
688  // update_boundary_forms does not
689  // make sense. On the other hand,
690  // it is necessary on a
691  // face. Currently,
692  // update_boundary_forms is simply
693  // ignored for the interior of a
694  // cell.
696  out |= update_boundary_forms;
697 
698  if (out &
702 
703  if (out &
708 
709  // The contravariant transformation is used in the Piola
710  // transformation, which requires the determinant of the Jacobi
711  // matrix of the transformation. Because we have no way of
712  // knowing here whether the finite element wants to use the
713  // contravariant or the Piola transforms, we add the JxW values
714  // to the list of flags to be updated for each cell.
716  out |= update_volume_elements;
717 
718  // the same is true when computing normal vectors: they require
719  // the determinant of the Jacobian
720  if (out & update_normal_vectors)
721  out |= update_volume_elements;
722  }
723 
724  return out;
725 }
726 
727 
728 
729 template <int dim, int spacedim>
730 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
732  const Quadrature<dim> &q) const
733 {
734  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
735  std::make_unique<InternalData>(polynomial_degree);
736  auto &data = dynamic_cast<InternalData &>(*data_ptr);
737  data.initialize(this->requires_update_flags(update_flags), q, q.size());
738 
739  return data_ptr;
740 }
741 
742 
743 
744 template <int dim, int spacedim>
745 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
747  const UpdateFlags update_flags,
748  const hp::QCollection<dim - 1> &quadrature) const
749 {
750  AssertDimension(quadrature.size(), 1);
751 
752  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
753  std::make_unique<InternalData>(polynomial_degree);
754  auto &data = dynamic_cast<InternalData &>(*data_ptr);
755  data.initialize_face(this->requires_update_flags(update_flags),
757  ReferenceCells::get_hypercube<dim>(), quadrature[0]),
758  quadrature[0].size());
759 
760  return data_ptr;
761 }
762 
763 
764 
765 template <int dim, int spacedim>
766 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
768  const UpdateFlags update_flags,
769  const Quadrature<dim - 1> &quadrature) const
770 {
771  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
772  std::make_unique<InternalData>(polynomial_degree);
773  auto &data = dynamic_cast<InternalData &>(*data_ptr);
774  data.initialize_face(this->requires_update_flags(update_flags),
776  ReferenceCells::get_hypercube<dim>(), quadrature),
777  quadrature.size());
778 
779  return data_ptr;
780 }
781 
782 
783 
784 template <int dim, int spacedim>
787  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
788  const CellSimilarity::Similarity cell_similarity,
789  const Quadrature<dim> &quadrature,
790  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
792  &output_data) const
793 {
794  // ensure that the following static_cast is really correct:
795  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
796  ExcInternalError());
797  const InternalData &data = static_cast<const InternalData &>(internal_data);
798  data.output_data = &output_data;
799 
800  const unsigned int n_q_points = quadrature.size();
801 
802  // recompute the support points of the transformation of this
803  // cell. we tried to be clever here in an earlier version of the
804  // library by checking whether the cell is the same as the one we
805  // had visited last, but it turns out to be difficult to determine
806  // that because a cell for the purposes of a mapping is
807  // characterized not just by its (triangulation, level, index)
808  // triple, but also by the locations of its vertices, the manifold
809  // object attached to the cell and all of its bounding faces/edges,
810  // etc. to reliably test that the "cell" we are on is, therefore,
811  // not easily done
812  data.mapping_support_points = this->compute_mapping_support_points(cell);
813  data.cell_of_current_support_points = cell;
814 
815  // if the order of the mapping is greater than 1, then do not reuse any cell
816  // similarity information. This is necessary because the cell similarity
817  // value is computed with just cell vertices and does not take into account
818  // cell curvature.
819  const CellSimilarity::Similarity computed_cell_similarity =
820  (polynomial_degree == 1 && this->preserves_vertex_locations() ?
821  cell_similarity :
823 
824  if (dim > 1 && data.tensor_product_quadrature)
825  {
826  internal::MappingQImplementation::
827  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
828  computed_cell_similarity,
829  data,
830  output_data.quadrature_points,
831  output_data.jacobians,
832  output_data.inverse_jacobians,
833  output_data.jacobian_grads);
834  }
835  else
836  {
838  computed_cell_similarity,
839  data,
840  make_array_view(quadrature.get_points()),
841  polynomials_1d,
842  renumber_lexicographic_to_hierarchic,
843  output_data.quadrature_points,
844  output_data.jacobians,
845  output_data.inverse_jacobians);
846 
848  spacedim>(
849  computed_cell_similarity,
850  data,
851  make_array_view(quadrature.get_points()),
852  polynomials_1d,
853  renumber_lexicographic_to_hierarchic,
854  output_data.jacobian_grads);
855  }
856 
858  dim,
859  spacedim>(computed_cell_similarity,
860  data,
861  make_array_view(quadrature.get_points()),
862  polynomials_1d,
863  renumber_lexicographic_to_hierarchic,
864  output_data.jacobian_pushed_forward_grads);
865 
867  dim,
868  spacedim>(computed_cell_similarity,
869  data,
870  make_array_view(quadrature.get_points()),
871  polynomials_1d,
872  renumber_lexicographic_to_hierarchic,
873  output_data.jacobian_2nd_derivatives);
874 
875  internal::MappingQImplementation::
876  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
877  computed_cell_similarity,
878  data,
879  make_array_view(quadrature.get_points()),
880  polynomials_1d,
881  renumber_lexicographic_to_hierarchic,
883 
885  dim,
886  spacedim>(computed_cell_similarity,
887  data,
888  make_array_view(quadrature.get_points()),
889  polynomials_1d,
890  renumber_lexicographic_to_hierarchic,
891  output_data.jacobian_3rd_derivatives);
892 
893  internal::MappingQImplementation::
894  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
895  computed_cell_similarity,
896  data,
897  make_array_view(quadrature.get_points()),
898  polynomials_1d,
899  renumber_lexicographic_to_hierarchic,
901 
902  const UpdateFlags update_flags = data.update_each;
903  const std::vector<double> &weights = quadrature.get_weights();
904 
905  // Multiply quadrature weights by absolute value of Jacobian determinants or
906  // the area element g=sqrt(DX^t DX) in case of codim > 0
907 
908  if (update_flags & (update_normal_vectors | update_JxW_values))
909  {
910  Assert(!(update_flags & update_JxW_values) ||
911  (output_data.JxW_values.size() == n_q_points),
912  ExcDimensionMismatch(output_data.JxW_values.size(), n_q_points));
913 
914  Assert(!(update_flags & update_normal_vectors) ||
915  (output_data.normal_vectors.size() == n_q_points),
916  ExcDimensionMismatch(output_data.normal_vectors.size(),
917  n_q_points));
918 
919 
920  if (computed_cell_similarity != CellSimilarity::translation)
921  for (unsigned int point = 0; point < n_q_points; ++point)
922  {
923  if (dim == spacedim)
924  {
925  const double det = data.volume_elements[point];
926 
927  // check for distorted cells.
928 
929  // TODO: this allows for anisotropies of up to 1e6 in 3d and
930  // 1e12 in 2d. might want to find a finer
931  // (dimension-independent) criterion
932  Assert(det >
933  1e-12 * Utilities::fixed_power<dim>(
934  cell->diameter() / std::sqrt(double(dim))),
936  cell->center(), det, point)));
937 
938  output_data.JxW_values[point] = weights[point] * det;
939  }
940  // if dim==spacedim, then there is no cell normal to
941  // compute. since this is for FEValues (and not FEFaceValues),
942  // there are also no face normals to compute
943  else // codim>0 case
944  {
945  Tensor<1, spacedim> DX_t[dim];
946  for (unsigned int i = 0; i < spacedim; ++i)
947  for (unsigned int j = 0; j < dim; ++j)
948  DX_t[j][i] = output_data.jacobians[point][i][j];
949 
950  Tensor<2, dim> G; // First fundamental form
951  for (unsigned int i = 0; i < dim; ++i)
952  for (unsigned int j = 0; j < dim; ++j)
953  G[i][j] = DX_t[i] * DX_t[j];
954 
955  if (update_flags & update_JxW_values)
956  output_data.JxW_values[point] =
957  std::sqrt(determinant(G)) * weights[point];
958 
959  if (computed_cell_similarity ==
961  {
962  // we only need to flip the normal
963  if (update_flags & update_normal_vectors)
964  output_data.normal_vectors[point] *= -1.;
965  }
966  else
967  {
968  if (update_flags & update_normal_vectors)
969  {
970  Assert(spacedim == dim + 1,
971  ExcMessage(
972  "There is no (unique) cell normal for " +
974  "-dimensional cells in " +
975  Utilities::int_to_string(spacedim) +
976  "-dimensional space. This only works if the "
977  "space dimension is one greater than the "
978  "dimensionality of the mesh cells."));
979 
980  if (dim == 1)
981  output_data.normal_vectors[point] =
982  cross_product_2d(-DX_t[0]);
983  else // dim == 2
984  output_data.normal_vectors[point] =
985  cross_product_3d(DX_t[0], DX_t[1]);
986 
987  output_data.normal_vectors[point] /=
988  output_data.normal_vectors[point].norm();
989 
990  if (cell->direction_flag() == false)
991  output_data.normal_vectors[point] *= -1.;
992  }
993  }
994  } // codim>0 case
995  }
996  }
997 
998  return computed_cell_similarity;
999 }
1000 
1001 
1002 
1003 template <int dim, int spacedim>
1004 void
1006  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1007  const unsigned int face_no,
1008  const hp::QCollection<dim - 1> &quadrature,
1009  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1011  &output_data) const
1012 {
1013  AssertDimension(quadrature.size(), 1);
1014 
1015  // ensure that the following cast is really correct:
1016  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1017  ExcInternalError());
1018  const InternalData &data = static_cast<const InternalData &>(internal_data);
1019  data.output_data = &output_data;
1020 
1021  // if necessary, recompute the support points of the transformation of this
1022  // cell (note that we need to first check the triangulation pointer, since
1023  // otherwise the second test might trigger an exception if the
1024  // triangulations are not the same)
1025  if ((data.mapping_support_points.empty()) ||
1026  (&cell->get_triangulation() !=
1028  (cell != data.cell_of_current_support_points))
1029  {
1030  data.mapping_support_points = this->compute_mapping_support_points(cell);
1031  data.cell_of_current_support_points = cell;
1032  }
1033 
1035  *this,
1036  cell,
1037  face_no,
1040  ReferenceCells::get_hypercube<dim>(),
1041  face_no,
1042  cell->face_orientation(face_no),
1043  cell->face_flip(face_no),
1044  cell->face_rotation(face_no),
1045  quadrature[0].size()),
1046  quadrature[0],
1047  data,
1048  polynomials_1d,
1049  renumber_lexicographic_to_hierarchic,
1050  output_data);
1051 }
1052 
1053 
1054 
1055 template <int dim, int spacedim>
1056 void
1058  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1059  const unsigned int face_no,
1060  const unsigned int subface_no,
1061  const Quadrature<dim - 1> &quadrature,
1062  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1064  &output_data) const
1065 {
1066  // ensure that the following cast is really correct:
1067  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1068  ExcInternalError());
1069  const InternalData &data = static_cast<const InternalData &>(internal_data);
1070  data.output_data = &output_data;
1071 
1072  // if necessary, recompute the support points of the transformation of this
1073  // cell (note that we need to first check the triangulation pointer, since
1074  // otherwise the second test might trigger an exception if the
1075  // triangulations are not the same)
1076  if ((data.mapping_support_points.empty()) ||
1077  (&cell->get_triangulation() !=
1079  (cell != data.cell_of_current_support_points))
1080  {
1081  data.mapping_support_points = this->compute_mapping_support_points(cell);
1082  data.cell_of_current_support_points = cell;
1083  }
1084 
1086  *this,
1087  cell,
1088  face_no,
1089  subface_no,
1091  ReferenceCells::get_hypercube<dim>(),
1092  face_no,
1093  subface_no,
1094  cell->face_orientation(face_no),
1095  cell->face_flip(face_no),
1096  cell->face_rotation(face_no),
1097  quadrature.size(),
1098  cell->subface_case(face_no)),
1099  quadrature,
1100  data,
1101  polynomials_1d,
1102  renumber_lexicographic_to_hierarchic,
1103  output_data);
1104 }
1105 
1106 
1107 
1108 template <int dim, int spacedim>
1109 void
1111  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1113  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1115  &output_data) const
1116 {
1117  Assert(dim == spacedim, ExcNotImplemented());
1118 
1119  // ensure that the following static_cast is really correct:
1120  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1121  ExcInternalError());
1122  const InternalData &data = static_cast<const InternalData &>(internal_data);
1123  data.output_data = &output_data;
1124 
1125  const unsigned int n_q_points = quadrature.size();
1126 
1127  data.mapping_support_points = this->compute_mapping_support_points(cell);
1128  data.cell_of_current_support_points = cell;
1129 
1132  data,
1133  make_array_view(quadrature.get_points()),
1134  polynomials_1d,
1135  renumber_lexicographic_to_hierarchic,
1136  output_data.quadrature_points,
1137  output_data.jacobians,
1138  output_data.inverse_jacobians);
1139 
1140  internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1142  data,
1143  make_array_view(quadrature.get_points()),
1144  polynomials_1d,
1145  renumber_lexicographic_to_hierarchic,
1146  output_data.jacobian_grads);
1147 
1149  dim,
1150  spacedim>(CellSimilarity::none,
1151  data,
1152  make_array_view(quadrature.get_points()),
1153  polynomials_1d,
1154  renumber_lexicographic_to_hierarchic,
1155  output_data.jacobian_pushed_forward_grads);
1156 
1158  dim,
1159  spacedim>(CellSimilarity::none,
1160  data,
1161  make_array_view(quadrature.get_points()),
1162  polynomials_1d,
1163  renumber_lexicographic_to_hierarchic,
1164  output_data.jacobian_2nd_derivatives);
1165 
1166  internal::MappingQImplementation::
1167  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1169  data,
1170  make_array_view(quadrature.get_points()),
1171  polynomials_1d,
1172  renumber_lexicographic_to_hierarchic,
1174 
1176  dim,
1177  spacedim>(CellSimilarity::none,
1178  data,
1179  make_array_view(quadrature.get_points()),
1180  polynomials_1d,
1181  renumber_lexicographic_to_hierarchic,
1182  output_data.jacobian_3rd_derivatives);
1183 
1184  internal::MappingQImplementation::
1185  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1187  data,
1188  make_array_view(quadrature.get_points()),
1189  polynomials_1d,
1190  renumber_lexicographic_to_hierarchic,
1192 
1193  const UpdateFlags update_flags = data.update_each;
1194  const std::vector<double> &weights = quadrature.get_weights();
1195 
1196  if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1197  {
1198  AssertDimension(output_data.JxW_values.size(), n_q_points);
1199 
1200  Assert(!(update_flags & update_normal_vectors) ||
1201  (output_data.normal_vectors.size() == n_q_points),
1202  ExcDimensionMismatch(output_data.normal_vectors.size(),
1203  n_q_points));
1204 
1205 
1206  for (unsigned int point = 0; point < n_q_points; ++point)
1207  {
1208  const double det = data.volume_elements[point];
1209 
1210  // check for distorted cells.
1211 
1212  // TODO: this allows for anisotropies of up to 1e6 in 3d and
1213  // 1e12 in 2d. might want to find a finer
1214  // (dimension-independent) criterion
1215  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1216  cell->diameter() / std::sqrt(double(dim))),
1218  cell->center(), det, point)));
1219 
1220  // The normals are n = J^{-T} * \hat{n} before normalizing.
1221  Tensor<1, spacedim> normal;
1222  for (unsigned int d = 0; d < spacedim; d++)
1223  normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
1224  quadrature.normal_vector(point);
1225 
1226  output_data.JxW_values[point] = weights[point] * det * normal.norm();
1227 
1228  if ((update_flags & update_normal_vectors) != 0u)
1229  {
1230  normal /= normal.norm();
1231  output_data.normal_vectors[point] = normal;
1232  }
1233  }
1234  }
1235 }
1236 
1237 
1238 
1239 template <int dim, int spacedim>
1240 void
1242  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1243  const ArrayView<const Point<dim>> &unit_points,
1244  const UpdateFlags update_flags,
1246  &output_data) const
1247 {
1248  if (update_flags == update_default)
1249  return;
1250 
1251  Assert(update_flags & update_inverse_jacobians ||
1252  update_flags & update_jacobians ||
1253  update_flags & update_quadrature_points,
1254  ExcNotImplemented());
1255 
1256  output_data.initialize(unit_points.size(), update_flags);
1257 
1258  auto internal_data =
1259  this->get_data(update_flags,
1260  Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
1261  unit_points.end())));
1262  const InternalData &data = static_cast<const InternalData &>(*internal_data);
1263  data.output_data = &output_data;
1264  data.mapping_support_points = this->compute_mapping_support_points(cell);
1265 
1268  data,
1269  unit_points,
1270  polynomials_1d,
1271  renumber_lexicographic_to_hierarchic,
1272  output_data.quadrature_points,
1273  output_data.jacobians,
1274  output_data.inverse_jacobians);
1275 }
1276 
1277 
1278 
1279 template <int dim, int spacedim>
1280 void
1282  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1283  const unsigned int face_no,
1284  const Quadrature<dim - 1> &face_quadrature,
1285  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1287  &output_data) const
1288 {
1289  if (face_quadrature.get_points().empty())
1290  return;
1291 
1292  // ensure that the following static_cast is really correct:
1293  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1294  ExcInternalError());
1295  const InternalData &data = static_cast<const InternalData &>(internal_data);
1296 
1297  data.mapping_support_points = this->compute_mapping_support_points(cell);
1298  data.output_data = &output_data;
1299 
1301  *this,
1302  cell,
1303  face_no,
1306  face_quadrature,
1307  data,
1308  polynomials_1d,
1309  renumber_lexicographic_to_hierarchic,
1310  output_data);
1311 }
1312 
1313 
1314 
1315 template <int dim, int spacedim>
1316 void
1318  const ArrayView<const Tensor<1, dim>> &input,
1319  const MappingKind mapping_kind,
1320  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1321  const ArrayView<Tensor<1, spacedim>> &output) const
1322 {
1324  mapping_kind,
1325  mapping_data,
1326  output);
1327 }
1328 
1329 
1330 
1331 template <int dim, int spacedim>
1332 void
1334  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1335  const MappingKind mapping_kind,
1336  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1337  const ArrayView<Tensor<2, spacedim>> &output) const
1338 {
1340  mapping_kind,
1341  mapping_data,
1342  output);
1343 }
1344 
1345 
1346 
1347 template <int dim, int spacedim>
1348 void
1350  const ArrayView<const Tensor<2, dim>> &input,
1351  const MappingKind mapping_kind,
1352  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1353  const ArrayView<Tensor<2, spacedim>> &output) const
1354 {
1355  switch (mapping_kind)
1356  {
1357  case mapping_contravariant:
1359  mapping_kind,
1360  mapping_data,
1361  output);
1362  return;
1363 
1368  mapping_kind,
1369  mapping_data,
1370  output);
1371  return;
1372  default:
1373  Assert(false, ExcNotImplemented());
1374  }
1375 }
1376 
1377 
1378 
1379 template <int dim, int spacedim>
1380 void
1382  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1383  const MappingKind mapping_kind,
1384  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1385  const ArrayView<Tensor<3, spacedim>> &output) const
1386 {
1387  AssertDimension(input.size(), output.size());
1388  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1389  ExcInternalError());
1391  &data = *static_cast<const InternalData &>(mapping_data).output_data;
1392 
1393  switch (mapping_kind)
1394  {
1396  {
1397  Assert(!data.inverse_jacobians.empty(),
1399  "update_covariant_transformation"));
1400 
1401  for (unsigned int q = 0; q < output.size(); ++q)
1402  for (unsigned int i = 0; i < spacedim; ++i)
1403  for (unsigned int j = 0; j < spacedim; ++j)
1404  {
1405  double tmp[dim];
1406  const DerivativeForm<1, dim, spacedim> covariant =
1407  data.inverse_jacobians[q].transpose();
1408  for (unsigned int K = 0; K < dim; ++K)
1409  {
1410  tmp[K] = covariant[j][0] * input[q][i][0][K];
1411  for (unsigned int J = 1; J < dim; ++J)
1412  tmp[K] += covariant[j][J] * input[q][i][J][K];
1413  }
1414  for (unsigned int k = 0; k < spacedim; ++k)
1415  {
1416  output[q][i][j][k] = covariant[k][0] * tmp[0];
1417  for (unsigned int K = 1; K < dim; ++K)
1418  output[q][i][j][k] += covariant[k][K] * tmp[K];
1419  }
1420  }
1421  return;
1422  }
1423 
1424  default:
1425  Assert(false, ExcNotImplemented());
1426  }
1427 }
1428 
1429 
1430 
1431 template <int dim, int spacedim>
1432 void
1434  const ArrayView<const Tensor<3, dim>> &input,
1435  const MappingKind mapping_kind,
1436  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1437  const ArrayView<Tensor<3, spacedim>> &output) const
1438 {
1439  switch (mapping_kind)
1440  {
1441  case mapping_piola_hessian:
1445  mapping_kind,
1446  mapping_data,
1447  output);
1448  return;
1449  default:
1450  Assert(false, ExcNotImplemented());
1451  }
1452 }
1453 
1454 
1455 
1456 template <int dim, int spacedim>
1457 void
1459  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1460  std::vector<Point<spacedim>> &a) const
1461 {
1462  // if we only need the midpoint, then ask for it.
1463  if (this->polynomial_degree == 2)
1464  {
1465  for (unsigned int line_no = 0;
1466  line_no < GeometryInfo<dim>::lines_per_cell;
1467  ++line_no)
1468  {
1469  const typename Triangulation<dim, spacedim>::line_iterator line =
1470  (dim == 1 ?
1471  static_cast<
1473  cell->line(line_no));
1474 
1475  const Manifold<dim, spacedim> &manifold =
1476  ((line->manifold_id() == numbers::flat_manifold_id) &&
1477  (dim < spacedim) ?
1478  cell->get_manifold() :
1479  line->get_manifold());
1480  a.push_back(manifold.get_new_point_on_line(line));
1481  }
1482  }
1483  else
1484  // otherwise call the more complicated functions and ask for inner points
1485  // from the manifold description
1486  {
1487  std::vector<Point<spacedim>> tmp_points;
1488  for (unsigned int line_no = 0;
1489  line_no < GeometryInfo<dim>::lines_per_cell;
1490  ++line_no)
1491  {
1492  const typename Triangulation<dim, spacedim>::line_iterator line =
1493  (dim == 1 ?
1494  static_cast<
1496  cell->line(line_no));
1497 
1498  const Manifold<dim, spacedim> &manifold =
1499  ((line->manifold_id() == numbers::flat_manifold_id) &&
1500  (dim < spacedim) ?
1501  cell->get_manifold() :
1502  line->get_manifold());
1503 
1504  const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1505  const std::array<Point<spacedim>, 2> vertices{
1506  {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
1507  cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
1508 
1509  const std::size_t n_rows =
1510  support_point_weights_perimeter_to_interior[0].size(0);
1511  a.resize(a.size() + n_rows);
1512  auto a_view = make_array_view(a.end() - n_rows, a.end());
1513  manifold.get_new_points(
1514  make_array_view(vertices.begin(), vertices.end()),
1515  support_point_weights_perimeter_to_interior[0],
1516  a_view);
1517  }
1518  }
1519 }
1520 
1521 
1522 
1523 template <>
1524 void
1527  std::vector<Point<3>> &a) const
1528 {
1529  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1530 
1531  // used if face quad at boundary or entirely in the interior of the domain
1532  std::vector<Point<3>> tmp_points;
1533 
1534  // loop over all faces and collect points on them
1535  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1536  {
1537  const Triangulation<3>::face_iterator face = cell->face(face_no);
1538 
1539 #ifdef DEBUG
1540  const bool face_orientation = cell->face_orientation(face_no),
1541  face_flip = cell->face_flip(face_no),
1542  face_rotation = cell->face_rotation(face_no);
1543  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1544  lines_per_face = GeometryInfo<3>::lines_per_face;
1545 
1546  // some sanity checks up front
1547  for (unsigned int i = 0; i < vertices_per_face; ++i)
1548  Assert(face->vertex_index(i) ==
1549  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1550  face_no, i, face_orientation, face_flip, face_rotation)),
1551  ExcInternalError());
1552 
1553  // indices of the lines that bound a face are given by GeometryInfo<3>::
1554  // face_to_cell_lines
1555  for (unsigned int i = 0; i < lines_per_face; ++i)
1556  Assert(face->line(i) ==
1558  face_no, i, face_orientation, face_flip, face_rotation)),
1559  ExcInternalError());
1560 #endif
1561  // extract the points surrounding a quad from the points
1562  // already computed. First get the 4 vertices and then the points on
1563  // the four lines
1564  boost::container::small_vector<Point<3>, 200> tmp_points(
1566  GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1567  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1568  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1569  if (polynomial_degree > 1)
1570  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1571  ++line)
1572  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1573  tmp_points[4 + line * (polynomial_degree - 1) + i] =
1575  (polynomial_degree - 1) *
1576  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
1577  i];
1578 
1579  const std::size_t n_rows =
1580  support_point_weights_perimeter_to_interior[1].size(0);
1581  a.resize(a.size() + n_rows);
1582  auto a_view = make_array_view(a.end() - n_rows, a.end());
1583  face->get_manifold().get_new_points(
1584  make_array_view(tmp_points.begin(), tmp_points.end()),
1585  support_point_weights_perimeter_to_interior[1],
1586  a_view);
1587  }
1588 }
1589 
1590 
1591 
1592 template <>
1593 void
1596  std::vector<Point<3>> &a) const
1597 {
1598  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1599  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1600  vertices[i] = cell->vertex(i);
1601 
1602  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1604  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1605  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1606  {
1607  Point<2> point(line_support_points[q1 + 1][0],
1608  line_support_points[q2 + 1][0]);
1609  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1610  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1611  }
1612 
1613  const std::size_t n_rows = weights.size(0);
1614  a.resize(a.size() + n_rows);
1615  auto a_view = make_array_view(a.end() - n_rows, a.end());
1616  cell->get_manifold().get_new_points(
1617  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1618 }
1619 
1620 
1621 
1622 template <int dim, int spacedim>
1623 void
1626  std::vector<Point<spacedim>> &) const
1627 {
1628  Assert(false, ExcInternalError());
1629 }
1630 
1631 
1632 
1633 template <int dim, int spacedim>
1634 std::vector<Point<spacedim>>
1636  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1637 {
1638  // get the vertices first
1639  std::vector<Point<spacedim>> a;
1640  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1641  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1642  a.push_back(cell->vertex(i));
1643 
1644  if (this->polynomial_degree > 1)
1645  {
1646  // check if all entities have the same manifold id which is when we can
1647  // simply ask the manifold for all points. the transfinite manifold can
1648  // do the interpolation better than this class, so if we detect that we
1649  // do not have to change anything here
1650  Assert(dim <= 3, ExcImpossibleInDim(dim));
1651  bool all_manifold_ids_are_equal = (dim == spacedim);
1652  if (all_manifold_ids_are_equal &&
1654  &cell->get_manifold()) == nullptr)
1655  {
1656  for (auto f : GeometryInfo<dim>::face_indices())
1657  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1658  all_manifold_ids_are_equal = false;
1659 
1660  if (dim == 3)
1661  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1662  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1663  all_manifold_ids_are_equal = false;
1664  }
1665 
1666  if (all_manifold_ids_are_equal)
1667  {
1668  const std::size_t n_rows = support_point_weights_cell.size(0);
1669  a.resize(a.size() + n_rows);
1670  auto a_view = make_array_view(a.end() - n_rows, a.end());
1671  cell->get_manifold().get_new_points(make_array_view(a.begin(),
1672  a.end() - n_rows),
1673  support_point_weights_cell,
1674  a_view);
1675  }
1676  else
1677  switch (dim)
1678  {
1679  case 1:
1680  add_line_support_points(cell, a);
1681  break;
1682  case 2:
1683  // in 2d, add the points on the four bounding lines to the
1684  // exterior (outer) points
1685  add_line_support_points(cell, a);
1686 
1687  // then get the interior support points
1688  if (dim != spacedim)
1689  add_quad_support_points(cell, a);
1690  else
1691  {
1692  const std::size_t n_rows =
1693  support_point_weights_perimeter_to_interior[1].size(0);
1694  a.resize(a.size() + n_rows);
1695  auto a_view = make_array_view(a.end() - n_rows, a.end());
1696  cell->get_manifold().get_new_points(
1697  make_array_view(a.begin(), a.end() - n_rows),
1698  support_point_weights_perimeter_to_interior[1],
1699  a_view);
1700  }
1701  break;
1702 
1703  case 3:
1704  // in 3d also add the points located on the boundary faces
1705  add_line_support_points(cell, a);
1706  add_quad_support_points(cell, a);
1707 
1708  // then compute the interior points
1709  {
1710  const std::size_t n_rows =
1711  support_point_weights_perimeter_to_interior[2].size(0);
1712  a.resize(a.size() + n_rows);
1713  auto a_view = make_array_view(a.end() - n_rows, a.end());
1714  cell->get_manifold().get_new_points(
1715  make_array_view(a.begin(), a.end() - n_rows),
1716  support_point_weights_perimeter_to_interior[2],
1717  a_view);
1718  }
1719  break;
1720 
1721  default:
1722  Assert(false, ExcNotImplemented());
1723  break;
1724  }
1725  }
1726 
1727  return a;
1728 }
1729 
1730 
1731 
1732 template <int dim, int spacedim>
1735  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1736 {
1737  return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1738 }
1739 
1740 
1741 
1742 template <int dim, int spacedim>
1743 bool
1745  const ReferenceCell &reference_cell) const
1746 {
1747  Assert(dim == reference_cell.get_dimension(),
1748  ExcMessage("The dimension of your mapping (" +
1749  Utilities::to_string(dim) +
1750  ") and the reference cell cell_type (" +
1751  Utilities::to_string(reference_cell.get_dimension()) +
1752  " ) do not agree."));
1753 
1754  return reference_cell.is_hyper_cube();
1755 }
1756 
1757 
1758 
1759 //--------------------------- Explicit instantiations -----------------------
1760 #include "mapping_q.inst"
1761 
1762 
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:838
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
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition: mapping_q.h:445
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:64
std::vector< Point< spacedim > > mapping_support_points
Definition: mapping_q.h:439
internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > * output_data
Definition: mapping_q.h:459
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:451
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:1458
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition: mapping_q.h:557
const Table< 2, double > support_point_weights_cell
Definition: mapping_q.h:605
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_q.cc:1635
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:731
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:1241
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:1281
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_q.cc:1734
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:767
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:1057
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:786
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:1110
const unsigned int polynomial_degree
Definition: mapping_q.h:533
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:1317
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:746
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:324
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition: mapping_q.h:591
const std::vector< Point< 1 > > line_support_points
Definition: mapping_q.h:543
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:290
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition: mapping_q.h:550
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:272
const std::vector< Point< dim > > unit_cell_support_points
Definition: mapping_q.h:569
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:593
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:473
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:675
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:1005
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition: mapping_q.cc:1744
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
Definition: mapping_q.cc:1624
unsigned int get_degree() const
Definition: mapping_q.cc:281
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:265
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
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:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
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_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
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:975
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:539
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)