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