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