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