Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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 
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 {}
58 
59 
60 
61 template <int dim, int spacedim>
62 std::size_t
64 {
65  return (
68  MemoryConsumption::memory_consumption(shape_derivatives) +
71  MemoryConsumption::memory_consumption(unit_tangentials) +
73  MemoryConsumption::memory_consumption(mapping_support_points) +
74  MemoryConsumption::memory_consumption(cell_of_current_support_points) +
75  MemoryConsumption::memory_consumption(volume_elements) +
77  MemoryConsumption::memory_consumption(n_shape_functions));
78 }
79 
80 
81 
82 template <int dim, int spacedim>
83 void
85  const UpdateFlags update_flags,
86  const Quadrature<dim> &q,
87  const unsigned int n_original_q_points)
88 {
89  // store the flags in the internal data object so we can access them
90  // in fill_fe_*_values()
91  this->update_each = update_flags;
92 
93  const unsigned int n_q_points = q.size();
94 
95  const bool needs_higher_order_terms =
96  this->update_each &
101 
102  if (this->update_each & update_covariant_transformation)
103  covariant.resize(n_original_q_points);
104 
105  if (this->update_each & update_contravariant_transformation)
106  contravariant.resize(n_original_q_points);
107 
108  if (this->update_each & update_volume_elements)
109  volume_elements.resize(n_original_q_points);
110 
111  tensor_product_quadrature = q.is_tensor_product();
112 
113  // use of MatrixFree only for higher order elements and with more than one
114  // point where tensor products do not make sense
115  if (polynomial_degree < 2 || n_q_points == 1)
116  tensor_product_quadrature = false;
117 
118  if (dim > 1)
119  {
120  // find out if the one-dimensional formula is the same
121  // in all directions
122  if (tensor_product_quadrature)
123  {
124  const std::array<Quadrature<1>, dim> &quad_array =
125  q.get_tensor_basis();
126  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
127  {
128  if (quad_array[i - 1].size() != quad_array[i].size())
129  {
130  tensor_product_quadrature = false;
131  break;
132  }
133  else
134  {
135  const std::vector<Point<1>> &points_1 =
136  quad_array[i - 1].get_points();
137  const std::vector<Point<1>> &points_2 =
138  quad_array[i].get_points();
139  const std::vector<double> &weights_1 =
140  quad_array[i - 1].get_weights();
141  const std::vector<double> &weights_2 =
142  quad_array[i].get_weights();
143  for (unsigned int j = 0; j < quad_array[i].size(); ++j)
144  {
145  if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
146  std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
147  {
148  tensor_product_quadrature = false;
149  break;
150  }
151  }
152  }
153  }
154 
155  if (tensor_product_quadrature)
156  {
157  // use a 1d FE_DGQ and adjust the hierarchic -> lexicographic
158  // numbering manually (building an FE_Q<dim> is relatively
159  // expensive due to constraints)
160  const FE_DGQ<1> fe(polynomial_degree);
161  shape_info.reinit(q.get_tensor_basis()[0], fe);
162  shape_info.lexicographic_numbering =
163  FETools::lexicographic_to_hierarchic_numbering<dim>(
165  shape_info.n_q_points = q.size();
166  shape_info.dofs_per_component_on_cell =
168  }
169  }
170  }
171 
172  // Only fill the big arrays on demand in case we cannot use the tensor
173  // product quadrature code path
174  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
175  {
176  // see if we need the (transformation) shape function values
177  // and/or gradients and resize the necessary arrays
178  if (this->update_each & update_quadrature_points)
179  shape_values.resize(n_shape_functions * n_q_points);
180 
181  if (this->update_each &
191  shape_derivatives.resize(n_shape_functions * n_q_points);
192 
193  if (this->update_each &
195  shape_second_derivatives.resize(n_shape_functions * n_q_points);
196 
197  if (this->update_each & (update_jacobian_2nd_derivatives |
199  shape_third_derivatives.resize(n_shape_functions * n_q_points);
200 
201  if (this->update_each & (update_jacobian_3rd_derivatives |
203  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
204 
205  // now also fill the various fields with their correct values
206  compute_shape_function_values(q.get_points());
207  }
208 }
209 
210 
211 
212 template <int dim, int spacedim>
213 void
215  const UpdateFlags update_flags,
216  const Quadrature<dim> &q,
217  const unsigned int n_original_q_points)
218 {
219  initialize(update_flags, q, n_original_q_points);
220 
221  if (dim > 1 && tensor_product_quadrature)
222  {
223  constexpr unsigned int facedim = dim - 1;
225  shape_info.reinit(q.get_tensor_basis()[0], fe);
226  shape_info.lexicographic_numbering =
227  FETools::lexicographic_to_hierarchic_numbering<facedim>(
229  shape_info.n_q_points = n_original_q_points;
230  shape_info.dofs_per_component_on_cell =
232  }
233 
234  if (dim > 1)
235  {
236  if (this->update_each &
239  {
240  aux.resize(dim - 1,
241  AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
242 
243  // Compute tangentials to the unit cell.
244  for (const unsigned int i : GeometryInfo<dim>::face_indices())
245  {
246  unit_tangentials[i].resize(n_original_q_points);
247  std::fill(unit_tangentials[i].begin(),
248  unit_tangentials[i].end(),
250  if (dim > 2)
251  {
252  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
253  .resize(n_original_q_points);
254  std::fill(
255  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
256  .begin(),
257  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
258  .end(),
260  }
261  }
262  }
263  }
264 }
265 
266 
267 
268 template <int dim, int spacedim>
269 void
271  const std::vector<Point<dim>> &unit_points)
272 {
273  const unsigned int n_points = unit_points.size();
274 
275  // Construct the tensor product polynomials used as shape functions for
276  // the Qp mapping of cells at the boundary.
277  const TensorProductPolynomials<dim> tensor_pols(
279  line_support_points.get_points()));
280  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
281 
282  // then also construct the mapping from lexicographic to the Qp shape
283  // function numbering
284  const std::vector<unsigned int> renumber =
285  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
286 
287  std::vector<double> values;
288  std::vector<Tensor<1, dim>> grads;
289  if (shape_values.size() != 0)
290  {
291  Assert(shape_values.size() == n_shape_functions * n_points,
293  values.resize(n_shape_functions);
294  }
295  if (shape_derivatives.size() != 0)
296  {
297  Assert(shape_derivatives.size() == n_shape_functions * n_points,
298  ExcInternalError());
299  grads.resize(n_shape_functions);
300  }
301 
302  std::vector<Tensor<2, dim>> grad2;
303  if (shape_second_derivatives.size() != 0)
304  {
305  Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
306  ExcInternalError());
307  grad2.resize(n_shape_functions);
308  }
309 
310  std::vector<Tensor<3, dim>> grad3;
311  if (shape_third_derivatives.size() != 0)
312  {
313  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
314  ExcInternalError());
315  grad3.resize(n_shape_functions);
316  }
317 
318  std::vector<Tensor<4, dim>> grad4;
319  if (shape_fourth_derivatives.size() != 0)
320  {
321  Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
323  grad4.resize(n_shape_functions);
324  }
325 
326 
327  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
328  shape_second_derivatives.size() != 0 ||
329  shape_third_derivatives.size() != 0 ||
330  shape_fourth_derivatives.size() != 0)
331  for (unsigned int point = 0; point < n_points; ++point)
332  {
333  tensor_pols.evaluate(
334  unit_points[point], values, grads, grad2, grad3, grad4);
335 
336  if (shape_values.size() != 0)
337  for (unsigned int i = 0; i < n_shape_functions; ++i)
338  shape(point, i) = values[renumber[i]];
339 
340  if (shape_derivatives.size() != 0)
341  for (unsigned int i = 0; i < n_shape_functions; ++i)
342  derivative(point, i) = grads[renumber[i]];
343 
344  if (shape_second_derivatives.size() != 0)
345  for (unsigned int i = 0; i < n_shape_functions; ++i)
346  second_derivative(point, i) = grad2[renumber[i]];
347 
348  if (shape_third_derivatives.size() != 0)
349  for (unsigned int i = 0; i < n_shape_functions; ++i)
350  third_derivative(point, i) = grad3[renumber[i]];
351 
352  if (shape_fourth_derivatives.size() != 0)
353  for (unsigned int i = 0; i < n_shape_functions; ++i)
354  fourth_derivative(point, i) = grad4[renumber[i]];
355  }
356 }
357 
358 
359 
360 template <int dim, int spacedim>
362  : polynomial_degree(p)
364  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
365  , polynomials_1d(
370  internal::MappingQImplementation::unit_support_points<dim>(
374  internal::MappingQImplementation::
376  this->polynomial_degree,
377  dim))
379  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
380  this->polynomial_degree))
381 {
382  Assert(p >= 1,
383  ExcMessage("It only makes sense to create polynomial mappings "
384  "with a polynomial degree greater or equal to one."));
385 }
386 
387 
388 
389 template <int dim, int spacedim>
390 MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
391  : polynomial_degree(p)
392  , line_support_points(
393  QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
394  , polynomials_1d(
395  Polynomials::generate_complete_Lagrange_basis(line_support_points))
396  , renumber_lexicographic_to_hierarchic(
398  , unit_cell_support_points(
399  internal::MappingQImplementation::unit_support_points<dim>(
400  line_support_points,
401  renumber_lexicographic_to_hierarchic))
402  , support_point_weights_perimeter_to_interior(
403  internal::MappingQImplementation::
405  this->polynomial_degree,
406  dim))
407  , support_point_weights_cell(
408  internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
409  this->polynomial_degree))
410 {
411  Assert(p >= 1,
412  ExcMessage("It only makes sense to create polynomial mappings "
413  "with a polynomial degree greater or equal to one."));
414 }
415 
416 
417 
418 template <int dim, int spacedim>
420  : polynomial_degree(mapping.polynomial_degree)
421  , line_support_points(mapping.line_support_points)
422  , polynomials_1d(mapping.polynomials_1d)
423  , renumber_lexicographic_to_hierarchic(
424  mapping.renumber_lexicographic_to_hierarchic)
425  , support_point_weights_perimeter_to_interior(
426  mapping.support_point_weights_perimeter_to_interior)
427  , support_point_weights_cell(mapping.support_point_weights_cell)
428 {}
429 
430 
431 
432 template <int dim, int spacedim>
433 std::unique_ptr<Mapping<dim, spacedim>>
435 {
436  return std::make_unique<MappingQ<dim, spacedim>>(*this);
437 }
438 
439 
440 
441 template <int dim, int spacedim>
442 unsigned int
444 {
445  return polynomial_degree;
446 }
447 
448 
449 
450 template <int dim, int spacedim>
453  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
454  const Point<dim> & p) const
455 {
457  polynomials_1d,
458  this->compute_mapping_support_points(cell),
459  p,
460  polynomials_1d.size() == 2,
461  renumber_lexicographic_to_hierarchic)
462  .first);
463 }
464 
465 
466 // In the code below, GCC tries to instantiate MappingQ<3,4> when
467 // seeing which of the overloaded versions of
468 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
469 // error messages and, generally, nothing very good. Avoid this by ensuring
470 // that this class exists, but does not have an inner InternalData
471 // type, thereby ruling out the codim-1 version of the function
472 // below when doing overload resolution.
473 template <>
474 class MappingQ<3, 4>
475 {};
476 
477 
478 
479 // visual studio freaks out when trying to determine if
480 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
481 // candidate. So instead of letting the compiler pick the correct overload, we
482 // use template specialization to make sure we pick up the right function to
483 // call:
484 
485 template <int dim, int spacedim>
489  const Point<spacedim> &,
490  const Point<dim> &) const
491 {
492  // default implementation (should never be called)
493  Assert(false, ExcInternalError());
494  return {};
495 }
496 
497 
498 
499 template <>
500 Point<1>
503  const Point<1> & p,
504  const Point<1> & initial_p_unit) const
505 {
506  // dispatch to the various specializations for spacedim=dim,
507  // spacedim=dim+1, etc
508  return internal::MappingQImplementation::
509  do_transform_real_to_unit_cell_internal<1>(
510  p,
511  initial_p_unit,
512  this->compute_mapping_support_points(cell),
513  polynomials_1d,
514  renumber_lexicographic_to_hierarchic);
515 }
516 
517 
518 
519 template <>
520 Point<2>
523  const Point<2> & p,
524  const Point<2> & initial_p_unit) const
525 {
526  return internal::MappingQImplementation::
527  do_transform_real_to_unit_cell_internal<2>(
528  p,
529  initial_p_unit,
530  this->compute_mapping_support_points(cell),
531  polynomials_1d,
532  renumber_lexicographic_to_hierarchic);
533 }
534 
535 
536 
537 template <>
538 Point<3>
541  const Point<3> & p,
542  const Point<3> & initial_p_unit) const
543 {
544  return internal::MappingQImplementation::
545  do_transform_real_to_unit_cell_internal<3>(
546  p,
547  initial_p_unit,
548  this->compute_mapping_support_points(cell),
549  polynomials_1d,
550  renumber_lexicographic_to_hierarchic);
551 }
552 
553 
554 
555 template <>
556 Point<1>
559  const Point<2> & p,
560  const Point<1> & initial_p_unit) const
561 {
562  const int dim = 1;
563  const int spacedim = 2;
564 
565  const Quadrature<dim> point_quadrature(initial_p_unit);
566 
568  if (spacedim > dim)
569  update_flags |= update_jacobian_grads;
570  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
571  get_data(update_flags, point_quadrature));
572 
573  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
574 
575  // dispatch to the various specializations for spacedim=dim,
576  // spacedim=dim+1, etc
577  return internal::MappingQImplementation::
578  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
579  p,
580  initial_p_unit,
581  *mdata);
582 }
583 
584 
585 
586 template <>
587 Point<2>
590  const Point<3> & p,
591  const Point<2> & initial_p_unit) const
592 {
593  const int dim = 2;
594  const int spacedim = 3;
595 
596  const Quadrature<dim> point_quadrature(initial_p_unit);
597 
599  if (spacedim > dim)
600  update_flags |= update_jacobian_grads;
601  auto mdata = Utilities::dynamic_unique_cast<InternalData>(
602  get_data(update_flags, point_quadrature));
603 
604  mdata->mapping_support_points = this->compute_mapping_support_points(cell);
605 
606  // dispatch to the various specializations for spacedim=dim,
607  // spacedim=dim+1, etc
608  return internal::MappingQImplementation::
609  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
610  p,
611  initial_p_unit,
612  *mdata);
613 }
614 
615 template <>
616 Point<1>
619  const Point<3> &,
620  const Point<1> &) const
621 {
622  Assert(false, ExcNotImplemented());
623  return {};
624 }
625 
626 
627 
628 template <int dim, int spacedim>
631  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
632  const Point<spacedim> & p) const
633 {
634  // Use an exact formula if one is available. this is only the case
635  // for Q1 mappings in 1d, and in 2d if dim==spacedim
636  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
637  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
638  {
639  // The dimension-dependent algorithms are much faster (about 25-45x in
640  // 2d) but fail most of the time when the given point (p) is not in the
641  // cell. The dimension-independent Newton algorithm given below is
642  // slower, but more robust (though it still sometimes fails). Therefore
643  // this function implements the following strategy based on the
644  // p's dimension:
645  //
646  // * In 1d this mapping is linear, so the mapping is always invertible
647  // (and the exact formula is known) as long as the cell has non-zero
648  // length.
649  // * In 2d the exact (quadratic) formula is called first. If either the
650  // exact formula does not succeed (negative discriminant in the
651  // quadratic formula) or succeeds but finds a solution outside of the
652  // unit cell, then the Newton solver is called. The rationale for the
653  // second choice is that the exact formula may provide two different
654  // answers when mapping a point outside of the real cell, but the
655  // Newton solver (if it converges) will only return one answer.
656  // Otherwise the exact formula successfully found a point in the unit
657  // cell and that value is returned.
658  // * In 3d there is no (known to the authors) exact formula, so the Newton
659  // algorithm is used.
660  const auto vertices_ = this->get_vertices(cell);
661 
662  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
663  vertices;
664  for (unsigned int i = 0; i < vertices.size(); ++i)
665  vertices[i] = vertices_[i];
666 
667  try
668  {
669  switch (dim)
670  {
671  case 1:
672  {
673  // formula not subject to any issues in 1d
674  if (spacedim == 1)
676  vertices, p);
677  else
678  break;
679  }
680 
681  case 2:
682  {
683  const Point<dim> point =
685  p);
686 
687  // formula not guaranteed to work for points outside of
688  // the cell. only take the computed point if it lies
689  // inside the reference cell
690  const double eps = 1e-15;
691  if (-eps <= point(1) && point(1) <= 1 + eps &&
692  -eps <= point(0) && point(0) <= 1 + eps)
693  {
694  return point;
695  }
696  else
697  break;
698  }
699 
700  default:
701  {
702  // we should get here, based on the if-condition at the top
703  Assert(false, ExcInternalError());
704  }
705  }
706  }
707  catch (
709  {
710  // simply fall through and continue on to the standard Newton code
711  }
712  }
713  else
714  {
715  // we can't use an explicit formula,
716  }
717 
718 
719  // Find the initial value for the Newton iteration by a normal
720  // projection to the least square plane determined by the vertices
721  // of the cell
722  Point<dim> initial_p_unit;
723  if (this->preserves_vertex_locations())
724  {
725  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
726  // in 1d with spacedim > 1 the affine approximation is exact
727  if (dim == 1 && polynomial_degree == 1)
728  return initial_p_unit;
729  }
730  else
731  {
732  // else, we simply use the mid point
733  for (unsigned int d = 0; d < dim; ++d)
734  initial_p_unit[d] = 0.5;
735  }
736 
737  // perform the Newton iteration and return the result. note that this
738  // statement may throw an exception, which we simply pass up to the caller
739  const Point<dim> p_unit =
740  this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
741  if (p_unit[0] == std::numeric_limits<double>::infinity())
742  AssertThrow(false,
744  return p_unit;
745 }
746 
747 
748 
749 template <int dim, int spacedim>
750 void
752  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
753  const ArrayView<const Point<spacedim>> & real_points,
754  const ArrayView<Point<dim>> & unit_points) const
755 {
756  // Go to base class functions for dim < spacedim because it is not yet
757  // implemented with optimized code.
758  if (dim < spacedim)
759  {
761  real_points,
762  unit_points);
763  return;
764  }
765 
766  AssertDimension(real_points.size(), unit_points.size());
767  const std::vector<Point<spacedim>> support_points =
768  this->compute_mapping_support_points(cell);
769 
770  // From the given (high-order) support points, now only pick the first
771  // 2^dim points and construct an affine approximation from those.
773  inverse_approximation(support_points, unit_cell_support_points);
774 
775  const unsigned int n_points = real_points.size();
776  const unsigned int n_lanes = VectorizedArray<double>::size();
777 
778  // Use the more heavy VectorizedArray code path if there is more than
779  // one point left to compute
780  for (unsigned int i = 0; i < n_points; i += n_lanes)
781  if (n_points - i > 1)
782  {
784  for (unsigned int j = 0; j < n_lanes; ++j)
785  if (i + j < n_points)
786  for (unsigned int d = 0; d < spacedim; ++d)
787  p_vec[d][j] = real_points[i + j][d];
788  else
789  for (unsigned int d = 0; d < spacedim; ++d)
790  p_vec[d][j] = real_points[i][d];
791 
793  internal::MappingQImplementation::
794  do_transform_real_to_unit_cell_internal<dim, spacedim>(
795  p_vec,
796  inverse_approximation.compute(p_vec),
797  support_points,
798  polynomials_1d,
799  renumber_lexicographic_to_hierarchic);
800 
801  // If the vectorized computation failed, it could be that only some of
802  // the lanes failed but others would have succeeded if we had let them
803  // compute alone without interference (like negative Jacobian
804  // determinants) from other SIMD lanes. Repeat the computation in this
805  // unlikely case with scalar arguments.
806  for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
807  if (unit_point[0][j] == std::numeric_limits<double>::infinity())
808  unit_points[i + j] = internal::MappingQImplementation::
809  do_transform_real_to_unit_cell_internal<dim, spacedim>(
810  real_points[i + j],
811  inverse_approximation.compute(real_points[i + j]),
812  support_points,
813  polynomials_1d,
814  renumber_lexicographic_to_hierarchic);
815  else
816  for (unsigned int d = 0; d < dim; ++d)
817  unit_points[i + j][d] = unit_point[d][j];
818  }
819  else
820  unit_points[i] = internal::MappingQImplementation::
821  do_transform_real_to_unit_cell_internal<dim, spacedim>(
822  real_points[i],
823  inverse_approximation.compute(real_points[i]),
824  support_points,
825  polynomials_1d,
826  renumber_lexicographic_to_hierarchic);
827 }
828 
829 
830 
831 template <int dim, int spacedim>
834 {
835  // add flags if the respective quantities are necessary to compute
836  // what we need. note that some flags appear in both the conditions
837  // and in subsequent set operations. this leads to some circular
838  // logic. the only way to treat this is to iterate. since there are
839  // 5 if-clauses in the loop, it will take at most 5 iterations to
840  // converge. do them:
841  UpdateFlags out = in;
842  for (unsigned int i = 0; i < 5; ++i)
843  {
844  // The following is a little incorrect:
845  // If not applied on a face,
846  // update_boundary_forms does not
847  // make sense. On the other hand,
848  // it is necessary on a
849  // face. Currently,
850  // update_boundary_forms is simply
851  // ignored for the interior of a
852  // cell.
854  out |= update_boundary_forms;
855 
860 
861  if (out &
866 
867  // The contravariant transformation is used in the Piola
868  // transformation, which requires the determinant of the Jacobi
869  // matrix of the transformation. Because we have no way of
870  // knowing here whether the finite element wants to use the
871  // contravariant or the Piola transforms, we add the JxW values
872  // to the list of flags to be updated for each cell.
874  out |= update_volume_elements;
875 
876  // the same is true when computing normal vectors: they require
877  // the determinant of the Jacobian
878  if (out & update_normal_vectors)
879  out |= update_volume_elements;
880  }
881 
882  return out;
883 }
884 
885 
886 
887 template <int dim, int spacedim>
888 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
890  const Quadrature<dim> &q) const
891 {
892  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
893  std::make_unique<InternalData>(polynomial_degree);
894  auto &data = dynamic_cast<InternalData &>(*data_ptr);
895  data.initialize(this->requires_update_flags(update_flags), q, q.size());
896 
897  return data_ptr;
898 }
899 
900 
901 
902 template <int dim, int spacedim>
903 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
905  const UpdateFlags update_flags,
906  const hp::QCollection<dim - 1> &quadrature) const
907 {
908  AssertDimension(quadrature.size(), 1);
909 
910  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
911  std::make_unique<InternalData>(polynomial_degree);
912  auto &data = dynamic_cast<InternalData &>(*data_ptr);
913  data.initialize_face(this->requires_update_flags(update_flags),
915  ReferenceCells::get_hypercube<dim>(), quadrature[0]),
916  quadrature[0].size());
917 
918  return data_ptr;
919 }
920 
921 
922 
923 template <int dim, int spacedim>
924 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
926  const UpdateFlags update_flags,
927  const Quadrature<dim - 1> &quadrature) const
928 {
929  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
930  std::make_unique<InternalData>(polynomial_degree);
931  auto &data = dynamic_cast<InternalData &>(*data_ptr);
932  data.initialize_face(this->requires_update_flags(update_flags),
934  ReferenceCells::get_hypercube<dim>(), quadrature),
935  quadrature.size());
936 
937  return data_ptr;
938 }
939 
940 
941 
942 template <int dim, int spacedim>
945  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
946  const CellSimilarity::Similarity cell_similarity,
947  const Quadrature<dim> & quadrature,
948  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
950  &output_data) const
951 {
952  // ensure that the following static_cast is really correct:
953  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
954  ExcInternalError());
955  const InternalData &data = static_cast<const InternalData &>(internal_data);
956 
957  const unsigned int n_q_points = quadrature.size();
958 
959  // recompute the support points of the transformation of this
960  // cell. we tried to be clever here in an earlier version of the
961  // library by checking whether the cell is the same as the one we
962  // had visited last, but it turns out to be difficult to determine
963  // that because a cell for the purposes of a mapping is
964  // characterized not just by its (triangulation, level, index)
965  // triple, but also by the locations of its vertices, the manifold
966  // object attached to the cell and all of its bounding faces/edges,
967  // etc. to reliably test that the "cell" we are on is, therefore,
968  // not easily done
969  data.mapping_support_points = this->compute_mapping_support_points(cell);
970  data.cell_of_current_support_points = cell;
971 
972  // if the order of the mapping is greater than 1, then do not reuse any cell
973  // similarity information. This is necessary because the cell similarity
974  // value is computed with just cell vertices and does not take into account
975  // cell curvature.
976  const CellSimilarity::Similarity computed_cell_similarity =
977  (polynomial_degree == 1 && this->preserves_vertex_locations() ?
978  cell_similarity :
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:704
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:123
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:352
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:63
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:214
InternalData(const unsigned int polynomial_degree)
Definition: mapping_q.cc:51
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
Definition: mapping_q.cc:270
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:84
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:889
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:925
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:944
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:904
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:487
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:452
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:434
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:751
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:630
MappingQ(const unsigned int polynomial_degree)
Definition: mapping_q.cc:361
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_q.cc:833
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:443
Abstract base class for mapping classes.
Definition: mapping.h:314
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:110
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 DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const override
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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:75
@ mapping_covariant_gradient
Definition: mapping.h:96
@ mapping_contravariant
Definition: mapping.h:90
@ mapping_contravariant_hessian
Definition: mapping.h:152
@ mapping_covariant_hessian
Definition: mapping.h:146
@ mapping_contravariant_gradient
Definition: mapping.h:102
@ mapping_piola_gradient
Definition: mapping.h:116
@ mapping_piola_hessian
Definition: mapping.h:158
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: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:966
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:213
const types::manifold_id flat_manifold_id
Definition: types.h:286
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)