Reference documentation for deal.II version GIT f3cf314fd8 2023-02-08 12:00: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 - 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 
38 #include <boost/container/small_vector.hpp>
40 
41 #include <algorithm>
42 #include <array>
43 #include <cmath>
44 #include <limits>
45 #include <memory>
46 #include <numeric>
47 
48 
50 
51 
52 template <int dim, int spacedim>
54  const unsigned int polynomial_degree)
55  : polynomial_degree(polynomial_degree)
56  , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
57  , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
58  , tensor_product_quadrature(false)
59 {}
60 
61 
62 
63 template <int dim, int spacedim>
64 std::size_t
66 {
67  return (
70  MemoryConsumption::memory_consumption(shape_derivatives) +
73  MemoryConsumption::memory_consumption(unit_tangentials) +
75  MemoryConsumption::memory_consumption(mapping_support_points) +
76  MemoryConsumption::memory_consumption(cell_of_current_support_points) +
77  MemoryConsumption::memory_consumption(volume_elements) +
79  MemoryConsumption::memory_consumption(n_shape_functions));
80 }
81 
82 
83 
84 template <int dim, int spacedim>
85 void
87  const UpdateFlags update_flags,
88  const Quadrature<dim> &q,
89  const unsigned int n_original_q_points)
90 {
91  // store the flags in the internal data object so we can access them
92  // in fill_fe_*_values()
93  this->update_each = update_flags;
94 
95  const unsigned int n_q_points = q.size();
96 
97  const bool needs_higher_order_terms =
98  this->update_each &
103 
104  if (this->update_each & update_covariant_transformation)
105  covariant.resize(n_original_q_points);
106 
107  if (this->update_each & update_contravariant_transformation)
108  contravariant.resize(n_original_q_points);
109 
110  if (this->update_each & update_volume_elements)
111  volume_elements.resize(n_original_q_points);
112 
113  tensor_product_quadrature = q.is_tensor_product();
114 
115  // use of MatrixFree only for higher order elements and with more than one
116  // point where tensor products do not make sense
117  if (polynomial_degree < 2 || n_q_points == 1)
118  tensor_product_quadrature = false;
119 
120  if (dim > 1)
121  {
122  // find out if the one-dimensional formula is the same
123  // in all directions
124  if (tensor_product_quadrature)
125  {
126  const std::array<Quadrature<1>, dim> &quad_array =
128  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
129  {
130  if (quad_array[i - 1].size() != quad_array[i].size())
131  {
132  tensor_product_quadrature = false;
133  break;
134  }
135  else
136  {
137  const std::vector<Point<1>> &points_1 =
138  quad_array[i - 1].get_points();
139  const std::vector<Point<1>> &points_2 =
140  quad_array[i].get_points();
141  const std::vector<double> &weights_1 =
142  quad_array[i - 1].get_weights();
143  const std::vector<double> &weights_2 =
144  quad_array[i].get_weights();
145  for (unsigned int j = 0; j < quad_array[i].size(); ++j)
146  {
147  if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
148  std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
149  {
150  tensor_product_quadrature = false;
151  break;
152  }
153  }
154  }
155  }
156 
157  if (tensor_product_quadrature)
158  {
159  // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
160  // numbering manually (building an FE_Q<dim> is relatively
161  // expensive due to constraints)
162  const FE_DGQ<1> fe(polynomial_degree);
163  shape_info.reinit(q.get_tensor_basis()[0], fe);
164  shape_info.lexicographic_numbering =
165  FETools::lexicographic_to_hierarchic_numbering<dim>(
167  shape_info.n_q_points = q.size();
168  shape_info.dofs_per_component_on_cell =
170  }
171  }
172  }
173 
174  // Only fill the big arrays on demand in case we cannot use the tensor
175  // product quadrature code path
176  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
177  {
178  // see if we need the (transformation) shape function values
179  // and/or gradients and resize the necessary arrays
180  if (this->update_each & update_quadrature_points)
181  shape_values.resize(n_shape_functions * n_q_points);
182 
183  if (this->update_each &
193  shape_derivatives.resize(n_shape_functions * n_q_points);
194 
195  if (this->update_each &
197  shape_second_derivatives.resize(n_shape_functions * n_q_points);
198 
199  if (this->update_each & (update_jacobian_2nd_derivatives |
201  shape_third_derivatives.resize(n_shape_functions * n_q_points);
202 
203  if (this->update_each & (update_jacobian_3rd_derivatives |
205  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
206 
207  // now also fill the various fields with their correct values
208  compute_shape_function_values(q.get_points());
209  }
210 }
211 
212 
213 
214 template <int dim, int spacedim>
215 void
217  const UpdateFlags update_flags,
218  const Quadrature<dim> &q,
219  const unsigned int n_original_q_points)
220 {
221  initialize(update_flags, q, n_original_q_points);
222 
223  if (dim > 1 && tensor_product_quadrature)
224  {
225  constexpr unsigned int facedim = dim - 1;
226  const FE_DGQ<1> fe(polynomial_degree);
227  shape_info.reinit(q.get_tensor_basis()[0], fe);
228  shape_info.lexicographic_numbering =
229  FETools::lexicographic_to_hierarchic_numbering<facedim>(
231  shape_info.n_q_points = n_original_q_points;
232  shape_info.dofs_per_component_on_cell =
234  }
235 
236  if (dim > 1)
237  {
238  if (this->update_each &
241  {
242  aux.resize(dim - 1,
243  AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
244 
245  // Compute tangentials to the unit cell.
246  for (const unsigned int i : GeometryInfo<dim>::face_indices())
247  {
248  unit_tangentials[i].resize(n_original_q_points);
249  std::fill(unit_tangentials[i].begin(),
250  unit_tangentials[i].end(),
252  if (dim > 2)
253  {
254  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
255  .resize(n_original_q_points);
256  std::fill(
257  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
258  .begin(),
259  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
260  .end(),
262  }
263  }
264  }
265  }
266 }
267 
268 
269 
270 template <int dim, int spacedim>
271 void
273  const std::vector<Point<dim>> &unit_points)
274 {
275  const unsigned int n_points = unit_points.size();
276 
277  // Construct the tensor product polynomials used as shape functions for
278  // the Qp mapping of cells at the boundary.
279  const TensorProductPolynomials<dim> tensor_pols(
281  line_support_points.get_points()));
282  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
283 
284  // then also construct the mapping from lexicographic to the Qp shape
285  // function numbering
286  const std::vector<unsigned int> renumber =
287  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
288 
289  std::vector<double> values;
290  std::vector<Tensor<1, dim>> grads;
291  if (shape_values.size() != 0)
292  {
293  Assert(shape_values.size() == n_shape_functions * n_points,
294  ExcInternalError());
295  values.resize(n_shape_functions);
296  }
297  if (shape_derivatives.size() != 0)
298  {
299  Assert(shape_derivatives.size() == n_shape_functions * n_points,
300  ExcInternalError());
301  grads.resize(n_shape_functions);
302  }
303 
304  std::vector<Tensor<2, dim>> grad2;
305  if (shape_second_derivatives.size() != 0)
306  {
307  Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
308  ExcInternalError());
309  grad2.resize(n_shape_functions);
310  }
311 
312  std::vector<Tensor<3, dim>> grad3;
313  if (shape_third_derivatives.size() != 0)
314  {
315  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
316  ExcInternalError());
317  grad3.resize(n_shape_functions);
318  }
319 
320  std::vector<Tensor<4, dim>> grad4;
321  if (shape_fourth_derivatives.size() != 0)
322  {
323  Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
324  ExcInternalError());
325  grad4.resize(n_shape_functions);
326  }
327 
328 
329  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
330  shape_second_derivatives.size() != 0 ||
331  shape_third_derivatives.size() != 0 ||
332  shape_fourth_derivatives.size() != 0)
333  for (unsigned int point = 0; point < n_points; ++point)
334  {
335  tensor_pols.evaluate(
336  unit_points[point], values, grads, grad2, grad3, grad4);
337 
338  if (shape_values.size() != 0)
339  for (unsigned int i = 0; i < n_shape_functions; ++i)
340  shape(point, i) = values[renumber[i]];
341 
342  if (shape_derivatives.size() != 0)
343  for (unsigned int i = 0; i < n_shape_functions; ++i)
344  derivative(point, i) = grads[renumber[i]];
345 
346  if (shape_second_derivatives.size() != 0)
347  for (unsigned int i = 0; i < n_shape_functions; ++i)
348  second_derivative(point, i) = grad2[renumber[i]];
349 
350  if (shape_third_derivatives.size() != 0)
351  for (unsigned int i = 0; i < n_shape_functions; ++i)
352  third_derivative(point, i) = grad3[renumber[i]];
353 
354  if (shape_fourth_derivatives.size() != 0)
355  for (unsigned int i = 0; i < n_shape_functions; ++i)
356  fourth_derivative(point, i) = grad4[renumber[i]];
357  }
358 }
359 
360 
361 
362 template <int dim, int spacedim>
364  : polynomial_degree(p)
366  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
367  , polynomials_1d(
372  internal::MappingQImplementation::unit_support_points<dim>(
376  internal::MappingQImplementation::
378  this->polynomial_degree,
379  dim))
381  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
382  this->polynomial_degree))
383 {
384  Assert(p >= 1,
385  ExcMessage("It only makes sense to create polynomial mappings "
386  "with a polynomial degree greater or equal to one."));
387 }
388 
389 
390 
391 template <int dim, int spacedim>
392 MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
393  : polynomial_degree(p)
394  , line_support_points(
395  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
396  , polynomials_1d(
397  Polynomials::generate_complete_Lagrange_basis(line_support_points))
398  , renumber_lexicographic_to_hierarchic(
400  , unit_cell_support_points(
401  internal::MappingQImplementation::unit_support_points<dim>(
402  line_support_points,
403  renumber_lexicographic_to_hierarchic))
404  , support_point_weights_perimeter_to_interior(
405  internal::MappingQImplementation::
407  this->polynomial_degree,
408  dim))
409  , support_point_weights_cell(
410  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
411  this->polynomial_degree))
412 {
413  Assert(p >= 1,
414  ExcMessage("It only makes sense to create polynomial mappings "
415  "with a polynomial degree greater or equal to one."));
416 }
417 
418 
419 
420 template <int dim, int spacedim>
422  : polynomial_degree(mapping.polynomial_degree)
423  , line_support_points(mapping.line_support_points)
424  , polynomials_1d(mapping.polynomials_1d)
425  , renumber_lexicographic_to_hierarchic(
426  mapping.renumber_lexicographic_to_hierarchic)
427  , support_point_weights_perimeter_to_interior(
428  mapping.support_point_weights_perimeter_to_interior)
429  , support_point_weights_cell(mapping.support_point_weights_cell)
430 {}
431 
432 
433 
434 template <int dim, int spacedim>
435 std::unique_ptr<Mapping<dim, spacedim>>
437 {
438  return std::make_unique<MappingQ<dim, spacedim>>(*this);
439 }
440 
441 
442 
443 template <int dim, int spacedim>
444 unsigned int
446 {
447  return polynomial_degree;
448 }
449 
450 
451 
452 template <int dim, int spacedim>
455  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
456  const Point<dim> & p) const
457 {
459  polynomials_1d,
460  this->compute_mapping_support_points(cell),
461  p,
462  polynomials_1d.size() == 2,
463  renumber_lexicographic_to_hierarchic)
464  .first);
465 }
466 
467 
468 // In the code below, GCC tries to instantiate MappingQ<3,4> when
469 // seeing which of the overloaded versions of
470 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
471 // error messages and, generally, nothing very good. Avoid this by ensuring
472 // that this class exists, but does not have an inner InternalData
473 // type, thereby ruling out the codim-1 version of the function
474 // below when doing overload resolution.
475 template <>
476 class MappingQ<3, 4>
477 {};
478 
479 
480 
481 // visual studio freaks out when trying to determine if
482 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
483 // candidate. So instead of letting the compiler pick the correct overload, we
484 // use template specialization to make sure we pick up the right function to
485 // call:
486 
487 template <int dim, int spacedim>
491  const Point<spacedim> &,
492  const Point<dim> &) const
493 {
494  // default implementation (should never be called)
495  Assert(false, ExcInternalError());
496  return {};
497 }
498 
499 
500 
501 template <>
502 Point<1>
505  const Point<1> & p,
506  const Point<1> & initial_p_unit) const
507 {
508  // dispatch to the various specializations for spacedim=dim,
509  // spacedim=dim+1, etc
510  return internal::MappingQImplementation::
511  do_transform_real_to_unit_cell_internal<1>(
512  p,
513  initial_p_unit,
514  this->compute_mapping_support_points(cell),
515  polynomials_1d,
516  renumber_lexicographic_to_hierarchic);
517 }
518 
519 
520 
521 template <>
522 Point<2>
525  const Point<2> & p,
526  const Point<2> & initial_p_unit) const
527 {
528  return internal::MappingQImplementation::
529  do_transform_real_to_unit_cell_internal<2>(
530  p,
531  initial_p_unit,
532  this->compute_mapping_support_points(cell),
533  polynomials_1d,
534  renumber_lexicographic_to_hierarchic);
535 }
536 
537 
538 
539 template <>
540 Point<3>
543  const Point<3> & p,
544  const Point<3> & initial_p_unit) const
545 {
546  return internal::MappingQImplementation::
547  do_transform_real_to_unit_cell_internal<3>(
548  p,
549  initial_p_unit,
550  this->compute_mapping_support_points(cell),
551  polynomials_1d,
552  renumber_lexicographic_to_hierarchic);
553 }
554 
555 
556 
557 template <>
558 Point<1>
561  const Point<2> & p,
562  const Point<1> & initial_p_unit) const
563 {
564  const int dim = 1;
565  const int spacedim = 2;
566 
567  const Quadrature<dim> point_quadrature(initial_p_unit);
568 
570  if (spacedim > dim)
571  update_flags |= update_jacobian_grads;
572  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
573  get_data(update_flags, point_quadrature));
574 
575  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
576 
577  // dispatch to the various specializations for spacedim=dim,
578  // spacedim=dim+1, etc
579  return internal::MappingQImplementation::
580  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
581  p,
582  initial_p_unit,
583  *mdata);
584 }
585 
586 
587 
588 template <>
589 Point<2>
592  const Point<3> & p,
593  const Point<2> & initial_p_unit) const
594 {
595  const int dim = 2;
596  const int spacedim = 3;
597 
598  const Quadrature<dim> point_quadrature(initial_p_unit);
599 
601  if (spacedim > dim)
602  update_flags |= update_jacobian_grads;
603  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
604  get_data(update_flags, point_quadrature));
605 
606  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
607 
608  // dispatch to the various specializations for spacedim=dim,
609  // spacedim=dim+1, etc
610  return internal::MappingQImplementation::
611  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
612  p,
613  initial_p_unit,
614  *mdata);
615 }
616 
617 template <>
618 Point<1>
621  const Point<3> &,
622  const Point<1> &) const
623 {
624  Assert(false, ExcNotImplemented());
625  return {};
626 }
627 
628 
629 
630 template <int dim, int spacedim>
633  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
634  const Point<spacedim> & p) const
635 {
636  // Use an exact formula if one is available. this is only the case
637  // for Q1 mappings in 1d, and in 2d if dim==spacedim
638  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
639  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
640  {
641  // The dimension-dependent algorithms are much faster (about 25-45x in
642  // 2D) but fail most of the time when the given point (p) is not in the
643  // cell. The dimension-independent Newton algorithm given below is
644  // slower, but more robust (though it still sometimes fails). Therefore
645  // this function implements the following strategy based on the
646  // p's dimension:
647  //
648  // * In 1D this mapping is linear, so the mapping is always invertible
649  // (and the exact formula is known) as long as the cell has non-zero
650  // length.
651  // * In 2D the exact (quadratic) formula is called first. If either the
652  // exact formula does not succeed (negative discriminant in the
653  // quadratic formula) or succeeds but finds a solution outside of the
654  // unit cell, then the Newton solver is called. The rationale for the
655  // second choice is that the exact formula may provide two different
656  // answers when mapping a point outside of the real cell, but the
657  // Newton solver (if it converges) will only return one answer.
658  // Otherwise the exact formula successfully found a point in the unit
659  // cell and that value is returned.
660  // * In 3D there is no (known to the authors) exact formula, so the Newton
661  // algorithm is used.
662  const auto vertices_ = this->get_vertices(cell);
663 
664  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
665  vertices;
666  for (unsigned int i = 0; i < vertices.size(); ++i)
667  vertices[i] = vertices_[i];
668 
669  try
670  {
671  switch (dim)
672  {
673  case 1:
674  {
675  // formula not subject to any issues in 1d
676  if (spacedim == 1)
678  vertices, p);
679  else
680  break;
681  }
682 
683  case 2:
684  {
685  const Point<dim> point =
687  p);
688 
689  // formula not guaranteed to work for points outside of
690  // the cell. only take the computed point if it lies
691  // inside the reference cell
692  const double eps = 1e-15;
693  if (-eps <= point(1) && point(1) <= 1 + eps &&
694  -eps <= point(0) && point(0) <= 1 + eps)
695  {
696  return point;
697  }
698  else
699  break;
700  }
701 
702  default:
703  {
704  // we should get here, based on the if-condition at the top
705  Assert(false, ExcInternalError());
706  }
707  }
708  }
709  catch (
711  {
712  // simply fall through and continue on to the standard Newton code
713  }
714  }
715  else
716  {
717  // we can't use an explicit formula,
718  }
719 
720 
721  // Find the initial value for the Newton iteration by a normal
722  // projection to the least square plane determined by the vertices
723  // of the cell
724  Point<dim> initial_p_unit;
725  if (this->preserves_vertex_locations())
726  {
727  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
728  // in 1d with spacedim > 1 the affine approximation is exact
729  if (dim == 1 && polynomial_degree == 1)
730  return initial_p_unit;
731  }
732  else
733  {
734  // else, we simply use the mid point
735  for (unsigned int d = 0; d < dim; ++d)
736  initial_p_unit[d] = 0.5;
737  }
738 
739  // perform the Newton iteration and return the result. note that this
740  // statement may throw an exception, which we simply pass up to the caller
741  const Point<dim> p_unit =
742  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
743  if (p_unit[0] == std::numeric_limits<double>::infinity())
744  AssertThrow(false,
746  return p_unit;
747 }
748 
749 
750 
751 template <int dim, int spacedim>
752 void
754  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
755  const ArrayView<const Point<spacedim>> & real_points,
756  const ArrayView<Point<dim>> & unit_points) const
757 {
758  // Go to base class functions for dim < spacedim because it is not yet
759  // implemented with optimized code.
760  if (dim < spacedim)
761  {
763  real_points,
764  unit_points);
765  return;
766  }
767 
768  AssertDimension(real_points.size(), unit_points.size());
769  const std::vector<Point<spacedim>> support_points =
770  this->compute_mapping_support_points(cell);
771 
772  // From the given (high-order) support points, now only pick the first
773  // 2^dim points and construct an affine approximation from those.
775  inverse_approximation(support_points, unit_cell_support_points);
776 
777  const unsigned int n_points = real_points.size();
778  const unsigned int n_lanes = VectorizedArray<double>::size();
779 
780  // Use the more heavy VectorizedArray code path if there is more than
781  // one point left to compute
782  for (unsigned int i = 0; i < n_points; i += n_lanes)
783  if (n_points - i > 1)
784  {
786  for (unsigned int j = 0; j < n_lanes; ++j)
787  if (i + j < n_points)
788  for (unsigned int d = 0; d < spacedim; ++d)
789  p_vec[d][j] = real_points[i + j][d];
790  else
791  for (unsigned int d = 0; d < spacedim; ++d)
792  p_vec[d][j] = real_points[i][d];
793 
795  internal::MappingQImplementation::
796  do_transform_real_to_unit_cell_internal<dim, spacedim>(
797  p_vec,
798  inverse_approximation.compute(p_vec),
799  support_points,
800  polynomials_1d,
801  renumber_lexicographic_to_hierarchic);
802 
803  // If the vectorized computation failed, it could be that only some of
804  // the lanes failed but others would have succeeded if we had let them
805  // compute alone without interference (like negative Jacobian
806  // determinants) from other SIMD lanes. Repeat the computation in this
807  // unlikely case with scalar arguments.
808  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
809  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
810  unit_points[i + j] = internal::MappingQImplementation::
811  do_transform_real_to_unit_cell_internal<dim, spacedim>(
812  real_points[i + j],
813  inverse_approximation.compute(real_points[i + j]),
814  support_points,
815  polynomials_1d,
816  renumber_lexicographic_to_hierarchic);
817  else
818  for (unsigned int d = 0; d < dim; ++d)
819  unit_points[i + j][d] = unit_point[d][j];
820  }
821  else
822  unit_points[i] = internal::MappingQImplementation::
823  do_transform_real_to_unit_cell_internal<dim, spacedim>(
824  real_points[i],
825  inverse_approximation.compute(real_points[i]),
826  support_points,
827  polynomials_1d,
828  renumber_lexicographic_to_hierarchic);
829 }
830 
831 
832 
833 template <int dim, int spacedim>
836 {
837  // add flags if the respective quantities are necessary to compute
838  // what we need. note that some flags appear in both the conditions
839  // and in subsequent set operations. this leads to some circular
840  // logic. the only way to treat this is to iterate. since there are
841  // 5 if-clauses in the loop, it will take at most 5 iterations to
842  // converge. do them:
843  UpdateFlags out = in;
844  for (unsigned int i = 0; i < 5; ++i)
845  {
846  // The following is a little incorrect:
847  // If not applied on a face,
848  // update_boundary_forms does not
849  // make sense. On the other hand,
850  // it is necessary on a
851  // face. Currently,
852  // update_boundary_forms is simply
853  // ignored for the interior of a
854  // cell.
856  out |= update_boundary_forms;
857 
862 
863  if (out &
868 
869  // The contravariant transformation is used in the Piola
870  // transformation, which requires the determinant of the Jacobi
871  // matrix of the transformation. Because we have no way of
872  // knowing here whether the finite element wants to use the
873  // contravariant or the Piola transforms, we add the JxW values
874  // to the list of flags to be updated for each cell.
876  out |= update_volume_elements;
877 
878  // the same is true when computing normal vectors: they require
879  // the determinant of the Jacobian
880  if (out & update_normal_vectors)
881  out |= update_volume_elements;
882  }
883 
884  return out;
885 }
886 
887 
888 
889 template <int dim, int spacedim>
890 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
892  const Quadrature<dim> &q) const
893 {
894  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
895  std::make_unique<InternalData>(polynomial_degree);
896  auto &data = dynamic_cast<InternalData &>(*data_ptr);
897  data.initialize(this->requires_update_flags(update_flags), q, q.size());
898 
899  return data_ptr;
900 }
901 
902 
903 
904 template <int dim, int spacedim>
905 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
907  const UpdateFlags update_flags,
908  const hp::QCollection<dim - 1> &quadrature) const
909 {
910  AssertDimension(quadrature.size(), 1);
911 
912  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
913  std::make_unique<InternalData>(polynomial_degree);
914  auto &data = dynamic_cast<InternalData &>(*data_ptr);
915  data.initialize_face(this->requires_update_flags(update_flags),
917  ReferenceCells::get_hypercube<dim>(), quadrature[0]),
918  quadrature[0].size());
919 
920  return data_ptr;
921 }
922 
923 
924 
925 template <int dim, int spacedim>
926 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
928  const UpdateFlags update_flags,
929  const Quadrature<dim - 1> &quadrature) const
930 {
931  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
932  std::make_unique<InternalData>(polynomial_degree);
933  auto &data = dynamic_cast<InternalData &>(*data_ptr);
934  data.initialize_face(this->requires_update_flags(update_flags),
936  ReferenceCells::get_hypercube<dim>(), quadrature),
937  quadrature.size());
938 
939  return data_ptr;
940 }
941 
942 
943 
944 template <int dim, int spacedim>
947  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
948  const CellSimilarity::Similarity cell_similarity,
949  const Quadrature<dim> & quadrature,
950  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
952  &output_data) const
953 {
954  // ensure that the following static_cast is really correct:
955  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
956  ExcInternalError());
957  const InternalData &data = static_cast<const InternalData &>(internal_data);
958 
959  const unsigned int n_q_points = quadrature.size();
960 
961  // recompute the support points of the transformation of this
962  // cell. we tried to be clever here in an earlier version of the
963  // library by checking whether the cell is the same as the one we
964  // had visited last, but it turns out to be difficult to determine
965  // that because a cell for the purposes of a mapping is
966  // characterized not just by its (triangulation, level, index)
967  // triple, but also by the locations of its vertices, the manifold
968  // object attached to the cell and all of its bounding faces/edges,
969  // etc. to reliably test that the "cell" we are on is, therefore,
970  // not easily done
971  data.mapping_support_points = this->compute_mapping_support_points(cell);
972  data.cell_of_current_support_points = cell;
973 
974  // if the order of the mapping is greater than 1, then do not reuse any cell
975  // similarity information. This is necessary because the cell similarity
976  // value is computed with just cell vertices and does not take into account
977  // cell curvature.
978  const CellSimilarity::Similarity computed_cell_similarity =
979  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
980 
981  if (dim > 1 && data.tensor_product_quadrature)
982  {
983  internal::MappingQImplementation::
984  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
985  computed_cell_similarity,
986  data,
987  output_data.quadrature_points,
988  output_data.jacobian_grads);
989  }
990  else
991  {
992  internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
994  data,
995  output_data.quadrature_points);
996 
997  internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
998  computed_cell_similarity,
1000  data);
1001 
1003  spacedim>(
1004  computed_cell_similarity,
1006  data,
1007  output_data.jacobian_grads);
1008  }
1009 
1011  dim,
1012  spacedim>(computed_cell_similarity,
1014  data,
1015  output_data.jacobian_pushed_forward_grads);
1016 
1018  dim,
1019  spacedim>(computed_cell_similarity,
1021  data,
1022  output_data.jacobian_2nd_derivatives);
1023 
1024  internal::MappingQImplementation::
1025  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1026  computed_cell_similarity,
1028  data,
1030 
1032  dim,
1033  spacedim>(computed_cell_similarity,
1035  data,
1036  output_data.jacobian_3rd_derivatives);
1037 
1038  internal::MappingQImplementation::
1039  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1040  computed_cell_similarity,
1042  data,
1044 
1045  const UpdateFlags update_flags = data.update_each;
1046  const std::vector<double> &weights = quadrature.get_weights();
1047 
1048  // Multiply quadrature weights by absolute value of Jacobian determinants or
1049  // the area element g=sqrt(DX^t DX) in case of codim > 0
1050 
1051  if (update_flags & (update_normal_vectors | update_JxW_values))
1052  {
1053  AssertDimension(output_data.JxW_values.size(), n_q_points);
1054 
1055  Assert(!(update_flags & update_normal_vectors) ||
1056  (output_data.normal_vectors.size() == n_q_points),
1057  ExcDimensionMismatch(output_data.normal_vectors.size(),
1058  n_q_points));
1059 
1060 
1061  if (computed_cell_similarity != CellSimilarity::translation)
1062  for (unsigned int point = 0; point < n_q_points; ++point)
1063  {
1064  if (dim == spacedim)
1065  {
1066  const double det = data.contravariant[point].determinant();
1067 
1068  // check for distorted cells.
1069 
1070  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1071  // 1e12 in 2D. might want to find a finer
1072  // (dimension-independent) criterion
1073  Assert(det >
1074  1e-12 * Utilities::fixed_power<dim>(
1075  cell->diameter() / std::sqrt(double(dim))),
1077  cell->center(), det, point)));
1078 
1079  output_data.JxW_values[point] = weights[point] * det;
1080  }
1081  // if dim==spacedim, then there is no cell normal to
1082  // compute. since this is for FEValues (and not FEFaceValues),
1083  // there are also no face normals to compute
1084  else // codim>0 case
1085  {
1086  Tensor<1, spacedim> DX_t[dim];
1087  for (unsigned int i = 0; i < spacedim; ++i)
1088  for (unsigned int j = 0; j < dim; ++j)
1089  DX_t[j][i] = data.contravariant[point][i][j];
1090 
1091  Tensor<2, dim> G; // First fundamental form
1092  for (unsigned int i = 0; i < dim; ++i)
1093  for (unsigned int j = 0; j < dim; ++j)
1094  G[i][j] = DX_t[i] * DX_t[j];
1095 
1096  output_data.JxW_values[point] =
1097  std::sqrt(determinant(G)) * weights[point];
1098 
1099  if (computed_cell_similarity ==
1101  {
1102  // we only need to flip the normal
1103  if (update_flags & update_normal_vectors)
1104  output_data.normal_vectors[point] *= -1.;
1105  }
1106  else
1107  {
1108  if (update_flags & update_normal_vectors)
1109  {
1110  Assert(spacedim == dim + 1,
1111  ExcMessage(
1112  "There is no (unique) cell normal for " +
1114  "-dimensional cells in " +
1115  Utilities::int_to_string(spacedim) +
1116  "-dimensional space. This only works if the "
1117  "space dimension is one greater than the "
1118  "dimensionality of the mesh cells."));
1119 
1120  if (dim == 1)
1121  output_data.normal_vectors[point] =
1122  cross_product_2d(-DX_t[0]);
1123  else // dim == 2
1124  output_data.normal_vectors[point] =
1125  cross_product_3d(DX_t[0], DX_t[1]);
1126 
1127  output_data.normal_vectors[point] /=
1128  output_data.normal_vectors[point].norm();
1129 
1130  if (cell->direction_flag() == false)
1131  output_data.normal_vectors[point] *= -1.;
1132  }
1133  }
1134  } // codim>0 case
1135  }
1136  }
1137 
1138 
1139 
1140  // copy values from InternalData to vector given by reference
1141  if (update_flags & update_jacobians)
1142  {
1143  AssertDimension(output_data.jacobians.size(), n_q_points);
1144  if (computed_cell_similarity != CellSimilarity::translation)
1145  for (unsigned int point = 0; point < n_q_points; ++point)
1146  output_data.jacobians[point] = data.contravariant[point];
1147  }
1148 
1149  // copy values from InternalData to vector given by reference
1150  if (update_flags & update_inverse_jacobians)
1151  {
1152  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1153  if (computed_cell_similarity != CellSimilarity::translation)
1154  for (unsigned int point = 0; point < n_q_points; ++point)
1155  output_data.inverse_jacobians[point] =
1156  data.covariant[point].transpose();
1157  }
1158 
1159  return computed_cell_similarity;
1160 }
1161 
1162 
1163 
1164 template <int dim, int spacedim>
1165 void
1167  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1168  const unsigned int face_no,
1169  const hp::QCollection<dim - 1> & quadrature,
1170  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1172  &output_data) const
1173 {
1174  AssertDimension(quadrature.size(), 1);
1175 
1176  // ensure that the following cast is really correct:
1177  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1178  ExcInternalError());
1179  const InternalData &data = static_cast<const InternalData &>(internal_data);
1180 
1181  // if necessary, recompute the support points of the transformation of this
1182  // cell (note that we need to first check the triangulation pointer, since
1183  // otherwise the second test might trigger an exception if the triangulations
1184  // are not the same)
1185  if ((data.mapping_support_points.size() == 0) ||
1186  (&cell->get_triangulation() !=
1188  (cell != data.cell_of_current_support_points))
1189  {
1190  data.mapping_support_points = this->compute_mapping_support_points(cell);
1191  data.cell_of_current_support_points = cell;
1192  }
1193 
1195  *this,
1196  cell,
1197  face_no,
1200  ReferenceCells::get_hypercube<dim>(),
1201  face_no,
1202  cell->face_orientation(face_no),
1203  cell->face_flip(face_no),
1204  cell->face_rotation(face_no),
1205  quadrature[0].size()),
1206  quadrature[0],
1207  data,
1208  output_data);
1209 }
1210 
1211 
1212 
1213 template <int dim, int spacedim>
1214 void
1216  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1217  const unsigned int face_no,
1218  const unsigned int subface_no,
1219  const Quadrature<dim - 1> & quadrature,
1220  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1222  &output_data) const
1223 {
1224  // ensure that the following cast is really correct:
1225  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1226  ExcInternalError());
1227  const InternalData &data = static_cast<const InternalData &>(internal_data);
1228 
1229  // if necessary, recompute the support points of the transformation of this
1230  // cell (note that we need to first check the triangulation pointer, since
1231  // otherwise the second test might trigger an exception if the triangulations
1232  // are not the same)
1233  if ((data.mapping_support_points.size() == 0) ||
1234  (&cell->get_triangulation() !=
1236  (cell != data.cell_of_current_support_points))
1237  {
1238  data.mapping_support_points = this->compute_mapping_support_points(cell);
1239  data.cell_of_current_support_points = cell;
1240  }
1241 
1243  *this,
1244  cell,
1245  face_no,
1246  subface_no,
1248  ReferenceCells::get_hypercube<dim>(),
1249  face_no,
1250  subface_no,
1251  cell->face_orientation(face_no),
1252  cell->face_flip(face_no),
1253  cell->face_rotation(face_no),
1254  quadrature.size(),
1255  cell->subface_case(face_no)),
1256  quadrature,
1257  data,
1258  output_data);
1259 }
1260 
1261 
1262 
1263 template <int dim, int spacedim>
1264 void
1266  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1268  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1270  &output_data) const
1271 {
1272  Assert(dim == spacedim, ExcNotImplemented());
1273 
1274  // ensure that the following static_cast is really correct:
1275  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1276  ExcInternalError());
1277  const InternalData &data = static_cast<const InternalData &>(internal_data);
1278 
1279  const unsigned int n_q_points = quadrature.size();
1280 
1281  data.mapping_support_points = this->compute_mapping_support_points(cell);
1282  data.cell_of_current_support_points = cell;
1283 
1284  internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
1286  data,
1287  output_data.quadrature_points);
1288 
1289  internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
1291 
1292  internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1295  data,
1296  output_data.jacobian_grads);
1297 
1299  dim,
1300  spacedim>(CellSimilarity::none,
1302  data,
1303  output_data.jacobian_pushed_forward_grads);
1304 
1306  dim,
1307  spacedim>(CellSimilarity::none,
1309  data,
1310  output_data.jacobian_2nd_derivatives);
1311 
1312  internal::MappingQImplementation::
1313  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1316  data,
1318 
1320  dim,
1321  spacedim>(CellSimilarity::none,
1323  data,
1324  output_data.jacobian_3rd_derivatives);
1325 
1326  internal::MappingQImplementation::
1327  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1330  data,
1332 
1333  const UpdateFlags update_flags = data.update_each;
1334  const std::vector<double> &weights = quadrature.get_weights();
1335 
1336  if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1337  {
1338  AssertDimension(output_data.JxW_values.size(), n_q_points);
1339 
1340  Assert(!(update_flags & update_normal_vectors) ||
1341  (output_data.normal_vectors.size() == n_q_points),
1342  ExcDimensionMismatch(output_data.normal_vectors.size(),
1343  n_q_points));
1344 
1345 
1346  for (unsigned int point = 0; point < n_q_points; ++point)
1347  {
1348  const double det = data.contravariant[point].determinant();
1349 
1350  // check for distorted cells.
1351 
1352  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1353  // 1e12 in 2D. might want to find a finer
1354  // (dimension-independent) criterion
1355  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1356  cell->diameter() / std::sqrt(double(dim))),
1358  cell->center(), det, point)));
1359 
1360  // The normals are n = J^{-T} * \hat{n} before normalizing.
1361  Tensor<1, spacedim> normal;
1362  for (unsigned int d = 0; d < spacedim; d++)
1363  normal[d] =
1364  data.covariant[point][d] * quadrature.normal_vector(point);
1365 
1366  output_data.JxW_values[point] = weights[point] * det * normal.norm();
1367 
1368  if ((update_flags & update_normal_vectors) != 0u)
1369  {
1370  normal /= normal.norm();
1371  output_data.normal_vectors[point] = normal;
1372  }
1373  }
1374  }
1375 
1376  // copy values from InternalData to vector given by reference
1377  if ((update_flags & update_jacobians) != 0u)
1378  {
1379  AssertDimension(output_data.jacobians.size(), n_q_points);
1380  for (unsigned int point = 0; point < n_q_points; ++point)
1381  output_data.jacobians[point] = data.contravariant[point];
1382  }
1383 
1384  // copy values from InternalData to vector given by reference
1385  if ((update_flags & update_inverse_jacobians) != 0u)
1386  {
1387  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1388  for (unsigned int point = 0; point < n_q_points; ++point)
1389  output_data.inverse_jacobians[point] =
1390  data.covariant[point].transpose();
1391  }
1392 }
1393 
1394 
1395 
1396 template <int dim, int spacedim>
1397 void
1399  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1400  const ArrayView<const Point<dim>> & unit_points,
1401  const UpdateFlags update_flags,
1403  &output_data) const
1404 {
1405  if (update_flags == update_default)
1406  return;
1407 
1408  Assert(update_flags & update_inverse_jacobians ||
1409  update_flags & update_jacobians ||
1410  update_flags & update_quadrature_points,
1411  ExcNotImplemented());
1412 
1413  output_data.initialize(unit_points.size(), update_flags);
1414  const std::vector<Point<spacedim>> support_points =
1415  this->compute_mapping_support_points(cell);
1416 
1417  const unsigned int n_points = unit_points.size();
1418  const unsigned int n_lanes = VectorizedArray<double>::size();
1419 
1420  // Use the more heavy VectorizedArray code path if there is more than
1421  // one point left to compute
1422  for (unsigned int i = 0; i < n_points; i += n_lanes)
1423  if (n_points - i > 1)
1424  {
1426  for (unsigned int j = 0; j < n_lanes; ++j)
1427  if (i + j < n_points)
1428  for (unsigned int d = 0; d < dim; ++d)
1429  p_vec[d][j] = unit_points[i + j][d];
1430  else
1431  for (unsigned int d = 0; d < dim; ++d)
1432  p_vec[d][j] = unit_points[i][d];
1433 
1434  const auto result =
1436  polynomials_1d,
1437  support_points,
1438  p_vec,
1439  polynomial_degree == 1,
1440  renumber_lexicographic_to_hierarchic);
1441 
1442  if (update_flags & update_quadrature_points)
1443  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1444  for (unsigned int d = 0; d < spacedim; ++d)
1445  output_data.quadrature_points[i + j][d] = result.first[d][j];
1446 
1447  if (update_flags & update_jacobians)
1448  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1449  for (unsigned int d = 0; d < spacedim; ++d)
1450  for (unsigned int e = 0; e < dim; ++e)
1451  output_data.jacobians[i + j][d][e] = result.second[e][d][j];
1452 
1453  if (update_flags & update_inverse_jacobians)
1454  {
1456  result.second);
1458  inv_jac = jac.covariant_form();
1459  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1460  for (unsigned int d = 0; d < dim; ++d)
1461  for (unsigned int e = 0; e < spacedim; ++e)
1462  output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j];
1463  }
1464  }
1465  else
1466  {
1467  const auto result =
1469  polynomials_1d,
1470  support_points,
1471  unit_points[i],
1472  polynomial_degree == 1,
1473  renumber_lexicographic_to_hierarchic);
1474 
1475  if (update_flags & update_quadrature_points)
1476  output_data.quadrature_points[i] = result.first;
1477 
1478  if (update_flags & update_jacobians)
1479  {
1480  DerivativeForm<1, spacedim, dim> jac = result.second;
1481  output_data.jacobians[i] = jac.transpose();
1482  }
1483 
1484  if (update_flags & update_inverse_jacobians)
1485  {
1486  DerivativeForm<1, spacedim, dim> jac(result.second);
1488  for (unsigned int d = 0; d < dim; ++d)
1489  for (unsigned int e = 0; e < spacedim; ++e)
1490  output_data.inverse_jacobians[i][d][e] = inv_jac[d][e];
1491  }
1492  }
1493 }
1494 
1495 
1496 
1497 template <int dim, int spacedim>
1498 void
1500  const ArrayView<const Tensor<1, dim>> & input,
1501  const MappingKind mapping_kind,
1502  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1503  const ArrayView<Tensor<1, spacedim>> & output) const
1504 {
1506  mapping_kind,
1507  mapping_data,
1508  output);
1509 }
1510 
1511 
1512 
1513 template <int dim, int spacedim>
1514 void
1516  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1517  const MappingKind mapping_kind,
1518  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1519  const ArrayView<Tensor<2, spacedim>> & output) const
1520 {
1522  mapping_kind,
1523  mapping_data,
1524  output);
1525 }
1526 
1527 
1528 
1529 template <int dim, int spacedim>
1530 void
1532  const ArrayView<const Tensor<2, dim>> & input,
1533  const MappingKind mapping_kind,
1534  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1535  const ArrayView<Tensor<2, spacedim>> & output) const
1536 {
1537  switch (mapping_kind)
1538  {
1539  case mapping_contravariant:
1541  mapping_kind,
1542  mapping_data,
1543  output);
1544  return;
1545 
1550  mapping_kind,
1551  mapping_data,
1552  output);
1553  return;
1554  default:
1555  Assert(false, ExcNotImplemented());
1556  }
1557 }
1558 
1559 
1560 
1561 template <int dim, int spacedim>
1562 void
1564  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1565  const MappingKind mapping_kind,
1566  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1567  const ArrayView<Tensor<3, spacedim>> & output) const
1568 {
1569  AssertDimension(input.size(), output.size());
1570  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1571  ExcInternalError());
1572  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1573 
1574  switch (mapping_kind)
1575  {
1577  {
1578  Assert(data.update_each & update_contravariant_transformation,
1580  "update_covariant_transformation"));
1581 
1582  for (unsigned int q = 0; q < output.size(); ++q)
1583  for (unsigned int i = 0; i < spacedim; ++i)
1584  for (unsigned int j = 0; j < spacedim; ++j)
1585  {
1586  double tmp[dim];
1587  for (unsigned int K = 0; K < dim; ++K)
1588  {
1589  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1590  for (unsigned int J = 1; J < dim; ++J)
1591  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1592  }
1593  for (unsigned int k = 0; k < spacedim; ++k)
1594  {
1595  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1596  for (unsigned int K = 1; K < dim; ++K)
1597  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1598  }
1599  }
1600  return;
1601  }
1602 
1603  default:
1604  Assert(false, ExcNotImplemented());
1605  }
1606 }
1607 
1608 
1609 
1610 template <int dim, int spacedim>
1611 void
1613  const ArrayView<const Tensor<3, dim>> & input,
1614  const MappingKind mapping_kind,
1615  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1616  const ArrayView<Tensor<3, spacedim>> & output) const
1617 {
1618  switch (mapping_kind)
1619  {
1620  case mapping_piola_hessian:
1624  mapping_kind,
1625  mapping_data,
1626  output);
1627  return;
1628  default:
1629  Assert(false, ExcNotImplemented());
1630  }
1631 }
1632 
1633 
1634 
1635 template <int dim, int spacedim>
1636 void
1638  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1639  std::vector<Point<spacedim>> & a) const
1640 {
1641  // if we only need the midpoint, then ask for it.
1642  if (this->polynomial_degree == 2)
1643  {
1644  for (unsigned int line_no = 0;
1645  line_no < GeometryInfo<dim>::lines_per_cell;
1646  ++line_no)
1647  {
1648  const typename Triangulation<dim, spacedim>::line_iterator line =
1649  (dim == 1 ?
1650  static_cast<
1652  cell->line(line_no));
1653 
1654  const Manifold<dim, spacedim> &manifold =
1655  ((line->manifold_id() == numbers::flat_manifold_id) &&
1656  (dim < spacedim) ?
1657  cell->get_manifold() :
1658  line->get_manifold());
1659  a.push_back(manifold.get_new_point_on_line(line));
1660  }
1661  }
1662  else
1663  // otherwise call the more complicated functions and ask for inner points
1664  // from the manifold description
1665  {
1666  std::vector<Point<spacedim>> tmp_points;
1667  for (unsigned int line_no = 0;
1668  line_no < GeometryInfo<dim>::lines_per_cell;
1669  ++line_no)
1670  {
1671  const typename Triangulation<dim, spacedim>::line_iterator line =
1672  (dim == 1 ?
1673  static_cast<
1675  cell->line(line_no));
1676 
1677  const Manifold<dim, spacedim> &manifold =
1678  ((line->manifold_id() == numbers::flat_manifold_id) &&
1679  (dim < spacedim) ?
1680  cell->get_manifold() :
1681  line->get_manifold());
1682 
1683  const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1684  const std::array<Point<spacedim>, 2> vertices{
1685  {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
1686  cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
1687 
1688  const std::size_t n_rows =
1689  support_point_weights_perimeter_to_interior[0].size(0);
1690  a.resize(a.size() + n_rows);
1691  auto a_view = make_array_view(a.end() - n_rows, a.end());
1692  manifold.get_new_points(
1693  make_array_view(vertices.begin(), vertices.end()),
1694  support_point_weights_perimeter_to_interior[0],
1695  a_view);
1696  }
1697  }
1698 }
1699 
1700 
1701 
1702 template <>
1703 void
1706  std::vector<Point<3>> & a) const
1707 {
1708  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1709 
1710  // used if face quad at boundary or entirely in the interior of the domain
1711  std::vector<Point<3>> tmp_points;
1712 
1713  // loop over all faces and collect points on them
1714  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1715  {
1716  const Triangulation<3>::face_iterator face = cell->face(face_no);
1717 
1718 #ifdef DEBUG
1719  const bool face_orientation = cell->face_orientation(face_no),
1720  face_flip = cell->face_flip(face_no),
1721  face_rotation = cell->face_rotation(face_no);
1722  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1723  lines_per_face = GeometryInfo<3>::lines_per_face;
1724 
1725  // some sanity checks up front
1726  for (unsigned int i = 0; i < vertices_per_face; ++i)
1727  Assert(face->vertex_index(i) ==
1728  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1729  face_no, i, face_orientation, face_flip, face_rotation)),
1730  ExcInternalError());
1731 
1732  // indices of the lines that bound a face are given by GeometryInfo<3>::
1733  // face_to_cell_lines
1734  for (unsigned int i = 0; i < lines_per_face; ++i)
1735  Assert(face->line(i) ==
1737  face_no, i, face_orientation, face_flip, face_rotation)),
1738  ExcInternalError());
1739 #endif
1740  // extract the points surrounding a quad from the points
1741  // already computed. First get the 4 vertices and then the points on
1742  // the four lines
1743  boost::container::small_vector<Point<3>, 200> tmp_points(
1745  GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1746  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1747  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1748  if (polynomial_degree > 1)
1749  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1750  ++line)
1751  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1752  tmp_points[4 + line * (polynomial_degree - 1) + i] =
1754  (polynomial_degree - 1) *
1755  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
1756  i];
1757 
1758  const std::size_t n_rows =
1759  support_point_weights_perimeter_to_interior[1].size(0);
1760  a.resize(a.size() + n_rows);
1761  auto a_view = make_array_view(a.end() - n_rows, a.end());
1762  face->get_manifold().get_new_points(
1763  make_array_view(tmp_points.begin(), tmp_points.end()),
1764  support_point_weights_perimeter_to_interior[1],
1765  a_view);
1766  }
1767 }
1768 
1769 
1770 
1771 template <>
1772 void
1775  std::vector<Point<3>> & a) const
1776 {
1777  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1778  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1779  vertices[i] = cell->vertex(i);
1780 
1781  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1783  for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1784  for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1785  {
1786  Point<2> point(line_support_points[q1 + 1][0],
1787  line_support_points[q2 + 1][0]);
1788  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1789  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1790  }
1791 
1792  const std::size_t n_rows = weights.size(0);
1793  a.resize(a.size() + n_rows);
1794  auto a_view = make_array_view(a.end() - n_rows, a.end());
1795  cell->get_manifold().get_new_points(
1796  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1797 }
1798 
1799 
1800 
1801 template <int dim, int spacedim>
1802 void
1805  std::vector<Point<spacedim>> &) const
1806 {
1807  Assert(false, ExcInternalError());
1808 }
1809 
1810 
1811 
1812 template <int dim, int spacedim>
1813 std::vector<Point<spacedim>>
1815  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1816 {
1817  // get the vertices first
1818  std::vector<Point<spacedim>> a;
1819  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1820  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1821  a.push_back(cell->vertex(i));
1822 
1823  if (this->polynomial_degree > 1)
1824  {
1825  // check if all entities have the same manifold id which is when we can
1826  // simply ask the manifold for all points. the transfinite manifold can
1827  // do the interpolation better than this class, so if we detect that we
1828  // do not have to change anything here
1829  Assert(dim <= 3, ExcImpossibleInDim(dim));
1830  bool all_manifold_ids_are_equal = (dim == spacedim);
1831  if (all_manifold_ids_are_equal &&
1833  &cell->get_manifold()) == nullptr)
1834  {
1835  for (auto f : GeometryInfo<dim>::face_indices())
1836  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1837  all_manifold_ids_are_equal = false;
1838 
1839  if (dim == 3)
1840  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1841  if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1842  all_manifold_ids_are_equal = false;
1843  }
1844 
1845  if (all_manifold_ids_are_equal)
1846  {
1847  const std::size_t n_rows = support_point_weights_cell.size(0);
1848  a.resize(a.size() + n_rows);
1849  auto a_view = make_array_view(a.end() - n_rows, a.end());
1850  cell->get_manifold().get_new_points(make_array_view(a.begin(),
1851  a.end() - n_rows),
1852  support_point_weights_cell,
1853  a_view);
1854  }
1855  else
1856  switch (dim)
1857  {
1858  case 1:
1859  add_line_support_points(cell, a);
1860  break;
1861  case 2:
1862  // in 2d, add the points on the four bounding lines to the
1863  // exterior (outer) points
1864  add_line_support_points(cell, a);
1865 
1866  // then get the interior support points
1867  if (dim != spacedim)
1868  add_quad_support_points(cell, a);
1869  else
1870  {
1871  const std::size_t n_rows =
1872  support_point_weights_perimeter_to_interior[1].size(0);
1873  a.resize(a.size() + n_rows);
1874  auto a_view = make_array_view(a.end() - n_rows, a.end());
1875  cell->get_manifold().get_new_points(
1876  make_array_view(a.begin(), a.end() - n_rows),
1877  support_point_weights_perimeter_to_interior[1],
1878  a_view);
1879  }
1880  break;
1881 
1882  case 3:
1883  // in 3d also add the points located on the boundary faces
1884  add_line_support_points(cell, a);
1885  add_quad_support_points(cell, a);
1886 
1887  // then compute the interior points
1888  {
1889  const std::size_t n_rows =
1890  support_point_weights_perimeter_to_interior[2].size(0);
1891  a.resize(a.size() + n_rows);
1892  auto a_view = make_array_view(a.end() - n_rows, a.end());
1893  cell->get_manifold().get_new_points(
1894  make_array_view(a.begin(), a.end() - n_rows),
1895  support_point_weights_perimeter_to_interior[2],
1896  a_view);
1897  }
1898  break;
1899 
1900  default:
1901  Assert(false, ExcNotImplemented());
1902  break;
1903  }
1904  }
1905 
1906  return a;
1907 }
1908 
1909 
1910 
1911 template <int dim, int spacedim>
1914  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1915 {
1916  return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1917 }
1918 
1919 
1920 
1921 template <int dim, int spacedim>
1922 bool
1924  const ReferenceCell &reference_cell) const
1925 {
1926  Assert(dim == reference_cell.get_dimension(),
1927  ExcMessage("The dimension of your mapping (" +
1928  Utilities::to_string(dim) +
1929  ") and the reference cell cell_type (" +
1930  Utilities::to_string(reference_cell.get_dimension()) +
1931  " ) do not agree."));
1932 
1933  return reference_cell.is_hyper_cube();
1934 }
1935 
1936 
1937 
1938 //--------------------------- Explicit instantiations -----------------------
1939 #include "mapping_q.inst"
1940 
1941 
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:685
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
DerivativeForm< 1, spacedim, dim, Number > transpose() const
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:125
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:354
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition: mapping_q.h:548
AlignedVector< DerivativeForm< 1, dim, spacedim > > covariant
Definition: mapping_q.h:522
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:65
std::vector< Point< spacedim > > mapping_support_points
Definition: mapping_q.h:542
AlignedVector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition: mapping_q.h:531
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:216
InternalData(const unsigned int polynomial_degree)
Definition: mapping_q.cc:53
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
Definition: mapping_q.cc:272
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:86
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
Definition: mapping_q.cc:1637
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition: mapping_q.h:652
const Table< 2, double > support_point_weights_cell
Definition: mapping_q.h:700
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_q.cc:1814
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:891
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:1398
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_q.cc:1913
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:927
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:1215
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:946
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:1265
const unsigned int polynomial_degree
Definition: mapping_q.h:628
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:1499
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:906
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:489
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition: mapping_q.h:686
const std::vector< Point< 1 > > line_support_points
Definition: mapping_q.h:638
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:454
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition: mapping_q.h:645
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:436
const std::vector< Point< dim > > unit_cell_support_points
Definition: mapping_q.h:664
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:753
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:632
MappingQ(const unsigned int polynomial_degree)
Definition: mapping_q.cc:363
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_q.cc:835
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:1166
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition: mapping_q.cc:1923
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
Definition: mapping_q.cc:1803
unsigned int get_degree() const
Definition: mapping_q.cc:445
Abstract base class for mapping classes.
Definition: mapping.h:312
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
Definition: point.h:111
static DataSetDescriptor cell()
Definition: qprojector.h:361
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:325
unsigned int size() const
constexpr 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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:517
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:1583
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1756
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1672
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_3rd_derivatives
@ 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.
@ update_jacobian_2nd_derivatives
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
MappingKind
Definition: mapping.h:73
@ mapping_covariant_gradient
Definition: mapping.h:94
@ mapping_contravariant
Definition: mapping.h:88
@ mapping_contravariant_hessian
Definition: mapping.h:150
@ mapping_covariant_hessian
Definition: mapping.h:144
@ mapping_contravariant_gradient
Definition: mapping.h:100
@ mapping_piola_gradient
Definition: mapping.h:114
@ mapping_piola_hessian
Definition: mapping.h:156
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
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:450
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
T fixed_power(const T t)
Definition: utilities.h:969
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)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
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 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, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
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 maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
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)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::manifold_id flat_manifold_id
Definition: types.h:279
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)