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