Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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_fe.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
23 #include <deal.II/base/table.h>
25 
26 #include <deal.II/fe/fe_poly.h>
27 #include <deal.II/fe/fe_values.h>
28 #include <deal.II/fe/mapping_fe.h>
29 
31 #include <deal.II/grid/tria.h>
33 
35 
36 #include <boost/container/small_vector.hpp>
37 
38 #include <algorithm>
39 #include <array>
40 #include <cmath>
41 #include <memory>
42 #include <numeric>
43 
44 
46 
47 
48 template <int dim, int spacedim>
51  : fe(fe)
52  , polynomial_degree(fe.tensor_degree())
53  , n_shape_functions(fe.n_dofs_per_cell())
54 {}
55 
56 
57 
58 template <int dim, int spacedim>
59 std::size_t
61 {
62  return (
65  MemoryConsumption::memory_consumption(shape_derivatives) +
68  MemoryConsumption::memory_consumption(unit_tangentials) +
70  MemoryConsumption::memory_consumption(mapping_support_points) +
71  MemoryConsumption::memory_consumption(cell_of_current_support_points) +
72  MemoryConsumption::memory_consumption(volume_elements) +
74  MemoryConsumption::memory_consumption(n_shape_functions));
75 }
76 
77 
78 template <int dim, int spacedim>
79 void
81  const UpdateFlags update_flags,
82  const Quadrature<dim> &q,
83  const unsigned int n_original_q_points)
84 {
85  // store the flags in the internal data object so we can access them
86  // in fill_fe_*_values()
87  this->update_each = update_flags;
88 
89  const unsigned int n_q_points = q.size();
90 
91  if (this->update_each & update_covariant_transformation)
92  covariant.resize(n_original_q_points);
93 
94  if (this->update_each & update_contravariant_transformation)
95  contravariant.resize(n_original_q_points);
96 
97  if (this->update_each & update_volume_elements)
98  volume_elements.resize(n_original_q_points);
99 
100  // see if we need the (transformation) shape function values
101  // and/or gradients and resize the necessary arrays
102  if (this->update_each & update_quadrature_points)
103  shape_values.resize(n_shape_functions * n_q_points);
104 
105  if (this->update_each &
113  shape_derivatives.resize(n_shape_functions * n_q_points);
114 
115  if (this->update_each &
117  shape_second_derivatives.resize(n_shape_functions * n_q_points);
118 
119  if (this->update_each & (update_jacobian_2nd_derivatives |
121  shape_third_derivatives.resize(n_shape_functions * n_q_points);
122 
123  if (this->update_each & (update_jacobian_3rd_derivatives |
125  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
126 
127  // now also fill the various fields with their correct values
128  compute_shape_function_values(q.get_points());
129 
130  // copy (projected) quadrature weights
131  quadrature_weights = q.get_weights();
132 }
133 
134 
135 
136 template <int dim, int spacedim>
137 void
139  const UpdateFlags update_flags,
140  const Quadrature<dim> &q,
141  const unsigned int n_original_q_points)
142 {
143  initialize(update_flags, q, n_original_q_points);
144 
145  if (this->update_each &
148  {
149  aux.resize(dim - 1,
150  std::vector<Tensor<1, spacedim>>(n_original_q_points));
151 
152  // Compute tangentials to the unit cell.
153  const auto reference_cell = this->fe.reference_cell();
154  const auto n_faces = reference_cell.n_faces();
155 
156  for (unsigned int i = 0; i < n_faces; ++i)
157  {
158  unit_tangentials[i].resize(n_original_q_points);
159  std::fill(unit_tangentials[i].begin(),
160  unit_tangentials[i].end(),
161  reference_cell.template unit_tangential_vectors<dim>(i, 0));
162  if (dim > 2)
163  {
164  unit_tangentials[n_faces + i].resize(n_original_q_points);
165  std::fill(
166  unit_tangentials[n_faces + i].begin(),
167  unit_tangentials[n_faces + i].end(),
168  reference_cell.template unit_tangential_vectors<dim>(i, 1));
169  }
170  }
171  }
172 }
173 
174 
175 
176 template <int dim, int spacedim>
177 void
179  const std::vector<Point<dim>> &unit_points)
180 {
181  const auto fe_poly = dynamic_cast<const FE_Poly<dim, spacedim> *>(&this->fe);
182 
183  Assert(fe_poly != nullptr, ExcNotImplemented());
184 
185  const auto &tensor_pols = fe_poly->get_poly_space();
186 
187  const unsigned int n_shape_functions = fe.n_dofs_per_cell();
188  const unsigned int n_points = unit_points.size();
189 
190  std::vector<double> values;
191  std::vector<Tensor<1, dim>> grads;
192  if (shape_values.size() != 0)
193  {
194  Assert(shape_values.size() == n_shape_functions * n_points,
195  ExcInternalError());
196  values.resize(n_shape_functions);
197  }
198  if (shape_derivatives.size() != 0)
199  {
200  Assert(shape_derivatives.size() == n_shape_functions * n_points,
201  ExcInternalError());
202  grads.resize(n_shape_functions);
203  }
204 
205  std::vector<Tensor<2, dim>> grad2;
206  if (shape_second_derivatives.size() != 0)
207  {
208  Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
209  ExcInternalError());
210  grad2.resize(n_shape_functions);
211  }
212 
213  std::vector<Tensor<3, dim>> grad3;
214  if (shape_third_derivatives.size() != 0)
215  {
216  Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
217  ExcInternalError());
218  grad3.resize(n_shape_functions);
219  }
220 
221  std::vector<Tensor<4, dim>> grad4;
222  if (shape_fourth_derivatives.size() != 0)
223  {
224  Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
225  ExcInternalError());
226  grad4.resize(n_shape_functions);
227  }
228 
229 
230  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
231  shape_second_derivatives.size() != 0 ||
232  shape_third_derivatives.size() != 0 ||
233  shape_fourth_derivatives.size() != 0)
234  for (unsigned int point = 0; point < n_points; ++point)
235  {
236  tensor_pols.evaluate(
237  unit_points[point], values, grads, grad2, grad3, grad4);
238 
239  if (shape_values.size() != 0)
240  for (unsigned int i = 0; i < n_shape_functions; ++i)
241  shape(point, i) = values[i];
242 
243  if (shape_derivatives.size() != 0)
244  for (unsigned int i = 0; i < n_shape_functions; ++i)
245  derivative(point, i) = grads[i];
246 
247  if (shape_second_derivatives.size() != 0)
248  for (unsigned int i = 0; i < n_shape_functions; ++i)
249  second_derivative(point, i) = grad2[i];
250 
251  if (shape_third_derivatives.size() != 0)
252  for (unsigned int i = 0; i < n_shape_functions; ++i)
253  third_derivative(point, i) = grad3[i];
254 
255  if (shape_fourth_derivatives.size() != 0)
256  for (unsigned int i = 0; i < n_shape_functions; ++i)
257  fourth_derivative(point, i) = grad4[i];
258  }
259 }
260 
261 
262 namespace internal
263 {
264  namespace MappingFEImplementation
265  {
266  namespace
267  {
274  template <int dim, int spacedim>
275  void
276  maybe_compute_q_points(
277  const typename QProjector<dim>::DataSetDescriptor data_set,
278  const typename ::MappingFE<dim, spacedim>::InternalData &data,
279  std::vector<Point<spacedim>> &quadrature_points,
280  const unsigned int n_q_points)
281  {
282  const UpdateFlags update_flags = data.update_each;
283 
284  if (update_flags & update_quadrature_points)
285  for (unsigned int point = 0; point < n_q_points; ++point)
286  {
287  const double *shape = &data.shape(point + data_set, 0);
288  Point<spacedim> result =
289  (shape[0] * data.mapping_support_points[0]);
290  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
291  for (unsigned int i = 0; i < spacedim; ++i)
292  result[i] += shape[k] * data.mapping_support_points[k][i];
293  quadrature_points[point] = result;
294  }
295  }
296 
297 
298 
307  template <int dim, int spacedim>
308  void
309  maybe_update_Jacobians(
310  const CellSimilarity::Similarity cell_similarity,
311  const typename ::QProjector<dim>::DataSetDescriptor data_set,
312  const typename ::MappingFE<dim, spacedim>::InternalData &data,
313  const unsigned int n_q_points)
314  {
315  const UpdateFlags update_flags = data.update_each;
316 
317  if (update_flags & update_contravariant_transformation)
318  // if the current cell is just a
319  // translation of the previous one, no
320  // need to recompute jacobians...
321  if (cell_similarity != CellSimilarity::translation)
322  {
323  std::fill(data.contravariant.begin(),
324  data.contravariant.end(),
326 
327  Assert(data.n_shape_functions > 0, ExcInternalError());
328 
329  for (unsigned int point = 0; point < n_q_points; ++point)
330  {
331  double result[spacedim][dim];
332 
333  // peel away part of sum to avoid zeroing the
334  // entries and adding for the first time
335  for (unsigned int i = 0; i < spacedim; ++i)
336  for (unsigned int j = 0; j < dim; ++j)
337  result[i][j] = data.derivative(point + data_set, 0)[j] *
338  data.mapping_support_points[0][i];
339  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
340  for (unsigned int i = 0; i < spacedim; ++i)
341  for (unsigned int j = 0; j < dim; ++j)
342  result[i][j] +=
343  data.derivative(point + data_set, k)[j] *
344  data.mapping_support_points[k][i];
345 
346  // write result into contravariant data. for
347  // j=dim in the case dim<spacedim, there will
348  // never be any nonzero data that arrives in
349  // here, so it is ok anyway because it was
350  // initialized to zero at the initialization
351  for (unsigned int i = 0; i < spacedim; ++i)
352  for (unsigned int j = 0; j < dim; ++j)
353  data.contravariant[point][i][j] = result[i][j];
354  }
355  }
356 
357  if (update_flags & update_covariant_transformation)
358  if (cell_similarity != CellSimilarity::translation)
359  {
360  for (unsigned int point = 0; point < n_q_points; ++point)
361  {
362  data.covariant[point] =
363  (data.contravariant[point]).covariant_form();
364  }
365  }
366 
367  if (update_flags & update_volume_elements)
368  if (cell_similarity != CellSimilarity::translation)
369  {
370  for (unsigned int point = 0; point < n_q_points; ++point)
371  data.volume_elements[point] =
372  data.contravariant[point].determinant();
373  }
374  }
375 
382  template <int dim, int spacedim>
383  void
385  const CellSimilarity::Similarity cell_similarity,
386  const typename QProjector<dim>::DataSetDescriptor data_set,
387  const typename ::MappingFE<dim, spacedim>::InternalData &data,
388  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads,
389  const unsigned int n_q_points)
390  {
391  const UpdateFlags update_flags = data.update_each;
392  if (update_flags & update_jacobian_grads)
393  {
394  AssertIndexRange(n_q_points, jacobian_grads.size() + 1);
395 
396  if (cell_similarity != CellSimilarity::translation)
397  for (unsigned int point = 0; point < n_q_points; ++point)
398  {
399  const Tensor<2, dim> *second =
400  &data.second_derivative(point + data_set, 0);
401  double result[spacedim][dim][dim];
402  for (unsigned int i = 0; i < spacedim; ++i)
403  for (unsigned int j = 0; j < dim; ++j)
404  for (unsigned int l = 0; l < dim; ++l)
405  result[i][j][l] =
406  (second[0][j][l] * data.mapping_support_points[0][i]);
407  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
408  for (unsigned int i = 0; i < spacedim; ++i)
409  for (unsigned int j = 0; j < dim; ++j)
410  for (unsigned int l = 0; l < dim; ++l)
411  result[i][j][l] +=
412  (second[k][j][l] *
413  data.mapping_support_points[k][i]);
414 
415  for (unsigned int i = 0; i < spacedim; ++i)
416  for (unsigned int j = 0; j < dim; ++j)
417  for (unsigned int l = 0; l < dim; ++l)
418  jacobian_grads[point][i][j][l] = result[i][j][l];
419  }
420  }
421  }
422 
429  template <int dim, int spacedim>
430  void
432  const CellSimilarity::Similarity cell_similarity,
433  const typename QProjector<dim>::DataSetDescriptor data_set,
434  const typename ::MappingFE<dim, spacedim>::InternalData &data,
435  std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads,
436  const unsigned int n_q_points)
437  {
438  const UpdateFlags update_flags = data.update_each;
439  if (update_flags & update_jacobian_pushed_forward_grads)
440  {
441  AssertIndexRange(n_q_points,
442  jacobian_pushed_forward_grads.size() + 1);
443 
444  if (cell_similarity != CellSimilarity::translation)
445  {
446  double tmp[spacedim][spacedim][spacedim];
447  for (unsigned int point = 0; point < n_q_points; ++point)
448  {
449  const Tensor<2, dim> *second =
450  &data.second_derivative(point + data_set, 0);
451  double result[spacedim][dim][dim];
452  for (unsigned int i = 0; i < spacedim; ++i)
453  for (unsigned int j = 0; j < dim; ++j)
454  for (unsigned int l = 0; l < dim; ++l)
455  result[i][j][l] = (second[0][j][l] *
456  data.mapping_support_points[0][i]);
457  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
458  for (unsigned int i = 0; i < spacedim; ++i)
459  for (unsigned int j = 0; j < dim; ++j)
460  for (unsigned int l = 0; l < dim; ++l)
461  result[i][j][l] +=
462  (second[k][j][l] *
463  data.mapping_support_points[k][i]);
464 
465  // first push forward the j-components
466  for (unsigned int i = 0; i < spacedim; ++i)
467  for (unsigned int j = 0; j < spacedim; ++j)
468  for (unsigned int l = 0; l < dim; ++l)
469  {
470  tmp[i][j][l] =
471  result[i][0][l] * data.covariant[point][j][0];
472  for (unsigned int jr = 1; jr < dim; ++jr)
473  {
474  tmp[i][j][l] += result[i][jr][l] *
475  data.covariant[point][j][jr];
476  }
477  }
478 
479  // now, pushing forward the l-components
480  for (unsigned int i = 0; i < spacedim; ++i)
481  for (unsigned int j = 0; j < spacedim; ++j)
482  for (unsigned int l = 0; l < spacedim; ++l)
483  {
484  jacobian_pushed_forward_grads[point][i][j][l] =
485  tmp[i][j][0] * data.covariant[point][l][0];
486  for (unsigned int lr = 1; lr < dim; ++lr)
487  {
488  jacobian_pushed_forward_grads[point][i][j][l] +=
489  tmp[i][j][lr] * data.covariant[point][l][lr];
490  }
491  }
492  }
493  }
494  }
495  }
496 
503  template <int dim, int spacedim>
504  void
506  const CellSimilarity::Similarity cell_similarity,
507  const typename QProjector<dim>::DataSetDescriptor data_set,
508  const typename ::MappingFE<dim, spacedim>::InternalData &data,
509  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives,
510  const unsigned int n_q_points)
511  {
512  const UpdateFlags update_flags = data.update_each;
513  if (update_flags & update_jacobian_2nd_derivatives)
514  {
515  AssertIndexRange(n_q_points, jacobian_2nd_derivatives.size() + 1);
516 
517  if (cell_similarity != CellSimilarity::translation)
518  {
519  for (unsigned int point = 0; point < n_q_points; ++point)
520  {
521  const Tensor<3, dim> *third =
522  &data.third_derivative(point + data_set, 0);
523  double result[spacedim][dim][dim][dim];
524  for (unsigned int i = 0; i < spacedim; ++i)
525  for (unsigned int j = 0; j < dim; ++j)
526  for (unsigned int l = 0; l < dim; ++l)
527  for (unsigned int m = 0; m < dim; ++m)
528  result[i][j][l][m] =
529  (third[0][j][l][m] *
530  data.mapping_support_points[0][i]);
531  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
532  for (unsigned int i = 0; i < spacedim; ++i)
533  for (unsigned int j = 0; j < dim; ++j)
534  for (unsigned int l = 0; l < dim; ++l)
535  for (unsigned int m = 0; m < dim; ++m)
536  result[i][j][l][m] +=
537  (third[k][j][l][m] *
538  data.mapping_support_points[k][i]);
539 
540  for (unsigned int i = 0; i < spacedim; ++i)
541  for (unsigned int j = 0; j < dim; ++j)
542  for (unsigned int l = 0; l < dim; ++l)
543  for (unsigned int m = 0; m < dim; ++m)
544  jacobian_2nd_derivatives[point][i][j][l][m] =
545  result[i][j][l][m];
546  }
547  }
548  }
549  }
550 
558  template <int dim, int spacedim>
559  void
561  const CellSimilarity::Similarity cell_similarity,
562  const typename QProjector<dim>::DataSetDescriptor data_set,
563  const typename ::MappingFE<dim, spacedim>::InternalData &data,
564  std::vector<Tensor<4, spacedim>>
565  &jacobian_pushed_forward_2nd_derivatives,
566  const unsigned int n_q_points)
567  {
568  const UpdateFlags update_flags = data.update_each;
570  {
571  AssertIndexRange(n_q_points,
572  jacobian_pushed_forward_2nd_derivatives.size() +
573  1);
574 
575  if (cell_similarity != CellSimilarity::translation)
576  {
577  double tmp[spacedim][spacedim][spacedim][spacedim];
578  for (unsigned int point = 0; point < n_q_points; ++point)
579  {
580  const Tensor<3, dim> *third =
581  &data.third_derivative(point + data_set, 0);
582  double result[spacedim][dim][dim][dim];
583  for (unsigned int i = 0; i < spacedim; ++i)
584  for (unsigned int j = 0; j < dim; ++j)
585  for (unsigned int l = 0; l < dim; ++l)
586  for (unsigned int m = 0; m < dim; ++m)
587  result[i][j][l][m] =
588  (third[0][j][l][m] *
589  data.mapping_support_points[0][i]);
590  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
591  for (unsigned int i = 0; i < spacedim; ++i)
592  for (unsigned int j = 0; j < dim; ++j)
593  for (unsigned int l = 0; l < dim; ++l)
594  for (unsigned int m = 0; m < dim; ++m)
595  result[i][j][l][m] +=
596  (third[k][j][l][m] *
597  data.mapping_support_points[k][i]);
598 
599  // push forward the j-coordinate
600  for (unsigned int i = 0; i < spacedim; ++i)
601  for (unsigned int j = 0; j < spacedim; ++j)
602  for (unsigned int l = 0; l < dim; ++l)
603  for (unsigned int m = 0; m < dim; ++m)
604  {
605  jacobian_pushed_forward_2nd_derivatives
606  [point][i][j][l][m] =
607  result[i][0][l][m] *
608  data.covariant[point][j][0];
609  for (unsigned int jr = 1; jr < dim; ++jr)
610  jacobian_pushed_forward_2nd_derivatives[point]
611  [i][j][l]
612  [m] +=
613  result[i][jr][l][m] *
614  data.covariant[point][j][jr];
615  }
616 
617  // push forward the l-coordinate
618  for (unsigned int i = 0; i < spacedim; ++i)
619  for (unsigned int j = 0; j < spacedim; ++j)
620  for (unsigned int l = 0; l < spacedim; ++l)
621  for (unsigned int m = 0; m < dim; ++m)
622  {
623  tmp[i][j][l][m] =
624  jacobian_pushed_forward_2nd_derivatives[point]
625  [i][j][0]
626  [m] *
627  data.covariant[point][l][0];
628  for (unsigned int lr = 1; lr < dim; ++lr)
629  tmp[i][j][l][m] +=
630  jacobian_pushed_forward_2nd_derivatives
631  [point][i][j][lr][m] *
632  data.covariant[point][l][lr];
633  }
634 
635  // push forward the m-coordinate
636  for (unsigned int i = 0; i < spacedim; ++i)
637  for (unsigned int j = 0; j < spacedim; ++j)
638  for (unsigned int l = 0; l < spacedim; ++l)
639  for (unsigned int m = 0; m < spacedim; ++m)
640  {
641  jacobian_pushed_forward_2nd_derivatives
642  [point][i][j][l][m] =
643  tmp[i][j][l][0] * data.covariant[point][m][0];
644  for (unsigned int mr = 1; mr < dim; ++mr)
645  jacobian_pushed_forward_2nd_derivatives[point]
646  [i][j][l]
647  [m] +=
648  tmp[i][j][l][mr] *
649  data.covariant[point][m][mr];
650  }
651  }
652  }
653  }
654  }
655 
662  template <int dim, int spacedim>
663  void
665  const CellSimilarity::Similarity cell_similarity,
666  const typename QProjector<dim>::DataSetDescriptor data_set,
667  const typename ::MappingFE<dim, spacedim>::InternalData &data,
668  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives,
669  const unsigned int n_q_points)
670  {
671  const UpdateFlags update_flags = data.update_each;
672  if (update_flags & update_jacobian_3rd_derivatives)
673  {
674  AssertIndexRange(n_q_points, jacobian_3rd_derivatives.size() + 1);
675 
676  if (cell_similarity != CellSimilarity::translation)
677  {
678  for (unsigned int point = 0; point < n_q_points; ++point)
679  {
680  const Tensor<4, dim> *fourth =
681  &data.fourth_derivative(point + data_set, 0);
682  double result[spacedim][dim][dim][dim][dim];
683  for (unsigned int i = 0; i < spacedim; ++i)
684  for (unsigned int j = 0; j < dim; ++j)
685  for (unsigned int l = 0; l < dim; ++l)
686  for (unsigned int m = 0; m < dim; ++m)
687  for (unsigned int n = 0; n < dim; ++n)
688  result[i][j][l][m][n] =
689  (fourth[0][j][l][m][n] *
690  data.mapping_support_points[0][i]);
691  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
692  for (unsigned int i = 0; i < spacedim; ++i)
693  for (unsigned int j = 0; j < dim; ++j)
694  for (unsigned int l = 0; l < dim; ++l)
695  for (unsigned int m = 0; m < dim; ++m)
696  for (unsigned int n = 0; n < dim; ++n)
697  result[i][j][l][m][n] +=
698  (fourth[k][j][l][m][n] *
699  data.mapping_support_points[k][i]);
700 
701  for (unsigned int i = 0; i < spacedim; ++i)
702  for (unsigned int j = 0; j < dim; ++j)
703  for (unsigned int l = 0; l < dim; ++l)
704  for (unsigned int m = 0; m < dim; ++m)
705  for (unsigned int n = 0; n < dim; ++n)
706  jacobian_3rd_derivatives[point][i][j][l][m][n] =
707  result[i][j][l][m][n];
708  }
709  }
710  }
711  }
712 
720  template <int dim, int spacedim>
721  void
723  const CellSimilarity::Similarity cell_similarity,
724  const typename QProjector<dim>::DataSetDescriptor data_set,
725  const typename ::MappingFE<dim, spacedim>::InternalData &data,
726  std::vector<Tensor<5, spacedim>>
727  &jacobian_pushed_forward_3rd_derivatives,
728  const unsigned int n_q_points)
729  {
730  const UpdateFlags update_flags = data.update_each;
732  {
733  AssertIndexRange(n_q_points,
734  jacobian_pushed_forward_3rd_derivatives.size() +
735  1);
736 
737  if (cell_similarity != CellSimilarity::translation)
738  {
739  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
740  for (unsigned int point = 0; point < n_q_points; ++point)
741  {
742  const Tensor<4, dim> *fourth =
743  &data.fourth_derivative(point + data_set, 0);
744  double result[spacedim][dim][dim][dim][dim];
745  for (unsigned int i = 0; i < spacedim; ++i)
746  for (unsigned int j = 0; j < dim; ++j)
747  for (unsigned int l = 0; l < dim; ++l)
748  for (unsigned int m = 0; m < dim; ++m)
749  for (unsigned int n = 0; n < dim; ++n)
750  result[i][j][l][m][n] =
751  (fourth[0][j][l][m][n] *
752  data.mapping_support_points[0][i]);
753  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
754  for (unsigned int i = 0; i < spacedim; ++i)
755  for (unsigned int j = 0; j < dim; ++j)
756  for (unsigned int l = 0; l < dim; ++l)
757  for (unsigned int m = 0; m < dim; ++m)
758  for (unsigned int n = 0; n < dim; ++n)
759  result[i][j][l][m][n] +=
760  (fourth[k][j][l][m][n] *
761  data.mapping_support_points[k][i]);
762 
763  // push-forward the j-coordinate
764  for (unsigned int i = 0; i < spacedim; ++i)
765  for (unsigned int j = 0; j < spacedim; ++j)
766  for (unsigned int l = 0; l < dim; ++l)
767  for (unsigned int m = 0; m < dim; ++m)
768  for (unsigned int n = 0; n < dim; ++n)
769  {
770  tmp[i][j][l][m][n] =
771  result[i][0][l][m][n] *
772  data.covariant[point][j][0];
773  for (unsigned int jr = 1; jr < dim; ++jr)
774  tmp[i][j][l][m][n] +=
775  result[i][jr][l][m][n] *
776  data.covariant[point][j][jr];
777  }
778 
779  // push-forward the l-coordinate
780  for (unsigned int i = 0; i < spacedim; ++i)
781  for (unsigned int j = 0; j < spacedim; ++j)
782  for (unsigned int l = 0; l < spacedim; ++l)
783  for (unsigned int m = 0; m < dim; ++m)
784  for (unsigned int n = 0; n < dim; ++n)
785  {
786  jacobian_pushed_forward_3rd_derivatives
787  [point][i][j][l][m][n] =
788  tmp[i][j][0][m][n] *
789  data.covariant[point][l][0];
790  for (unsigned int lr = 1; lr < dim; ++lr)
791  jacobian_pushed_forward_3rd_derivatives
792  [point][i][j][l][m][n] +=
793  tmp[i][j][lr][m][n] *
794  data.covariant[point][l][lr];
795  }
796 
797  // push-forward the m-coordinate
798  for (unsigned int i = 0; i < spacedim; ++i)
799  for (unsigned int j = 0; j < spacedim; ++j)
800  for (unsigned int l = 0; l < spacedim; ++l)
801  for (unsigned int m = 0; m < spacedim; ++m)
802  for (unsigned int n = 0; n < dim; ++n)
803  {
804  tmp[i][j][l][m][n] =
805  jacobian_pushed_forward_3rd_derivatives
806  [point][i][j][l][0][n] *
807  data.covariant[point][m][0];
808  for (unsigned int mr = 1; mr < dim; ++mr)
809  tmp[i][j][l][m][n] +=
810  jacobian_pushed_forward_3rd_derivatives
811  [point][i][j][l][mr][n] *
812  data.covariant[point][m][mr];
813  }
814 
815  // push-forward the n-coordinate
816  for (unsigned int i = 0; i < spacedim; ++i)
817  for (unsigned int j = 0; j < spacedim; ++j)
818  for (unsigned int l = 0; l < spacedim; ++l)
819  for (unsigned int m = 0; m < spacedim; ++m)
820  for (unsigned int n = 0; n < spacedim; ++n)
821  {
822  jacobian_pushed_forward_3rd_derivatives
823  [point][i][j][l][m][n] =
824  tmp[i][j][l][m][0] *
825  data.covariant[point][n][0];
826  for (unsigned int nr = 1; nr < dim; ++nr)
827  jacobian_pushed_forward_3rd_derivatives
828  [point][i][j][l][m][n] +=
829  tmp[i][j][l][m][nr] *
830  data.covariant[point][n][nr];
831  }
832  }
833  }
834  }
835  }
836  } // namespace
837  } // namespace MappingFEImplementation
838 } // namespace internal
839 
840 
841 
842 template <int dim, int spacedim>
844  : fe(fe.clone())
845  , polynomial_degree(fe.tensor_degree())
846 {
848  ExcMessage("It only makes sense to create polynomial mappings "
849  "with a polynomial degree greater or equal to one."));
850  Assert(fe.n_components() == 1, ExcNotImplemented());
851 
852  Assert(fe.has_support_points(), ExcNotImplemented());
853 
854  const auto &mapping_support_points = fe.get_unit_support_points();
855 
856  const auto reference_cell = fe.reference_cell();
857 
858  const unsigned int n_points = mapping_support_points.size();
859  const unsigned int n_shape_functions = reference_cell.n_vertices();
860 
862  Table<2, double>(n_points, n_shape_functions);
863 
864  for (unsigned int point = 0; point < n_points; ++point)
865  for (unsigned int i = 0; i < n_shape_functions; ++i)
867  reference_cell.d_linear_shape_function(mapping_support_points[point],
868  i);
869 }
870 
871 
872 
873 template <int dim, int spacedim>
875  : fe(mapping.fe->clone())
876  , polynomial_degree(mapping.polynomial_degree)
877  , mapping_support_point_weights(mapping.mapping_support_point_weights)
878 {}
879 
880 
881 
882 template <int dim, int spacedim>
883 std::unique_ptr<Mapping<dim, spacedim>>
885 {
886  return std::make_unique<MappingFE<dim, spacedim>>(*this);
887 }
888 
889 
890 
891 template <int dim, int spacedim>
892 unsigned int
894 {
895  return polynomial_degree;
896 }
897 
898 
899 
900 template <int dim, int spacedim>
903  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
904  const Point<dim> &p) const
905 {
906  const std::vector<Point<spacedim>> support_points =
907  this->compute_mapping_support_points(cell);
908 
909  Point<spacedim> mapped_point;
910 
911  for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
912  mapped_point += support_points[i] * this->fe->shape_value(i, p);
913 
914  return mapped_point;
915 }
916 
917 
918 
919 template <int dim, int spacedim>
922  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
923  const Point<spacedim> &p) const
924 {
925  const std::vector<Point<spacedim>> support_points =
926  this->compute_mapping_support_points(cell);
927 
928  const double eps = 1.e-12 * cell->diameter();
929  const unsigned int loop_limit = 10;
930 
931  Point<dim> p_unit;
932 
933  unsigned int loop = 0;
934 
935  // This loop solves the following problem:
936  // grad_F^T residual + (grad_F^T grad_F + grad_F^T hess_F^T dp) dp = 0
937  // where the term
938  // (grad_F^T hess_F dp) is approximated by (-hess_F * residual)
939  // This is basically a second order approximation of Newton method, where the
940  // Jacobian is corrected with a higher order term coming from the hessian.
941  do
942  {
943  Point<spacedim> mapped_point;
944 
945  // Transpose of the gradient map
948 
949  for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
950  {
951  mapped_point += support_points[i] * this->fe->shape_value(i, p_unit);
952  const auto grad_F_i = this->fe->shape_grad(i, p_unit);
953  const auto hessian_F_i = this->fe->shape_grad_grad(i, p_unit);
954  for (unsigned int j = 0; j < dim; ++j)
955  {
956  grad_FT[j] += grad_F_i[j] * support_points[i];
957  for (unsigned int l = 0; l < dim; ++l)
958  hess_FT[j][l] += hessian_F_i[j][l] * support_points[i];
959  }
960  }
961 
962  // Residual
963  const auto residual = p - mapped_point;
964  // Project the residual on the reference coordinate system
965  // to compute the error, and to filter components orthogonal to the
966  // manifold, and compute a 2nd order correction of the metric tensor
967  const auto grad_FT_residual = apply_transformation(grad_FT, residual);
968 
969  // Do not invert nor compute the metric if not necessary.
970  if (grad_FT_residual.norm() <= eps)
971  break;
972 
973  // Now compute the (corrected) metric tensor
974  Tensor<2, dim> corrected_metric_tensor;
975  for (unsigned int j = 0; j < dim; ++j)
976  for (unsigned int l = 0; l < dim; ++l)
977  corrected_metric_tensor[j][l] =
978  -grad_FT[j] * grad_FT[l] + hess_FT[j][l] * residual;
979 
980  // And compute the update
981  const auto g_inverse = invert(corrected_metric_tensor);
982  p_unit -= Point<dim>(g_inverse * grad_FT_residual);
983 
984  ++loop;
985  }
986  while (loop < loop_limit);
987 
988  // Here we check that in the last execution of while the first
989  // condition was already wrong, meaning the residual was below
990  // eps. Only if the first condition failed, loop will have been
991  // increased and tested, and thus have reached the limit.
992  AssertThrow(loop < loop_limit,
994 
995  return p_unit;
996 }
997 
998 
999 
1000 template <int dim, int spacedim>
1003 {
1004  // add flags if the respective quantities are necessary to compute
1005  // what we need. note that some flags appear in both the conditions
1006  // and in subsequent set operations. this leads to some circular
1007  // logic. the only way to treat this is to iterate. since there are
1008  // 5 if-clauses in the loop, it will take at most 5 iterations to
1009  // converge. do them:
1010  UpdateFlags out = in;
1011  for (unsigned int i = 0; i < 5; ++i)
1012  {
1013  // The following is a little incorrect:
1014  // If not applied on a face,
1015  // update_boundary_forms does not
1016  // make sense. On the other hand,
1017  // it is necessary on a
1018  // face. Currently,
1019  // update_boundary_forms is simply
1020  // ignored for the interior of a
1021  // cell.
1023  out |= update_boundary_forms;
1024 
1029 
1030  if (out &
1035 
1036  // The contravariant transformation is used in the Piola
1037  // transformation, which requires the determinant of the Jacobi
1038  // matrix of the transformation. Because we have no way of
1039  // knowing here whether the finite element wants to use the
1040  // contravariant or the Piola transforms, we add the JxW values
1041  // to the list of flags to be updated for each cell.
1043  out |= update_volume_elements;
1044 
1045  // the same is true when computing normal vectors: they require
1046  // the determinant of the Jacobian
1047  if (out & update_normal_vectors)
1048  out |= update_volume_elements;
1049  }
1050 
1051  return out;
1052 }
1053 
1054 
1055 
1056 template <int dim, int spacedim>
1057 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1059  const Quadrature<dim> &q) const
1060 {
1061  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1062  std::make_unique<InternalData>(*this->fe);
1063  auto &data = dynamic_cast<InternalData &>(*data_ptr);
1064  data.initialize(this->requires_update_flags(update_flags), q, q.size());
1065 
1066  return data_ptr;
1067 }
1068 
1069 
1070 
1071 template <int dim, int spacedim>
1072 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1074  const UpdateFlags update_flags,
1075  const hp::QCollection<dim - 1> &quadrature) const
1076 {
1077  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1078  std::make_unique<InternalData>(*this->fe);
1079  auto &data = dynamic_cast<InternalData &>(*data_ptr);
1080  data.initialize_face(this->requires_update_flags(update_flags),
1082  this->fe->reference_cell(), quadrature),
1083  quadrature.max_n_quadrature_points());
1084 
1085  return data_ptr;
1086 }
1087 
1088 
1089 
1090 template <int dim, int spacedim>
1091 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
1093  const UpdateFlags update_flags,
1094  const Quadrature<dim - 1> &quadrature) const
1095 {
1096  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
1097  std::make_unique<InternalData>(*this->fe);
1098  auto &data = dynamic_cast<InternalData &>(*data_ptr);
1099  data.initialize_face(this->requires_update_flags(update_flags),
1101  this->fe->reference_cell(), quadrature),
1102  quadrature.size());
1103 
1104  return data_ptr;
1105 }
1106 
1107 
1108 
1109 template <int dim, int spacedim>
1112  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1113  const CellSimilarity::Similarity cell_similarity,
1114  const Quadrature<dim> &quadrature,
1115  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1117  &output_data) const
1118 {
1119  // ensure that the following static_cast is really correct:
1120  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1121  ExcInternalError());
1122  const InternalData &data = static_cast<const InternalData &>(internal_data);
1123 
1124  const unsigned int n_q_points = quadrature.size();
1125 
1126  // recompute the support points of the transformation of this
1127  // cell. we tried to be clever here in an earlier version of the
1128  // library by checking whether the cell is the same as the one we
1129  // had visited last, but it turns out to be difficult to determine
1130  // that because a cell for the purposes of a mapping is
1131  // characterized not just by its (triangulation, level, index)
1132  // triple, but also by the locations of its vertices, the manifold
1133  // object attached to the cell and all of its bounding faces/edges,
1134  // etc. to reliably test that the "cell" we are on is, therefore,
1135  // not easily done
1136  data.mapping_support_points = this->compute_mapping_support_points(cell);
1137  data.cell_of_current_support_points = cell;
1138 
1139  // if the order of the mapping is greater than 1, then do not reuse any cell
1140  // similarity information. This is necessary because the cell similarity
1141  // value is computed with just cell vertices and does not take into account
1142  // cell curvature.
1143  const CellSimilarity::Similarity computed_cell_similarity =
1144  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
1145 
1146  internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
1148  data,
1149  output_data.quadrature_points,
1150  n_q_points);
1151 
1152  internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
1153  computed_cell_similarity,
1155  data,
1156  n_q_points);
1157 
1158  internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1159  computed_cell_similarity,
1161  data,
1162  output_data.jacobian_grads,
1163  n_q_points);
1164 
1166  dim,
1167  spacedim>(computed_cell_similarity,
1169  data,
1170  output_data.jacobian_pushed_forward_grads,
1171  n_q_points);
1172 
1174  dim,
1175  spacedim>(computed_cell_similarity,
1177  data,
1178  output_data.jacobian_2nd_derivatives,
1179  n_q_points);
1180 
1181  internal::MappingFEImplementation::
1182  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1183  computed_cell_similarity,
1185  data,
1187  n_q_points);
1188 
1190  dim,
1191  spacedim>(computed_cell_similarity,
1193  data,
1194  output_data.jacobian_3rd_derivatives,
1195  n_q_points);
1196 
1197  internal::MappingFEImplementation::
1198  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1199  computed_cell_similarity,
1201  data,
1203  n_q_points);
1204 
1205  const UpdateFlags update_flags = data.update_each;
1206  const std::vector<double> &weights = quadrature.get_weights();
1207 
1208  // Multiply quadrature weights by absolute value of Jacobian determinants or
1209  // the area element g=sqrt(DX^t DX) in case of codim > 0
1210 
1211  if (update_flags & (update_normal_vectors | update_JxW_values))
1212  {
1213  AssertDimension(output_data.JxW_values.size(), n_q_points);
1214 
1215  Assert(!(update_flags & update_normal_vectors) ||
1216  (output_data.normal_vectors.size() == n_q_points),
1217  ExcDimensionMismatch(output_data.normal_vectors.size(),
1218  n_q_points));
1219 
1220 
1221  if (computed_cell_similarity != CellSimilarity::translation)
1222  for (unsigned int point = 0; point < n_q_points; ++point)
1223  {
1224  if (dim == spacedim)
1225  {
1226  const double det = data.contravariant[point].determinant();
1227 
1228  // check for distorted cells.
1229 
1230  // TODO: this allows for anisotropies of up to 1e6 in 3d and
1231  // 1e12 in 2d. might want to find a finer
1232  // (dimension-independent) criterion
1233  Assert(det >
1234  1e-12 * Utilities::fixed_power<dim>(
1235  cell->diameter() / std::sqrt(double(dim))),
1237  cell->center(), det, point)));
1238 
1239  output_data.JxW_values[point] = weights[point] * det;
1240  }
1241  // if dim==spacedim, then there is no cell normal to
1242  // compute. since this is for FEValues (and not FEFaceValues),
1243  // there are also no face normals to compute
1244  else // codim>0 case
1245  {
1246  Tensor<1, spacedim> DX_t[dim];
1247  for (unsigned int i = 0; i < spacedim; ++i)
1248  for (unsigned int j = 0; j < dim; ++j)
1249  DX_t[j][i] = data.contravariant[point][i][j];
1250 
1251  Tensor<2, dim> G; // First fundamental form
1252  for (unsigned int i = 0; i < dim; ++i)
1253  for (unsigned int j = 0; j < dim; ++j)
1254  G[i][j] = DX_t[i] * DX_t[j];
1255 
1256  output_data.JxW_values[point] =
1257  std::sqrt(determinant(G)) * weights[point];
1258 
1259  if (computed_cell_similarity ==
1261  {
1262  // we only need to flip the normal
1263  if (update_flags & update_normal_vectors)
1264  output_data.normal_vectors[point] *= -1.;
1265  }
1266  else
1267  {
1268  if (update_flags & update_normal_vectors)
1269  {
1270  Assert(spacedim == dim + 1,
1271  ExcMessage(
1272  "There is no (unique) cell normal for " +
1274  "-dimensional cells in " +
1275  Utilities::int_to_string(spacedim) +
1276  "-dimensional space. This only works if the "
1277  "space dimension is one greater than the "
1278  "dimensionality of the mesh cells."));
1279 
1280  if (dim == 1)
1281  output_data.normal_vectors[point] =
1282  cross_product_2d(-DX_t[0]);
1283  else // dim == 2
1284  output_data.normal_vectors[point] =
1285  cross_product_3d(DX_t[0], DX_t[1]);
1286 
1287  output_data.normal_vectors[point] /=
1288  output_data.normal_vectors[point].norm();
1289 
1290  if (cell->direction_flag() == false)
1291  output_data.normal_vectors[point] *= -1.;
1292  }
1293  }
1294  } // codim>0 case
1295  }
1296  }
1297 
1298 
1299 
1300  // copy values from InternalData to vector given by reference
1301  if (update_flags & update_jacobians)
1302  {
1303  AssertDimension(output_data.jacobians.size(), n_q_points);
1304  if (computed_cell_similarity != CellSimilarity::translation)
1305  for (unsigned int point = 0; point < n_q_points; ++point)
1306  output_data.jacobians[point] = data.contravariant[point];
1307  }
1308 
1309  // copy values from InternalData to vector given by reference
1310  if (update_flags & update_inverse_jacobians)
1311  {
1312  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1313  if (computed_cell_similarity != CellSimilarity::translation)
1314  for (unsigned int point = 0; point < n_q_points; ++point)
1315  output_data.inverse_jacobians[point] =
1316  data.covariant[point].transpose();
1317  }
1318 
1319  return computed_cell_similarity;
1320 }
1321 
1322 
1323 
1324 namespace internal
1325 {
1326  namespace MappingFEImplementation
1327  {
1328  namespace
1329  {
1340  template <int dim, int spacedim>
1341  void
1343  const ::MappingFE<dim, spacedim> &mapping,
1344  const typename ::Triangulation<dim, spacedim>::cell_iterator
1345  &cell,
1346  const unsigned int face_no,
1347  const unsigned int subface_no,
1348  const unsigned int n_q_points,
1349  const typename QProjector<dim>::DataSetDescriptor data_set,
1350  const typename ::MappingFE<dim, spacedim>::InternalData &data,
1352  &output_data)
1353  {
1354  const UpdateFlags update_flags = data.update_each;
1355 
1356  if (update_flags &
1359  {
1360  if (update_flags & update_boundary_forms)
1361  AssertIndexRange(n_q_points,
1362  output_data.boundary_forms.size() + 1);
1363  if (update_flags & update_normal_vectors)
1364  AssertIndexRange(n_q_points,
1365  output_data.normal_vectors.size() + 1);
1366  if (update_flags & update_JxW_values)
1367  AssertIndexRange(n_q_points, output_data.JxW_values.size() + 1);
1368 
1369  Assert(data.aux.size() + 1 >= dim, ExcInternalError());
1370 
1371  // first compute some common data that is used for evaluating
1372  // all of the flags below
1373 
1374  // map the unit tangentials to the real cell. checking for
1375  // d!=dim-1 eliminates compiler warnings regarding unsigned int
1376  // expressions < 0.
1377  for (unsigned int d = 0; d != dim - 1; ++d)
1378  {
1379  Assert(face_no + cell->n_faces() * d <
1380  data.unit_tangentials.size(),
1381  ExcInternalError());
1382  Assert(
1383  data.aux[d].size() <=
1384  data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1385  ExcInternalError());
1386 
1387  mapping.transform(
1389  data.unit_tangentials[face_no + cell->n_faces() * d]),
1391  data,
1392  make_array_view(data.aux[d]));
1393  }
1394 
1395  if (update_flags & update_boundary_forms)
1396  {
1397  // if dim==spacedim, we can use the unit tangentials to
1398  // compute the boundary form by simply taking the cross
1399  // product
1400  if (dim == spacedim)
1401  {
1402  for (unsigned int i = 0; i < n_q_points; ++i)
1403  switch (dim)
1404  {
1405  case 1:
1406  // in 1d, we don't have access to any of the
1407  // data.aux fields (because it has only dim-1
1408  // components), but we can still compute the
1409  // boundary form by simply looking at the number
1410  // of the face
1411  output_data.boundary_forms[i][0] =
1412  (face_no == 0 ? -1 : +1);
1413  break;
1414  case 2:
1415  output_data.boundary_forms[i] =
1416  cross_product_2d(data.aux[0][i]);
1417  break;
1418  case 3:
1419  output_data.boundary_forms[i] =
1420  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1421  break;
1422  default:
1423  Assert(false, ExcNotImplemented());
1424  }
1425  }
1426  else //(dim < spacedim)
1427  {
1428  // in the codim-one case, the boundary form results from
1429  // the cross product of all the face tangential vectors
1430  // and the cell normal vector
1431  //
1432  // to compute the cell normal, use the same method used in
1433  // fill_fe_values for cells above
1434  AssertIndexRange(n_q_points, data.contravariant.size() + 1);
1435 
1436  for (unsigned int point = 0; point < n_q_points; ++point)
1437  {
1438  if (dim == 1)
1439  {
1440  // J is a tangent vector
1441  output_data.boundary_forms[point] =
1442  data.contravariant[point].transpose()[0];
1443  output_data.boundary_forms[point] /=
1444  (face_no == 0 ? -1. : +1.) *
1445  output_data.boundary_forms[point].norm();
1446  }
1447 
1448  if (dim == 2)
1449  {
1451  data.contravariant[point].transpose();
1452 
1453  Tensor<1, spacedim> cell_normal =
1454  cross_product_3d(DX_t[0], DX_t[1]);
1455  cell_normal /= cell_normal.norm();
1456 
1457  // then compute the face normal from the face
1458  // tangent and the cell normal:
1459  output_data.boundary_forms[point] =
1460  cross_product_3d(data.aux[0][point], cell_normal);
1461  }
1462  }
1463  }
1464  }
1465 
1466  if (update_flags & update_JxW_values)
1467  for (unsigned int i = 0; i < n_q_points; ++i)
1468  {
1469  output_data.JxW_values[i] =
1470  output_data.boundary_forms[i].norm() *
1471  data.quadrature_weights[i + data_set];
1472 
1473  if (subface_no != numbers::invalid_unsigned_int)
1474  {
1475 #if false
1476  const double area_ratio =
1478  cell->subface_case(face_no), subface_no);
1479  output_data.JxW_values[i] *= area_ratio;
1480 #else
1481  Assert(false, ExcNotImplemented());
1482 #endif
1483  }
1484  }
1485 
1486  if (update_flags & update_normal_vectors)
1487  for (unsigned int i = 0; i < n_q_points; ++i)
1488  output_data.normal_vectors[i] =
1489  Point<spacedim>(output_data.boundary_forms[i] /
1490  output_data.boundary_forms[i].norm());
1491 
1492  if (update_flags & update_jacobians)
1493  for (unsigned int point = 0; point < n_q_points; ++point)
1494  output_data.jacobians[point] = data.contravariant[point];
1495 
1496  if (update_flags & update_inverse_jacobians)
1497  for (unsigned int point = 0; point < n_q_points; ++point)
1498  output_data.inverse_jacobians[point] =
1499  data.covariant[point].transpose();
1500  }
1501  }
1502 
1503 
1510  template <int dim, int spacedim>
1511  void
1513  const ::MappingFE<dim, spacedim> &mapping,
1514  const typename ::Triangulation<dim, spacedim>::cell_iterator
1515  &cell,
1516  const unsigned int face_no,
1517  const unsigned int subface_no,
1518  const typename QProjector<dim>::DataSetDescriptor data_set,
1519  const Quadrature<dim - 1> &quadrature,
1520  const typename ::MappingFE<dim, spacedim>::InternalData &data,
1522  &output_data)
1523  {
1524  const unsigned int n_q_points = quadrature.size();
1525 
1526  maybe_compute_q_points<dim, spacedim>(data_set,
1527  data,
1528  output_data.quadrature_points,
1529  n_q_points);
1530  maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
1531  data_set,
1532  data,
1533  n_q_points);
1534  maybe_update_jacobian_grads<dim, spacedim>(CellSimilarity::none,
1535  data_set,
1536  data,
1537  output_data.jacobian_grads,
1538  n_q_points);
1539  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
1541  data_set,
1542  data,
1543  output_data.jacobian_pushed_forward_grads,
1544  n_q_points);
1545  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
1547  data_set,
1548  data,
1549  output_data.jacobian_2nd_derivatives,
1550  n_q_points);
1551  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1553  data_set,
1554  data,
1556  n_q_points);
1557  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
1559  data_set,
1560  data,
1561  output_data.jacobian_3rd_derivatives,
1562  n_q_points);
1563  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1565  data_set,
1566  data,
1568  n_q_points);
1569 
1570  maybe_compute_face_data(mapping,
1571  cell,
1572  face_no,
1573  subface_no,
1574  n_q_points,
1575  data_set,
1576  data,
1577  output_data);
1578  }
1579  } // namespace
1580  } // namespace MappingFEImplementation
1581 } // namespace internal
1582 
1583 
1584 
1585 template <int dim, int spacedim>
1586 void
1588  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1589  const unsigned int face_no,
1590  const hp::QCollection<dim - 1> &quadrature,
1591  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1593  &output_data) const
1594 {
1595  // ensure that the following cast is really correct:
1596  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1597  ExcInternalError());
1598  const InternalData &data = static_cast<const InternalData &>(internal_data);
1599 
1600  // if necessary, recompute the support points of the transformation of this
1601  // cell (note that we need to first check the triangulation pointer, since
1602  // otherwise the second test might trigger an exception if the
1603  // triangulations are not the same)
1604  if ((data.mapping_support_points.empty()) ||
1605  (&cell->get_triangulation() !=
1607  (cell != data.cell_of_current_support_points))
1608  {
1609  data.mapping_support_points = this->compute_mapping_support_points(cell);
1610  data.cell_of_current_support_points = cell;
1611  }
1612 
1614  *this,
1615  cell,
1616  face_no,
1618  QProjector<dim>::DataSetDescriptor::face(this->fe->reference_cell(),
1619  face_no,
1620  cell->face_orientation(face_no),
1621  cell->face_flip(face_no),
1622  cell->face_rotation(face_no),
1623  quadrature),
1624  quadrature[quadrature.size() == 1 ? 0 : face_no],
1625  data,
1626  output_data);
1627 }
1628 
1629 
1630 
1631 template <int dim, int spacedim>
1632 void
1634  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1635  const unsigned int face_no,
1636  const unsigned int subface_no,
1637  const Quadrature<dim - 1> &quadrature,
1638  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1640  &output_data) const
1641 {
1642  // ensure that the following cast is really correct:
1643  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1644  ExcInternalError());
1645  const InternalData &data = static_cast<const InternalData &>(internal_data);
1646 
1647  // if necessary, recompute the support points of the transformation of this
1648  // cell (note that we need to first check the triangulation pointer, since
1649  // otherwise the second test might trigger an exception if the
1650  // triangulations are not the same)
1651  if ((data.mapping_support_points.empty()) ||
1652  (&cell->get_triangulation() !=
1654  (cell != data.cell_of_current_support_points))
1655  {
1656  data.mapping_support_points = this->compute_mapping_support_points(cell);
1657  data.cell_of_current_support_points = cell;
1658  }
1659 
1661  *this,
1662  cell,
1663  face_no,
1664  subface_no,
1665  QProjector<dim>::DataSetDescriptor::subface(this->fe->reference_cell(),
1666  face_no,
1667  subface_no,
1668  cell->face_orientation(face_no),
1669  cell->face_flip(face_no),
1670  cell->face_rotation(face_no),
1671  quadrature.size(),
1672  cell->subface_case(face_no)),
1673  quadrature,
1674  data,
1675  output_data);
1676 }
1677 
1678 
1679 
1680 namespace internal
1681 {
1682  namespace MappingFEImplementation
1683  {
1684  namespace
1685  {
1686  template <int dim, int spacedim, int rank>
1687  void
1689  const ArrayView<const Tensor<rank, dim>> &input,
1690  const MappingKind mapping_kind,
1691  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1692  const ArrayView<Tensor<rank, spacedim>> &output)
1693  {
1694  // In the case of wedges and pyramids, faces might have different
1695  // numbers of quadrature points on each face with the result
1696  // that input and output have different sizes, since input has
1697  // the correct size but the size of output is the maximum of
1698  // all possible sizes.
1699  AssertIndexRange(input.size(), output.size() + 1);
1700 
1701  Assert(
1702  (dynamic_cast<
1703  const typename ::MappingFE<dim, spacedim>::InternalData *>(
1704  &mapping_data) != nullptr),
1705  ExcInternalError());
1706  const typename ::MappingFE<dim, spacedim>::InternalData &data =
1707  static_cast<
1708  const typename ::MappingFE<dim, spacedim>::InternalData &>(
1709  mapping_data);
1710 
1711  switch (mapping_kind)
1712  {
1713  case mapping_contravariant:
1714  {
1715  Assert(
1716  data.update_each & update_contravariant_transformation,
1718  "update_contravariant_transformation"));
1719 
1720  for (unsigned int i = 0; i < input.size(); ++i)
1721  output[i] =
1722  apply_transformation(data.contravariant[i], input[i]);
1723 
1724  return;
1725  }
1726 
1727  case mapping_piola:
1728  {
1729  Assert(
1730  data.update_each & update_contravariant_transformation,
1732  "update_contravariant_transformation"));
1733  Assert(
1734  data.update_each & update_volume_elements,
1736  "update_volume_elements"));
1737  Assert(rank == 1, ExcMessage("Only for rank 1"));
1738  if (rank != 1)
1739  return;
1740 
1741  for (unsigned int i = 0; i < input.size(); ++i)
1742  {
1743  output[i] =
1744  apply_transformation(data.contravariant[i], input[i]);
1745  output[i] /= data.volume_elements[i];
1746  }
1747  return;
1748  }
1749  // We still allow this operation as in the
1750  // reference cell Derivatives are Tensor
1751  // rather than DerivativeForm
1752  case mapping_covariant:
1753  {
1754  Assert(
1755  data.update_each & update_contravariant_transformation,
1757  "update_covariant_transformation"));
1758 
1759  for (unsigned int i = 0; i < input.size(); ++i)
1760  output[i] = apply_transformation(data.covariant[i], input[i]);
1761 
1762  return;
1763  }
1764 
1765  default:
1766  Assert(false, ExcNotImplemented());
1767  }
1768  }
1769 
1770 
1771  template <int dim, int spacedim, int rank>
1772  void
1774  const ArrayView<const Tensor<rank, dim>> &input,
1775  const MappingKind mapping_kind,
1776  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1777  const ArrayView<Tensor<rank, spacedim>> &output)
1778  {
1779  AssertDimension(input.size(), output.size());
1780  Assert(
1781  (dynamic_cast<
1782  const typename ::MappingFE<dim, spacedim>::InternalData *>(
1783  &mapping_data) != nullptr),
1784  ExcInternalError());
1785  const typename ::MappingFE<dim, spacedim>::InternalData &data =
1786  static_cast<
1787  const typename ::MappingFE<dim, spacedim>::InternalData &>(
1788  mapping_data);
1789 
1790  switch (mapping_kind)
1791  {
1793  {
1794  Assert(
1795  data.update_each & update_covariant_transformation,
1797  "update_covariant_transformation"));
1798  Assert(
1799  data.update_each & update_contravariant_transformation,
1801  "update_contravariant_transformation"));
1802  Assert(rank == 2, ExcMessage("Only for rank 2"));
1803 
1804  for (unsigned int i = 0; i < output.size(); ++i)
1805  {
1807  apply_transformation(data.contravariant[i],
1808  transpose(input[i]));
1809  output[i] =
1810  apply_transformation(data.covariant[i], A.transpose());
1811  }
1812 
1813  return;
1814  }
1815 
1817  {
1818  Assert(
1819  data.update_each & update_covariant_transformation,
1821  "update_covariant_transformation"));
1822  Assert(rank == 2, ExcMessage("Only for rank 2"));
1823 
1824  for (unsigned int i = 0; i < output.size(); ++i)
1825  {
1827  apply_transformation(data.covariant[i],
1828  transpose(input[i]));
1829  output[i] =
1830  apply_transformation(data.covariant[i], A.transpose());
1831  }
1832 
1833  return;
1834  }
1835 
1837  {
1838  Assert(
1839  data.update_each & update_covariant_transformation,
1841  "update_covariant_transformation"));
1842  Assert(
1843  data.update_each & update_contravariant_transformation,
1845  "update_contravariant_transformation"));
1846  Assert(
1847  data.update_each & update_volume_elements,
1849  "update_volume_elements"));
1850  Assert(rank == 2, ExcMessage("Only for rank 2"));
1851 
1852  for (unsigned int i = 0; i < output.size(); ++i)
1853  {
1855  apply_transformation(data.covariant[i], input[i]);
1856  const Tensor<2, spacedim> T =
1857  apply_transformation(data.contravariant[i],
1858  A.transpose());
1859 
1860  output[i] = transpose(T);
1861  output[i] /= data.volume_elements[i];
1862  }
1863 
1864  return;
1865  }
1866 
1867  default:
1868  Assert(false, ExcNotImplemented());
1869  }
1870  }
1871 
1872 
1873 
1874  template <int dim, int spacedim>
1875  void
1877  const ArrayView<const Tensor<3, dim>> &input,
1878  const MappingKind mapping_kind,
1879  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1880  const ArrayView<Tensor<3, spacedim>> &output)
1881  {
1882  AssertDimension(input.size(), output.size());
1883  Assert(
1884  (dynamic_cast<
1885  const typename ::MappingFE<dim, spacedim>::InternalData *>(
1886  &mapping_data) != nullptr),
1887  ExcInternalError());
1888  const typename ::MappingFE<dim, spacedim>::InternalData &data =
1889  static_cast<
1890  const typename ::MappingFE<dim, spacedim>::InternalData &>(
1891  mapping_data);
1892 
1893  switch (mapping_kind)
1894  {
1896  {
1897  Assert(
1898  data.update_each & update_covariant_transformation,
1900  "update_covariant_transformation"));
1901  Assert(
1902  data.update_each & update_contravariant_transformation,
1904  "update_contravariant_transformation"));
1905 
1906  for (unsigned int q = 0; q < output.size(); ++q)
1907  for (unsigned int i = 0; i < spacedim; ++i)
1908  {
1909  double tmp1[dim][dim];
1910  for (unsigned int J = 0; J < dim; ++J)
1911  for (unsigned int K = 0; K < dim; ++K)
1912  {
1913  tmp1[J][K] =
1914  data.contravariant[q][i][0] * input[q][0][J][K];
1915  for (unsigned int I = 1; I < dim; ++I)
1916  tmp1[J][K] +=
1917  data.contravariant[q][i][I] * input[q][I][J][K];
1918  }
1919  for (unsigned int j = 0; j < spacedim; ++j)
1920  {
1921  double tmp2[dim];
1922  for (unsigned int K = 0; K < dim; ++K)
1923  {
1924  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1925  for (unsigned int J = 1; J < dim; ++J)
1926  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1927  }
1928  for (unsigned int k = 0; k < spacedim; ++k)
1929  {
1930  output[q][i][j][k] =
1931  data.covariant[q][k][0] * tmp2[0];
1932  for (unsigned int K = 1; K < dim; ++K)
1933  output[q][i][j][k] +=
1934  data.covariant[q][k][K] * tmp2[K];
1935  }
1936  }
1937  }
1938  return;
1939  }
1940 
1942  {
1943  Assert(
1944  data.update_each & update_covariant_transformation,
1946  "update_covariant_transformation"));
1947 
1948  for (unsigned int q = 0; q < output.size(); ++q)
1949  for (unsigned int i = 0; i < spacedim; ++i)
1950  {
1951  double tmp1[dim][dim];
1952  for (unsigned int J = 0; J < dim; ++J)
1953  for (unsigned int K = 0; K < dim; ++K)
1954  {
1955  tmp1[J][K] =
1956  data.covariant[q][i][0] * input[q][0][J][K];
1957  for (unsigned int I = 1; I < dim; ++I)
1958  tmp1[J][K] +=
1959  data.covariant[q][i][I] * input[q][I][J][K];
1960  }
1961  for (unsigned int j = 0; j < spacedim; ++j)
1962  {
1963  double tmp2[dim];
1964  for (unsigned int K = 0; K < dim; ++K)
1965  {
1966  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1967  for (unsigned int J = 1; J < dim; ++J)
1968  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1969  }
1970  for (unsigned int k = 0; k < spacedim; ++k)
1971  {
1972  output[q][i][j][k] =
1973  data.covariant[q][k][0] * tmp2[0];
1974  for (unsigned int K = 1; K < dim; ++K)
1975  output[q][i][j][k] +=
1976  data.covariant[q][k][K] * tmp2[K];
1977  }
1978  }
1979  }
1980 
1981  return;
1982  }
1983 
1984  case mapping_piola_hessian:
1985  {
1986  Assert(
1987  data.update_each & update_covariant_transformation,
1989  "update_covariant_transformation"));
1990  Assert(
1991  data.update_each & update_contravariant_transformation,
1993  "update_contravariant_transformation"));
1994  Assert(
1995  data.update_each & update_volume_elements,
1997  "update_volume_elements"));
1998 
1999  for (unsigned int q = 0; q < output.size(); ++q)
2000  for (unsigned int i = 0; i < spacedim; ++i)
2001  {
2002  double factor[dim];
2003  for (unsigned int I = 0; I < dim; ++I)
2004  factor[I] =
2005  data.contravariant[q][i][I] / data.volume_elements[q];
2006  double tmp1[dim][dim];
2007  for (unsigned int J = 0; J < dim; ++J)
2008  for (unsigned int K = 0; K < dim; ++K)
2009  {
2010  tmp1[J][K] = factor[0] * input[q][0][J][K];
2011  for (unsigned int I = 1; I < dim; ++I)
2012  tmp1[J][K] += factor[I] * input[q][I][J][K];
2013  }
2014  for (unsigned int j = 0; j < spacedim; ++j)
2015  {
2016  double tmp2[dim];
2017  for (unsigned int K = 0; K < dim; ++K)
2018  {
2019  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
2020  for (unsigned int J = 1; J < dim; ++J)
2021  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
2022  }
2023  for (unsigned int k = 0; k < spacedim; ++k)
2024  {
2025  output[q][i][j][k] =
2026  data.covariant[q][k][0] * tmp2[0];
2027  for (unsigned int K = 1; K < dim; ++K)
2028  output[q][i][j][k] +=
2029  data.covariant[q][k][K] * tmp2[K];
2030  }
2031  }
2032  }
2033 
2034  return;
2035  }
2036 
2037  default:
2038  Assert(false, ExcNotImplemented());
2039  }
2040  }
2041 
2042 
2043 
2044  template <int dim, int spacedim, int rank>
2045  void
2047  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2048  const MappingKind mapping_kind,
2049  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2050  const ArrayView<Tensor<rank + 1, spacedim>> &output)
2051  {
2052  AssertDimension(input.size(), output.size());
2053  Assert(
2054  (dynamic_cast<
2055  const typename ::MappingFE<dim, spacedim>::InternalData *>(
2056  &mapping_data) != nullptr),
2057  ExcInternalError());
2058  const typename ::MappingFE<dim, spacedim>::InternalData &data =
2059  static_cast<
2060  const typename ::MappingFE<dim, spacedim>::InternalData &>(
2061  mapping_data);
2062 
2063  switch (mapping_kind)
2064  {
2065  case mapping_covariant:
2066  {
2067  Assert(
2068  data.update_each & update_contravariant_transformation,
2070  "update_covariant_transformation"));
2071 
2072  for (unsigned int i = 0; i < output.size(); ++i)
2073  output[i] = apply_transformation(data.covariant[i], input[i]);
2074 
2075  return;
2076  }
2077  default:
2078  Assert(false, ExcNotImplemented());
2079  }
2080  }
2081  } // namespace
2082  } // namespace MappingFEImplementation
2083 } // namespace internal
2084 
2085 
2086 
2087 template <int dim, int spacedim>
2088 void
2090  const ArrayView<const Tensor<1, dim>> &input,
2091  const MappingKind mapping_kind,
2092  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2093  const ArrayView<Tensor<1, spacedim>> &output) const
2094 {
2096  mapping_kind,
2097  mapping_data,
2098  output);
2099 }
2100 
2101 
2102 
2103 template <int dim, int spacedim>
2104 void
2106  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2107  const MappingKind mapping_kind,
2108  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2109  const ArrayView<Tensor<2, spacedim>> &output) const
2110 {
2112  mapping_kind,
2113  mapping_data,
2114  output);
2115 }
2116 
2117 
2118 
2119 template <int dim, int spacedim>
2120 void
2122  const ArrayView<const Tensor<2, dim>> &input,
2123  const MappingKind mapping_kind,
2124  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2125  const ArrayView<Tensor<2, spacedim>> &output) const
2126 {
2127  switch (mapping_kind)
2128  {
2129  case mapping_contravariant:
2131  mapping_kind,
2132  mapping_data,
2133  output);
2134  return;
2135 
2140  mapping_kind,
2141  mapping_data,
2142  output);
2143  return;
2144  default:
2145  Assert(false, ExcNotImplemented());
2146  }
2147 }
2148 
2149 
2150 
2151 template <int dim, int spacedim>
2152 void
2154  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2155  const MappingKind mapping_kind,
2156  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2157  const ArrayView<Tensor<3, spacedim>> &output) const
2158 {
2159  AssertDimension(input.size(), output.size());
2160  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2161  ExcInternalError());
2162  const InternalData &data = static_cast<const InternalData &>(mapping_data);
2163 
2164  switch (mapping_kind)
2165  {
2167  {
2168  Assert(data.update_each & update_contravariant_transformation,
2170  "update_covariant_transformation"));
2171 
2172  for (unsigned int q = 0; q < output.size(); ++q)
2173  for (unsigned int i = 0; i < spacedim; ++i)
2174  for (unsigned int j = 0; j < spacedim; ++j)
2175  {
2176  double tmp[dim];
2177  for (unsigned int K = 0; K < dim; ++K)
2178  {
2179  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
2180  for (unsigned int J = 1; J < dim; ++J)
2181  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
2182  }
2183  for (unsigned int k = 0; k < spacedim; ++k)
2184  {
2185  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
2186  for (unsigned int K = 1; K < dim; ++K)
2187  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
2188  }
2189  }
2190  return;
2191  }
2192 
2193  default:
2194  Assert(false, ExcNotImplemented());
2195  }
2196 }
2197 
2198 
2199 
2200 template <int dim, int spacedim>
2201 void
2203  const ArrayView<const Tensor<3, dim>> &input,
2204  const MappingKind mapping_kind,
2205  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2206  const ArrayView<Tensor<3, spacedim>> &output) const
2207 {
2208  switch (mapping_kind)
2209  {
2210  case mapping_piola_hessian:
2214  mapping_kind,
2215  mapping_data,
2216  output);
2217  return;
2218  default:
2219  Assert(false, ExcNotImplemented());
2220  }
2221 }
2222 
2223 
2224 
2225 namespace
2226 {
2227  template <int spacedim>
2228  bool
2229  check_all_manifold_ids_identical(
2231  {
2232  return true;
2233  }
2234 
2235 
2236 
2237  template <int spacedim>
2238  bool
2239  check_all_manifold_ids_identical(
2241  {
2242  const auto m_id = cell->manifold_id();
2243 
2244  for (const auto f : cell->face_indices())
2245  if (m_id != cell->face(f)->manifold_id())
2246  return false;
2247 
2248  return true;
2249  }
2250 
2251 
2252 
2253  template <int spacedim>
2254  bool
2255  check_all_manifold_ids_identical(
2257  {
2258  const auto m_id = cell->manifold_id();
2259 
2260  for (const auto f : cell->face_indices())
2261  if (m_id != cell->face(f)->manifold_id())
2262  return false;
2263 
2264  for (const auto l : cell->line_indices())
2265  if (m_id != cell->line(l)->manifold_id())
2266  return false;
2267 
2268  return true;
2269  }
2270 } // namespace
2271 
2272 
2273 
2274 template <int dim, int spacedim>
2275 std::vector<Point<spacedim>>
2277  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2278 {
2279  Assert(
2280  check_all_manifold_ids_identical(cell),
2281  ExcMessage(
2282  "All entities of a cell need to have the same manifold id as the cell has."));
2283 
2284  std::vector<Point<spacedim>> vertices(cell->n_vertices());
2285 
2286  for (const unsigned int i : cell->vertex_indices())
2287  vertices[i] = cell->vertex(i);
2288 
2289  std::vector<Point<spacedim>> mapping_support_points(
2290  fe->get_unit_support_points().size());
2291 
2292  cell->get_manifold().get_new_points(vertices,
2293  mapping_support_point_weights,
2294  mapping_support_points);
2295 
2296  return mapping_support_points;
2297 }
2298 
2299 
2300 
2301 template <int dim, int spacedim>
2304  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
2305 {
2306  return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
2307 }
2308 
2309 
2310 
2311 template <int dim, int spacedim>
2312 bool
2314  const ReferenceCell &reference_cell) const
2315 {
2316  Assert(dim == reference_cell.get_dimension(),
2317  ExcMessage("The dimension of your mapping (" +
2318  Utilities::to_string(dim) +
2319  ") and the reference cell cell_type (" +
2320  Utilities::to_string(reference_cell.get_dimension()) +
2321  " ) do not agree."));
2322 
2323  return fe->reference_cell() == reference_cell;
2324 }
2325 
2326 
2327 
2328 //--------------------------- Explicit instantiations -----------------------
2329 #include "mapping_fe.inst"
2330 
2331 
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:950
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
DerivativeForm< 1, spacedim, dim, Number > transpose() const
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition: mapping_fe.h:388
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_fe.cc:138
std::vector< Point< spacedim > > mapping_support_points
Definition: mapping_fe.h:382
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_fe.cc:80
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition: mapping_fe.h:371
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Definition: mapping_fe.h:362
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
Definition: mapping_fe.cc:178
virtual std::size_t memory_consumption() const override
Definition: mapping_fe.cc:60
InternalData(const FiniteElement< dim, spacedim > &fe)
Definition: mapping_fe.cc:49
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_fe.cc:1111
const unsigned int polynomial_degree
Definition: mapping_fe.h:467
MappingFE(const FiniteElement< dim, spacedim > &fe)
Definition: mapping_fe.cc:843
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_fe.cc:1073
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
Definition: mapping_fe.cc:921
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_fe.cc:1002
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_fe.cc:2089
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_fe.cc:1587
Table< 2, double > mapping_support_point_weights
Definition: mapping_fe.h:477
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_fe.cc:2276
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_fe.cc:2303
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_fe.cc:1058
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition: mapping_fe.cc:1092
unsigned int get_degree() const
Definition: mapping_fe.cc:893
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition: mapping_fe.cc:2313
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_fe.cc:1633
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_fe.cc:884
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
Definition: mapping_fe.cc:902
const std::unique_ptr< FiniteElement< dim, spacedim > > fe
Definition: mapping_fe.h:461
Abstract base class for mapping classes.
Definition: mapping.h:317
Definition: point.h:112
static DataSetDescriptor cell()
Definition: qprojector.h:394
const std::vector< Point< dim > > & get_points() const
const std::vector< double > & get_weights() const
unsigned int size() const
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
unsigned int size() const
Definition: collection.h:265
unsigned int max_n_quadrature_points() const
Definition: q_collection.h:174
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
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< Tensor< 1, spacedim > > boundary_forms
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Point< 2 > second
Definition: grid_out.cc:4615
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
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_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:78
@ mapping_piola
Definition: mapping.h:113
@ mapping_covariant_gradient
Definition: mapping.h:99
@ mapping_covariant
Definition: mapping.h:88
@ mapping_contravariant
Definition: mapping.h:93
@ mapping_contravariant_hessian
Definition: mapping.h:155
@ mapping_covariant_hessian
Definition: mapping.h:149
@ mapping_contravariant_gradient
Definition: mapping.h:105
@ mapping_piola_gradient
Definition: mapping.h:119
@ mapping_piola_hessian
Definition: mapping.h:161
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char A
static const char T
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:490
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
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
void transform_gradients(const ArrayView< const Tensor< rank, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim >> &output)
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim >> &output)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 4, spacedim >> &jacobian_pushed_forward_2nd_derivatives)
void maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 5, spacedim >> &jacobian_pushed_forward_3rd_derivatives)
void maybe_compute_face_data(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int n_q_points, const std::vector< double > &weights, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
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)
static const unsigned int invalid_unsigned_int
Definition: types.h:221
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)