Reference documentation for deal.II version 9.2.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_field.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 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 #include <deal.II/base/array_view.h>
17 #include <deal.II/base/memory_consumption.h>
18 #include <deal.II/base/polynomial.h>
19 #include <deal.II/base/qprojector.h>
20 #include <deal.II/base/quadrature.h>
21 #include <deal.II/base/quadrature_lib.h>
22 #include <deal.II/base/std_cxx14/memory.h>
23 #include <deal.II/base/tensor_product_polynomials.h>
24 #include <deal.II/base/utilities.h>
25 
26 #include <deal.II/dofs/dof_accessor.h>
27 
28 #include <deal.II/fe/fe_q.h>
29 #include <deal.II/fe/fe_system.h>
30 #include <deal.II/fe/fe_tools.h>
31 #include <deal.II/fe/fe_values.h>
32 #include <deal.II/fe/mapping.h>
33 #include <deal.II/fe/mapping_fe_field.h>
34 #include <deal.II/fe/mapping_q1.h>
35 
36 #include <deal.II/grid/tria_iterator.h>
37 
38 #include <deal.II/lac/block_vector.h>
39 #include <deal.II/lac/full_matrix.h>
40 #include <deal.II/lac/la_parallel_block_vector.h>
41 #include <deal.II/lac/la_parallel_vector.h>
42 #include <deal.II/lac/la_vector.h>
43 #include <deal.II/lac/petsc_block_vector.h>
44 #include <deal.II/lac/petsc_vector.h>
45 #include <deal.II/lac/trilinos_epetra_vector.h>
46 #include <deal.II/lac/trilinos_parallel_block_vector.h>
47 #include <deal.II/lac/trilinos_tpetra_vector.h>
48 #include <deal.II/lac/trilinos_vector.h>
49 #include <deal.II/lac/vector.h>
50 
51 #include <deal.II/numerics/vector_tools.h>
52 
53 #include <fstream>
54 #include <memory>
55 #include <numeric>
56 
57 
58 
59 DEAL_II_NAMESPACE_OPEN
60 
61 
62 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
63 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::InternalData::
64  InternalData(const FiniteElement<dim, spacedim> &fe,
65  const ComponentMask & mask)
66  : unit_tangentials()
67  , n_shape_functions(fe.dofs_per_cell)
68  , mask(mask)
69  , local_dof_indices(fe.dofs_per_cell)
70  , local_dof_values(fe.dofs_per_cell)
71 {}
72 
73 
74 
75 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
76 std::size_t
77 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::InternalData::
78  memory_consumption() const
79 {
80  Assert(false, ExcNotImplemented());
81  return 0;
82 }
83 
84 
85 
86 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
87 double &
88 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::shape(
89  const unsigned int qpoint,
90  const unsigned int shape_nr)
91 {
92  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
93  return shape_values[qpoint * n_shape_functions + shape_nr];
94 }
95 
96 
97 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
98 const Tensor<1, dim> &
99 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
100  derivative(const unsigned int qpoint, const unsigned int shape_nr) const
101 {
102  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
103  shape_derivatives.size());
104  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
105 }
106 
107 
108 
109 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
110 Tensor<1, dim> &
111 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
112  derivative(const unsigned int qpoint, const unsigned int shape_nr)
113 {
114  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
115  shape_derivatives.size());
116  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
117 }
118 
119 
120 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
121 const Tensor<2, dim> &
122 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
123  second_derivative(const unsigned int qpoint,
124  const unsigned int shape_nr) const
125 {
126  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
127  shape_second_derivatives.size());
128  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
129 }
130 
131 
132 
133 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
134 Tensor<2, dim> &
135 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
136  second_derivative(const unsigned int qpoint, const unsigned int shape_nr)
137 {
138  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
139  shape_second_derivatives.size());
140  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
141 }
142 
143 
144 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
145 const Tensor<3, dim> &
146 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
147  third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
148 {
149  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
150  shape_third_derivatives.size());
151  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
152 }
153 
154 
155 
156 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
157 Tensor<3, dim> &
158 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
159  third_derivative(const unsigned int qpoint, const unsigned int shape_nr)
160 {
161  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
162  shape_third_derivatives.size());
163  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
164 }
165 
166 
167 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
168 const Tensor<4, dim> &
169 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
170  fourth_derivative(const unsigned int qpoint,
171  const unsigned int shape_nr) const
172 {
173  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
174  shape_fourth_derivatives.size());
175  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
176 }
177 
178 
179 
180 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
181 Tensor<4, dim> &
182 MappingFEField<dim, spacedim, DoFHandlerType, VectorType>::InternalData::
183  fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr)
184 {
185  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
186  shape_fourth_derivatives.size());
187  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
188 }
189 
190 
191 
192 namespace
193 {
194  // Construct a quadrature formula containing the vertices of the reference
195  // cell in dimension dim (with invalid weights)
196  template <int dim>
197  Quadrature<dim> &
198  get_vertex_quadrature()
199  {
200  static Quadrature<dim> quad;
201  if (quad.size() == 0)
202  {
203  std::vector<Point<dim>> points(GeometryInfo<dim>::vertices_per_cell);
204  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
205  points[i] = GeometryInfo<dim>::unit_cell_vertex(i);
206  quad = Quadrature<dim>(points);
207  }
208  return quad;
209  }
210 } // namespace
211 
212 
213 
214 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
215 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::MappingFEField(
216  const DoFHandlerType &euler_dof_handler,
217  const VectorType & euler_vector,
218  const ComponentMask & mask)
219  : uses_level_dofs(false)
220  , euler_vector({&euler_vector})
221  , euler_dof_handler(&euler_dof_handler)
222  , fe_mask(mask.size() ?
223  mask :
224  ComponentMask(
225  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
226  true))
227  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
228  , fe_values(this->euler_dof_handler->get_fe(),
229  get_vertex_quadrature<dim>(),
230  update_values)
231 {
232  unsigned int size = 0;
233  for (unsigned int i = 0; i < fe_mask.size(); ++i)
234  {
235  if (fe_mask[i])
236  fe_to_real[i] = size++;
237  }
238  AssertDimension(size, spacedim);
239 }
240 
241 
242 
243 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
244 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::MappingFEField(
245  const DoFHandlerType & euler_dof_handler,
246  const std::vector<VectorType> &euler_vector,
247  const ComponentMask & mask)
248  : uses_level_dofs(true)
249  , euler_dof_handler(&euler_dof_handler)
250  , fe_mask(mask.size() ?
251  mask :
252  ComponentMask(
253  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
254  true))
255  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
256  , fe_values(this->euler_dof_handler->get_fe(),
257  get_vertex_quadrature<dim>(),
258  update_values)
259 {
260  unsigned int size = 0;
261  for (unsigned int i = 0; i < fe_mask.size(); ++i)
262  {
263  if (fe_mask[i])
264  fe_to_real[i] = size++;
265  }
266  AssertDimension(size, spacedim);
267 
268  Assert(euler_dof_handler.has_level_dofs(),
269  ExcMessage("The underlying DoFHandler object did not call "
270  "distribute_mg_dofs(). In this case, the construction via "
271  "level vectors does not make sense."));
272  AssertDimension(euler_vector.size(),
273  euler_dof_handler.get_triangulation().n_global_levels());
274  this->euler_vector.clear();
275  this->euler_vector.resize(euler_vector.size());
276  for (unsigned int i = 0; i < euler_vector.size(); ++i)
277  this->euler_vector[i] = &euler_vector[i];
278 }
279 
280 
281 
282 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
283 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::MappingFEField(
284  const DoFHandlerType & euler_dof_handler,
285  const MGLevelObject<VectorType> &euler_vector,
286  const ComponentMask & mask)
287  : uses_level_dofs(true)
288  , euler_dof_handler(&euler_dof_handler)
289  , fe_mask(mask.size() ?
290  mask :
291  ComponentMask(
292  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
293  true))
294  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
295  , fe_values(this->euler_dof_handler->get_fe(),
296  get_vertex_quadrature<dim>(),
297  update_values)
298 {
299  unsigned int size = 0;
300  for (unsigned int i = 0; i < fe_mask.size(); ++i)
301  {
302  if (fe_mask[i])
303  fe_to_real[i] = size++;
304  }
305  AssertDimension(size, spacedim);
306 
307  Assert(euler_dof_handler.has_level_dofs(),
308  ExcMessage("The underlying DoFHandler object did not call "
309  "distribute_mg_dofs(). In this case, the construction via "
310  "level vectors does not make sense."));
311  AssertDimension(euler_vector.max_level() + 1,
312  euler_dof_handler.get_triangulation().n_global_levels());
313  this->euler_vector.clear();
314  this->euler_vector.resize(
315  euler_dof_handler.get_triangulation().n_global_levels());
316  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
317  ++i)
318  this->euler_vector[i] = &euler_vector[i];
319 }
320 
321 
322 
323 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
324 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::MappingFEField(
325  const MappingFEField<dim, spacedim, VectorType, DoFHandlerType> &mapping)
326  : uses_level_dofs(mapping.uses_level_dofs)
327  , euler_vector(mapping.euler_vector)
328  , euler_dof_handler(mapping.euler_dof_handler)
329  , fe_mask(mapping.fe_mask)
330  , fe_to_real(mapping.fe_to_real)
331  , fe_values(mapping.euler_dof_handler->get_fe(),
332  get_vertex_quadrature<dim>(),
333  update_values)
334 {}
335 
336 
337 
338 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
339 inline const double &
340 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::InternalData::shape(
341  const unsigned int qpoint,
342  const unsigned int shape_nr) const
343 {
344  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
345  return shape_values[qpoint * n_shape_functions + shape_nr];
346 }
347 
348 
349 
350 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
351 bool
352 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
353  preserves_vertex_locations() const
354 {
355  return false;
356 }
357 
358 
359 
360 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
361 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
362 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_vertices(
363  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
364 {
365  // we transform our tria iterator into a dof iterator so we can access
366  // data not associated with triangulations
367  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
368  *cell, euler_dof_handler);
369 
370  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
371  AssertDimension(GeometryInfo<dim>::vertices_per_cell,
372  fe_values.n_quadrature_points);
373  AssertDimension(fe_to_real.size(),
374  euler_dof_handler->get_fe().n_components());
375  if (uses_level_dofs)
376  {
377  AssertIndexRange(cell->level(), euler_vector.size());
378  AssertDimension(euler_vector[cell->level()]->size(),
379  euler_dof_handler->n_dofs(cell->level()));
380  }
381  else
382  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
383 
384  {
385  std::lock_guard<std::mutex> lock(fe_values_mutex);
386  fe_values.reinit(dof_cell);
387  }
388  const unsigned int dofs_per_cell = euler_dof_handler->get_fe().dofs_per_cell;
389  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
390  if (uses_level_dofs)
391  dof_cell->get_mg_dof_indices(dof_indices);
392  else
393  dof_cell->get_dof_indices(dof_indices);
394 
395  const VectorType &vector =
396  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
397 
398  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
399  for (unsigned int i = 0; i < dofs_per_cell; ++i)
400  {
401  const unsigned int comp = fe_to_real
402  [euler_dof_handler->get_fe().system_to_component_index(i).first];
403  if (comp != numbers::invalid_unsigned_int)
404  {
405  typename VectorType::value_type value =
406  internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
407  if (euler_dof_handler->get_fe().is_primitive(i))
408  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
409  vertices[v][comp] += fe_values.shape_value(i, v) * value;
410  else
411  Assert(false, ExcNotImplemented());
412  }
413  }
414 
415  return vertices;
416 }
417 
418 
419 
420 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
421 void
422 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
423  compute_shapes_virtual(
424  const std::vector<Point<dim>> &unit_points,
425  typename MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
426  InternalData &data) const
427 {
428  const auto fe = &euler_dof_handler->get_fe();
429  const unsigned int n_points = unit_points.size();
430 
431  for (unsigned int point = 0; point < n_points; ++point)
432  {
433  if (data.shape_values.size() != 0)
434  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
435  data.shape(point, i) = fe->shape_value(i, unit_points[point]);
436 
437  if (data.shape_derivatives.size() != 0)
438  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
439  data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
440 
441  if (data.shape_second_derivatives.size() != 0)
442  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
443  data.second_derivative(point, i) =
444  fe->shape_grad_grad(i, unit_points[point]);
445 
446  if (data.shape_third_derivatives.size() != 0)
447  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
448  data.third_derivative(point, i) =
449  fe->shape_3rd_derivative(i, unit_points[point]);
450 
451  if (data.shape_fourth_derivatives.size() != 0)
452  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
453  data.fourth_derivative(point, i) =
454  fe->shape_4th_derivative(i, unit_points[point]);
455  }
456 }
457 
458 
459 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
460 UpdateFlags
461 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
462  requires_update_flags(const UpdateFlags in) const
463 {
464  // add flags if the respective quantities are necessary to compute
465  // what we need. note that some flags appear in both conditions and
466  // in subsequent set operations. this leads to some circular
467  // logic. the only way to treat this is to iterate. since there are
468  // 5 if-clauses in the loop, it will take at most 4 iterations to
469  // converge. do them:
470  UpdateFlags out = in;
471  for (unsigned int i = 0; i < 5; ++i)
472  {
473  // The following is a little incorrect:
474  // If not applied on a face,
475  // update_boundary_forms does not
476  // make sense. On the other hand,
477  // it is necessary on a
478  // face. Currently,
479  // update_boundary_forms is simply
480  // ignored for the interior of a
481  // cell.
482  if (out & (update_JxW_values | update_normal_vectors))
483  out |= update_boundary_forms;
484 
485  if (out & (update_covariant_transformation | update_JxW_values |
486  update_jacobians | update_jacobian_grads |
487  update_boundary_forms | update_normal_vectors))
488  out |= update_contravariant_transformation;
489 
490  if (out &
491  (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
492  update_jacobian_pushed_forward_2nd_derivatives |
493  update_jacobian_pushed_forward_3rd_derivatives))
494  out |= update_covariant_transformation;
495 
496  // The contravariant transformation
497  // is a Piola transformation, which
498  // requires the determinant of the
499  // Jacobi matrix of the transformation.
500  // Therefore these values have to be
501  // updated for each cell.
502  if (out & update_contravariant_transformation)
503  out |= update_JxW_values;
504 
505  if (out & update_normal_vectors)
506  out |= update_JxW_values;
507  }
508 
509  return out;
510 }
511 
512 
513 
514 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
515 void
516 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::compute_data(
517  const UpdateFlags update_flags,
518  const Quadrature<dim> &q,
519  const unsigned int n_original_q_points,
520  InternalData & data) const
521 {
522  // store the flags in the internal data object so we can access them
523  // in fill_fe_*_values(). use the transitive hull of the required
524  // flags
525  data.update_each = requires_update_flags(update_flags);
526 
527  const unsigned int n_q_points = q.size();
528 
529  // see if we need the (transformation) shape function values
530  // and/or gradients and resize the necessary arrays
531  if (data.update_each & update_quadrature_points)
532  data.shape_values.resize(data.n_shape_functions * n_q_points);
533 
534  if (data.update_each &
535  (update_covariant_transformation | update_contravariant_transformation |
536  update_JxW_values | update_boundary_forms | update_normal_vectors |
537  update_jacobians | update_jacobian_grads | update_inverse_jacobians))
538  data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
539 
540  if (data.update_each & update_covariant_transformation)
541  data.covariant.resize(n_original_q_points);
542 
543  if (data.update_each & update_contravariant_transformation)
544  data.contravariant.resize(n_original_q_points);
545 
546  if (data.update_each & update_volume_elements)
547  data.volume_elements.resize(n_original_q_points);
548 
549  if (data.update_each &
550  (update_jacobian_grads | update_jacobian_pushed_forward_grads))
551  data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
552 
553  if (data.update_each & (update_jacobian_2nd_derivatives |
554  update_jacobian_pushed_forward_2nd_derivatives))
555  data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
556 
557  if (data.update_each & (update_jacobian_3rd_derivatives |
558  update_jacobian_pushed_forward_3rd_derivatives))
559  data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
560 
561  compute_shapes_virtual(q.get_points(), data);
562 }
563 
564 
565 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
566 void
567 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::compute_face_data(
568  const UpdateFlags update_flags,
569  const Quadrature<dim> &q,
570  const unsigned int n_original_q_points,
571  InternalData & data) const
572 {
573  compute_data(update_flags, q, n_original_q_points, data);
574 
575  if (dim > 1)
576  {
577  if (data.update_each & update_boundary_forms)
578  {
579  data.aux.resize(
580  dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
581 
582  // Compute tangentials to the unit cell.
583  for (const unsigned int i : GeometryInfo<dim>::face_indices())
584  {
585  data.unit_tangentials[i].resize(n_original_q_points);
586  std::fill(data.unit_tangentials[i].begin(),
587  data.unit_tangentials[i].end(),
588  GeometryInfo<dim>::unit_tangential_vectors[i][0]);
589  if (dim > 2)
590  {
591  data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
592  .resize(n_original_q_points);
593  std::fill(
594  data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
595  .begin(),
596  data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
597  .end(),
598  GeometryInfo<dim>::unit_tangential_vectors[i][1]);
599  }
600  }
601  }
602  }
603 }
604 
605 
606 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
607 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
608 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_data(
609  const UpdateFlags update_flags,
610  const Quadrature<dim> &quadrature) const
611 {
612  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
613  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
614  auto &data = dynamic_cast<InternalData &>(*data_ptr);
615  this->compute_data(update_flags, quadrature, quadrature.size(), data);
616 
617  return data_ptr;
618 }
619 
620 
621 
622 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
623 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
624 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_face_data(
625  const UpdateFlags update_flags,
626  const Quadrature<dim - 1> &quadrature) const
627 {
628  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
629  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
630  auto & data = dynamic_cast<InternalData &>(*data_ptr);
631  const Quadrature<dim> q(QProjector<dim>::project_to_all_faces(quadrature));
632  this->compute_face_data(update_flags, q, quadrature.size(), data);
633 
634  return data_ptr;
635 }
636 
637 
638 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
639 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
640 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_subface_data(
641  const UpdateFlags update_flags,
642  const Quadrature<dim - 1> &quadrature) const
643 {
644  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
645  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
646  auto & data = dynamic_cast<InternalData &>(*data_ptr);
647  const Quadrature<dim> q(QProjector<dim>::project_to_all_subfaces(quadrature));
648  this->compute_face_data(update_flags, q, quadrature.size(), data);
649 
650  return data_ptr;
651 }
652 
653 
654 
655 namespace internal
656 {
657  namespace MappingFEFieldImplementation
658  {
659  namespace
660  {
667  template <int dim,
668  int spacedim,
669  typename VectorType,
670  typename DoFHandlerType>
671  void
672  maybe_compute_q_points(
673  const typename ::QProjector<dim>::DataSetDescriptor data_set,
674  const typename ::
675  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
676  InternalData & data,
677  const FiniteElement<dim, spacedim> &fe,
678  const ComponentMask & fe_mask,
679  const std::vector<unsigned int> & fe_to_real,
680  std::vector<Point<spacedim>> & quadrature_points)
681  {
682  const UpdateFlags update_flags = data.update_each;
683 
684  if (update_flags & update_quadrature_points)
685  {
686  for (unsigned int point = 0; point < quadrature_points.size();
687  ++point)
688  {
689  Point<spacedim> result;
690  const double * shape = &data.shape(point + data_set, 0);
691 
692  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
693  {
694  const unsigned int comp_k =
695  fe.system_to_component_index(k).first;
696  if (fe_mask[comp_k])
697  result[fe_to_real[comp_k]] +=
698  data.local_dof_values[k] * shape[k];
699  }
700 
701  quadrature_points[point] = result;
702  }
703  }
704  }
705 
713  template <int dim,
714  int spacedim,
715  typename VectorType,
716  typename DoFHandlerType>
717  void
718  maybe_update_Jacobians(
719  const typename ::QProjector<dim>::DataSetDescriptor data_set,
720  const typename ::
721  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
722  InternalData & data,
723  const FiniteElement<dim, spacedim> &fe,
724  const ComponentMask & fe_mask,
725  const std::vector<unsigned int> & fe_to_real)
726  {
727  const UpdateFlags update_flags = data.update_each;
728 
729  // then Jacobians
730  if (update_flags & update_contravariant_transformation)
731  {
732  const unsigned int n_q_points = data.contravariant.size();
733 
734  Assert(data.n_shape_functions > 0, ExcInternalError());
735 
736  for (unsigned int point = 0; point < n_q_points; ++point)
737  {
738  const Tensor<1, dim> *data_derv =
739  &data.derivative(point + data_set, 0);
740 
741  Tensor<1, dim> result[spacedim];
742 
743  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
744  {
745  const unsigned int comp_k =
746  fe.system_to_component_index(k).first;
747  if (fe_mask[comp_k])
748  result[fe_to_real[comp_k]] +=
749  data.local_dof_values[k] * data_derv[k];
750  }
751 
752  // write result into contravariant data
753  for (unsigned int i = 0; i < spacedim; ++i)
754  {
755  data.contravariant[point][i] = result[i];
756  }
757  }
758  }
759 
760  if (update_flags & update_covariant_transformation)
761  {
762  AssertDimension(data.covariant.size(), data.contravariant.size());
763  for (unsigned int point = 0; point < data.contravariant.size();
764  ++point)
765  data.covariant[point] =
766  (data.contravariant[point]).covariant_form();
767  }
768 
769  if (update_flags & update_volume_elements)
770  {
771  AssertDimension(data.contravariant.size(),
772  data.volume_elements.size());
773  for (unsigned int point = 0; point < data.contravariant.size();
774  ++point)
775  data.volume_elements[point] =
776  data.contravariant[point].determinant();
777  }
778  }
779 
786  template <int dim,
787  int spacedim,
788  typename VectorType,
789  typename DoFHandlerType>
790  void
791  maybe_update_jacobian_grads(
792  const typename ::QProjector<dim>::DataSetDescriptor data_set,
793  const typename ::
794  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
795  InternalData & data,
796  const FiniteElement<dim, spacedim> & fe,
797  const ComponentMask & fe_mask,
798  const std::vector<unsigned int> & fe_to_real,
799  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
800  {
801  const UpdateFlags update_flags = data.update_each;
802  if (update_flags & update_jacobian_grads)
803  {
804  const unsigned int n_q_points = jacobian_grads.size();
805 
806  for (unsigned int point = 0; point < n_q_points; ++point)
807  {
808  const Tensor<2, dim> *second =
809  &data.second_derivative(point + data_set, 0);
810 
811  DerivativeForm<2, dim, spacedim> result;
812 
813  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
814  {
815  const unsigned int comp_k =
816  fe.system_to_component_index(k).first;
817  if (fe_mask[comp_k])
818  for (unsigned int j = 0; j < dim; ++j)
819  for (unsigned int l = 0; l < dim; ++l)
820  result[fe_to_real[comp_k]][j][l] +=
821  (second[k][j][l] * data.local_dof_values[k]);
822  }
823 
824  // never touch any data for j=dim in case dim<spacedim, so
825  // it will always be zero as it was initialized
826  for (unsigned int i = 0; i < spacedim; ++i)
827  for (unsigned int j = 0; j < dim; ++j)
828  for (unsigned int l = 0; l < dim; ++l)
829  jacobian_grads[point][i][j][l] = result[i][j][l];
830  }
831  }
832  }
833 
840  template <int dim,
841  int spacedim,
842  typename VectorType,
843  typename DoFHandlerType>
844  void
845  maybe_update_jacobian_pushed_forward_grads(
846  const typename ::QProjector<dim>::DataSetDescriptor data_set,
847  const typename ::
848  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
849  InternalData & data,
850  const FiniteElement<dim, spacedim> &fe,
851  const ComponentMask & fe_mask,
852  const std::vector<unsigned int> & fe_to_real,
853  std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
854  {
855  const UpdateFlags update_flags = data.update_each;
856  if (update_flags & update_jacobian_pushed_forward_grads)
857  {
858  const unsigned int n_q_points =
859  jacobian_pushed_forward_grads.size();
860 
861  double tmp[spacedim][spacedim][spacedim];
862  for (unsigned int point = 0; point < n_q_points; ++point)
863  {
864  const Tensor<2, dim> *second =
865  &data.second_derivative(point + data_set, 0);
866 
867  DerivativeForm<2, dim, spacedim> result;
868 
869  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
870  {
871  const unsigned int comp_k =
872  fe.system_to_component_index(k).first;
873  if (fe_mask[comp_k])
874  for (unsigned int j = 0; j < dim; ++j)
875  for (unsigned int l = 0; l < dim; ++l)
876  result[fe_to_real[comp_k]][j][l] +=
877  (second[k][j][l] * data.local_dof_values[k]);
878  }
879 
880  // first push forward the j-components
881  for (unsigned int i = 0; i < spacedim; ++i)
882  for (unsigned int j = 0; j < spacedim; ++j)
883  for (unsigned int l = 0; l < dim; ++l)
884  {
885  tmp[i][j][l] =
886  result[i][0][l] * data.covariant[point][j][0];
887  for (unsigned int jr = 1; jr < dim; ++jr)
888  {
889  tmp[i][j][l] +=
890  result[i][jr][l] * data.covariant[point][j][jr];
891  }
892  }
893 
894  // now, pushing forward the l-components
895  for (unsigned int i = 0; i < spacedim; ++i)
896  for (unsigned int j = 0; j < spacedim; ++j)
897  for (unsigned int l = 0; l < spacedim; ++l)
898  {
899  jacobian_pushed_forward_grads[point][i][j][l] =
900  tmp[i][j][0] * data.covariant[point][l][0];
901  for (unsigned int lr = 1; lr < dim; ++lr)
902  {
903  jacobian_pushed_forward_grads[point][i][j][l] +=
904  tmp[i][j][lr] * data.covariant[point][l][lr];
905  }
906  }
907  }
908  }
909  }
910 
917  template <int dim,
918  int spacedim,
919  typename VectorType,
920  typename DoFHandlerType>
921  void
922  maybe_update_jacobian_2nd_derivatives(
923  const typename ::QProjector<dim>::DataSetDescriptor data_set,
924  const typename ::
925  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
926  InternalData & data,
927  const FiniteElement<dim, spacedim> & fe,
928  const ComponentMask & fe_mask,
929  const std::vector<unsigned int> & fe_to_real,
930  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
931  {
932  const UpdateFlags update_flags = data.update_each;
933  if (update_flags & update_jacobian_2nd_derivatives)
934  {
935  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
936 
937  for (unsigned int point = 0; point < n_q_points; ++point)
938  {
939  const Tensor<3, dim> *third =
940  &data.third_derivative(point + data_set, 0);
941 
942  DerivativeForm<3, dim, spacedim> result;
943 
944  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
945  {
946  const unsigned int comp_k =
947  fe.system_to_component_index(k).first;
948  if (fe_mask[comp_k])
949  for (unsigned int j = 0; j < dim; ++j)
950  for (unsigned int l = 0; l < dim; ++l)
951  for (unsigned int m = 0; m < dim; ++m)
952  result[fe_to_real[comp_k]][j][l][m] +=
953  (third[k][j][l][m] * data.local_dof_values[k]);
954  }
955 
956  // never touch any data for j=dim in case dim<spacedim, so
957  // it will always be zero as it was initialized
958  for (unsigned int i = 0; i < spacedim; ++i)
959  for (unsigned int j = 0; j < dim; ++j)
960  for (unsigned int l = 0; l < dim; ++l)
961  for (unsigned int m = 0; m < dim; ++m)
962  jacobian_2nd_derivatives[point][i][j][l][m] =
963  result[i][j][l][m];
964  }
965  }
966  }
967 
975  template <int dim,
976  int spacedim,
977  typename VectorType,
978  typename DoFHandlerType>
979  void
980  maybe_update_jacobian_pushed_forward_2nd_derivatives(
981  const typename ::QProjector<dim>::DataSetDescriptor data_set,
982  const typename ::
983  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
984  InternalData & data,
985  const FiniteElement<dim, spacedim> &fe,
986  const ComponentMask & fe_mask,
987  const std::vector<unsigned int> & fe_to_real,
988  std::vector<Tensor<4, spacedim>>
989  &jacobian_pushed_forward_2nd_derivatives)
990  {
991  const UpdateFlags update_flags = data.update_each;
992  if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
993  {
994  const unsigned int n_q_points =
995  jacobian_pushed_forward_2nd_derivatives.size();
996 
997  double tmp[spacedim][spacedim][spacedim][spacedim];
998  for (unsigned int point = 0; point < n_q_points; ++point)
999  {
1000  const Tensor<3, dim> *third =
1001  &data.third_derivative(point + data_set, 0);
1002 
1003  DerivativeForm<3, dim, spacedim> result;
1004 
1005  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1006  {
1007  const unsigned int comp_k =
1008  fe.system_to_component_index(k).first;
1009  if (fe_mask[comp_k])
1010  for (unsigned int j = 0; j < dim; ++j)
1011  for (unsigned int l = 0; l < dim; ++l)
1012  for (unsigned int m = 0; m < dim; ++m)
1013  result[fe_to_real[comp_k]][j][l][m] +=
1014  (third[k][j][l][m] * data.local_dof_values[k]);
1015  }
1016 
1017  // push forward the j-coordinate
1018  for (unsigned int i = 0; i < spacedim; ++i)
1019  for (unsigned int j = 0; j < spacedim; ++j)
1020  for (unsigned int l = 0; l < dim; ++l)
1021  for (unsigned int m = 0; m < dim; ++m)
1022  {
1023  jacobian_pushed_forward_2nd_derivatives
1024  [point][i][j][l][m] =
1025  result[i][0][l][m] * data.covariant[point][j][0];
1026  for (unsigned int jr = 1; jr < dim; ++jr)
1027  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1028  [l][m] +=
1029  result[i][jr][l][m] *
1030  data.covariant[point][j][jr];
1031  }
1032 
1033  // push forward the l-coordinate
1034  for (unsigned int i = 0; i < spacedim; ++i)
1035  for (unsigned int j = 0; j < spacedim; ++j)
1036  for (unsigned int l = 0; l < spacedim; ++l)
1037  for (unsigned int m = 0; m < dim; ++m)
1038  {
1039  tmp[i][j][l][m] =
1040  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1041  [0][m] *
1042  data.covariant[point][l][0];
1043  for (unsigned int lr = 1; lr < dim; ++lr)
1044  tmp[i][j][l][m] +=
1045  jacobian_pushed_forward_2nd_derivatives[point][i]
1046  [j][lr]
1047  [m] *
1048  data.covariant[point][l][lr];
1049  }
1050 
1051  // push forward the m-coordinate
1052  for (unsigned int i = 0; i < spacedim; ++i)
1053  for (unsigned int j = 0; j < spacedim; ++j)
1054  for (unsigned int l = 0; l < spacedim; ++l)
1055  for (unsigned int m = 0; m < spacedim; ++m)
1056  {
1057  jacobian_pushed_forward_2nd_derivatives
1058  [point][i][j][l][m] =
1059  tmp[i][j][l][0] * data.covariant[point][m][0];
1060  for (unsigned int mr = 1; mr < dim; ++mr)
1061  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1062  [l][m] +=
1063  tmp[i][j][l][mr] * data.covariant[point][m][mr];
1064  }
1065  }
1066  }
1067  }
1068 
1075  template <int dim,
1076  int spacedim,
1077  typename VectorType,
1078  typename DoFHandlerType>
1079  void
1080  maybe_update_jacobian_3rd_derivatives(
1081  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1082  const typename ::
1083  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1084  InternalData & data,
1085  const FiniteElement<dim, spacedim> & fe,
1086  const ComponentMask & fe_mask,
1087  const std::vector<unsigned int> & fe_to_real,
1088  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1089  {
1090  const UpdateFlags update_flags = data.update_each;
1091  if (update_flags & update_jacobian_3rd_derivatives)
1092  {
1093  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1094 
1095  for (unsigned int point = 0; point < n_q_points; ++point)
1096  {
1097  const Tensor<4, dim> *fourth =
1098  &data.fourth_derivative(point + data_set, 0);
1099 
1100  DerivativeForm<4, dim, spacedim> result;
1101 
1102  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1103  {
1104  const unsigned int comp_k =
1105  fe.system_to_component_index(k).first;
1106  if (fe_mask[comp_k])
1107  for (unsigned int j = 0; j < dim; ++j)
1108  for (unsigned int l = 0; l < dim; ++l)
1109  for (unsigned int m = 0; m < dim; ++m)
1110  for (unsigned int n = 0; n < dim; ++n)
1111  result[fe_to_real[comp_k]][j][l][m][n] +=
1112  (fourth[k][j][l][m][n] *
1113  data.local_dof_values[k]);
1114  }
1115 
1116  // never touch any data for j,l,m,n=dim in case
1117  // dim<spacedim, so it will always be zero as it was
1118  // initialized
1119  for (unsigned int i = 0; i < spacedim; ++i)
1120  for (unsigned int j = 0; j < dim; ++j)
1121  for (unsigned int l = 0; l < dim; ++l)
1122  for (unsigned int m = 0; m < dim; ++m)
1123  for (unsigned int n = 0; n < dim; ++n)
1124  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1125  result[i][j][l][m][n];
1126  }
1127  }
1128  }
1129 
1137  template <int dim,
1138  int spacedim,
1139  typename VectorType,
1140  typename DoFHandlerType>
1141  void
1142  maybe_update_jacobian_pushed_forward_3rd_derivatives(
1143  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1144  const typename ::
1145  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1146  InternalData & data,
1147  const FiniteElement<dim, spacedim> &fe,
1148  const ComponentMask & fe_mask,
1149  const std::vector<unsigned int> & fe_to_real,
1150  std::vector<Tensor<5, spacedim>>
1151  &jacobian_pushed_forward_3rd_derivatives)
1152  {
1153  const UpdateFlags update_flags = data.update_each;
1154  if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
1155  {
1156  const unsigned int n_q_points =
1157  jacobian_pushed_forward_3rd_derivatives.size();
1158 
1159  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1160  for (unsigned int point = 0; point < n_q_points; ++point)
1161  {
1162  const Tensor<4, dim> *fourth =
1163  &data.fourth_derivative(point + data_set, 0);
1164 
1165  DerivativeForm<4, dim, spacedim> result;
1166 
1167  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1168  {
1169  const unsigned int comp_k =
1170  fe.system_to_component_index(k).first;
1171  if (fe_mask[comp_k])
1172  for (unsigned int j = 0; j < dim; ++j)
1173  for (unsigned int l = 0; l < dim; ++l)
1174  for (unsigned int m = 0; m < dim; ++m)
1175  for (unsigned int n = 0; n < dim; ++n)
1176  result[fe_to_real[comp_k]][j][l][m][n] +=
1177  (fourth[k][j][l][m][n] *
1178  data.local_dof_values[k]);
1179  }
1180 
1181  // push-forward the j-coordinate
1182  for (unsigned int i = 0; i < spacedim; ++i)
1183  for (unsigned int j = 0; j < spacedim; ++j)
1184  for (unsigned int l = 0; l < dim; ++l)
1185  for (unsigned int m = 0; m < dim; ++m)
1186  for (unsigned int n = 0; n < dim; ++n)
1187  {
1188  tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1189  data.covariant[point][j][0];
1190  for (unsigned int jr = 1; jr < dim; ++jr)
1191  tmp[i][j][l][m][n] +=
1192  result[i][jr][l][m][n] *
1193  data.covariant[point][j][jr];
1194  }
1195 
1196  // push-forward the l-coordinate
1197  for (unsigned int i = 0; i < spacedim; ++i)
1198  for (unsigned int j = 0; j < spacedim; ++j)
1199  for (unsigned int l = 0; l < spacedim; ++l)
1200  for (unsigned int m = 0; m < dim; ++m)
1201  for (unsigned int n = 0; n < dim; ++n)
1202  {
1203  jacobian_pushed_forward_3rd_derivatives
1204  [point][i][j][l][m][n] =
1205  tmp[i][j][0][m][n] *
1206  data.covariant[point][l][0];
1207  for (unsigned int lr = 1; lr < dim; ++lr)
1208  jacobian_pushed_forward_3rd_derivatives[point][i]
1209  [j][l][m]
1210  [n] +=
1211  tmp[i][j][lr][m][n] *
1212  data.covariant[point][l][lr];
1213  }
1214 
1215  // push-forward the m-coordinate
1216  for (unsigned int i = 0; i < spacedim; ++i)
1217  for (unsigned int j = 0; j < spacedim; ++j)
1218  for (unsigned int l = 0; l < spacedim; ++l)
1219  for (unsigned int m = 0; m < spacedim; ++m)
1220  for (unsigned int n = 0; n < dim; ++n)
1221  {
1222  tmp[i][j][l][m][n] =
1223  jacobian_pushed_forward_3rd_derivatives[point][i]
1224  [j][l][0]
1225  [n] *
1226  data.covariant[point][m][0];
1227  for (unsigned int mr = 1; mr < dim; ++mr)
1228  tmp[i][j][l][m][n] +=
1229  jacobian_pushed_forward_3rd_derivatives[point]
1230  [i][j][l]
1231  [mr][n] *
1232  data.covariant[point][m][mr];
1233  }
1234 
1235  // push-forward the n-coordinate
1236  for (unsigned int i = 0; i < spacedim; ++i)
1237  for (unsigned int j = 0; j < spacedim; ++j)
1238  for (unsigned int l = 0; l < spacedim; ++l)
1239  for (unsigned int m = 0; m < spacedim; ++m)
1240  for (unsigned int n = 0; n < spacedim; ++n)
1241  {
1242  jacobian_pushed_forward_3rd_derivatives
1243  [point][i][j][l][m][n] =
1244  tmp[i][j][l][m][0] *
1245  data.covariant[point][n][0];
1246  for (unsigned int nr = 1; nr < dim; ++nr)
1247  jacobian_pushed_forward_3rd_derivatives[point][i]
1248  [j][l][m]
1249  [n] +=
1250  tmp[i][j][l][m][nr] *
1251  data.covariant[point][n][nr];
1252  }
1253  }
1254  }
1255  }
1256 
1257 
1267  template <int dim,
1268  int spacedim,
1269  typename VectorType,
1270  typename DoFHandlerType>
1271  void
1272  maybe_compute_face_data(
1273  const ::Mapping<dim, spacedim> &mapping,
1274  const typename ::Triangulation<dim, spacedim>::cell_iterator
1275  & cell,
1276  const unsigned int face_no,
1277  const unsigned int subface_no,
1278  const std::vector<double> &weights,
1279  const typename ::
1280  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1281  InternalData &data,
1282  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1283  &output_data)
1284  {
1285  const UpdateFlags update_flags = data.update_each;
1286 
1287  if (update_flags & update_boundary_forms)
1288  {
1289  const unsigned int n_q_points = output_data.boundary_forms.size();
1290  if (update_flags & update_normal_vectors)
1291  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1292  if (update_flags & update_JxW_values)
1293  AssertDimension(output_data.JxW_values.size(), n_q_points);
1294 
1295  // map the unit tangentials to the real cell. checking for d!=dim-1
1296  // eliminates compiler warnings regarding unsigned int expressions <
1297  // 0.
1298  for (unsigned int d = 0; d != dim - 1; ++d)
1299  {
1300  Assert(face_no + GeometryInfo<dim>::faces_per_cell * d <
1301  data.unit_tangentials.size(),
1302  ExcInternalError());
1303  Assert(
1304  data.aux[d].size() <=
1305  data
1306  .unit_tangentials[face_no +
1307  GeometryInfo<dim>::faces_per_cell * d]
1308  .size(),
1309  ExcInternalError());
1310 
1311  mapping.transform(
1312  make_array_view(
1313  data
1314  .unit_tangentials[face_no +
1315  GeometryInfo<dim>::faces_per_cell * d]),
1316  mapping_contravariant,
1317  data,
1318  make_array_view(data.aux[d]));
1319  }
1320 
1321  // if dim==spacedim, we can use the unit tangentials to compute the
1322  // boundary form by simply taking the cross product
1323  if (dim == spacedim)
1324  {
1325  for (unsigned int i = 0; i < n_q_points; ++i)
1326  switch (dim)
1327  {
1328  case 1:
1329  // in 1d, we don't have access to any of the data.aux
1330  // fields (because it has only dim-1 components), but we
1331  // can still compute the boundary form by simply looking
1332  // at the number of the face
1333  output_data.boundary_forms[i][0] =
1334  (face_no == 0 ? -1 : +1);
1335  break;
1336  case 2:
1337  output_data.boundary_forms[i] =
1338  cross_product_2d(data.aux[0][i]);
1339  break;
1340  case 3:
1341  output_data.boundary_forms[i] =
1342  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1343  break;
1344  default:
1345  Assert(false, ExcNotImplemented());
1346  }
1347  }
1348  else //(dim < spacedim)
1349  {
1350  // in the codim-one case, the boundary form results from the
1351  // cross product of all the face tangential vectors and the cell
1352  // normal vector
1353  //
1354  // to compute the cell normal, use the same method used in
1355  // fill_fe_values for cells above
1356  AssertDimension(data.contravariant.size(), n_q_points);
1357 
1358  for (unsigned int point = 0; point < n_q_points; ++point)
1359  {
1360  if (dim == 1)
1361  {
1362  // J is a tangent vector
1363  output_data.boundary_forms[point] =
1364  data.contravariant[point].transpose()[0];
1365  output_data.boundary_forms[point] /=
1366  (face_no == 0 ? -1. : +1.) *
1367  output_data.boundary_forms[point].norm();
1368  }
1369 
1370  if (dim == 2)
1371  {
1372  const DerivativeForm<1, spacedim, dim> DX_t =
1373  data.contravariant[point].transpose();
1374 
1375  Tensor<1, spacedim> cell_normal =
1376  cross_product_3d(DX_t[0], DX_t[1]);
1377  cell_normal /= cell_normal.norm();
1378 
1379  // then compute the face normal from the face tangent
1380  // and the cell normal:
1381  output_data.boundary_forms[point] =
1382  cross_product_3d(data.aux[0][point], cell_normal);
1383  }
1384  }
1385  }
1386 
1387  if (update_flags & (update_normal_vectors | update_JxW_values))
1388  for (unsigned int i = 0; i < output_data.boundary_forms.size();
1389  ++i)
1390  {
1391  if (update_flags & update_JxW_values)
1392  {
1393  output_data.JxW_values[i] =
1394  output_data.boundary_forms[i].norm() * weights[i];
1395 
1396  if (subface_no != numbers::invalid_unsigned_int)
1397  {
1398  const double area_ratio =
1399  GeometryInfo<dim>::subface_ratio(
1400  cell->subface_case(face_no), subface_no);
1401  output_data.JxW_values[i] *= area_ratio;
1402  }
1403  }
1404 
1405  if (update_flags & update_normal_vectors)
1406  output_data.normal_vectors[i] =
1407  Point<spacedim>(output_data.boundary_forms[i] /
1408  output_data.boundary_forms[i].norm());
1409  }
1410 
1411  if (update_flags & update_jacobians)
1412  for (unsigned int point = 0; point < n_q_points; ++point)
1413  output_data.jacobians[point] = data.contravariant[point];
1414 
1415  if (update_flags & update_inverse_jacobians)
1416  for (unsigned int point = 0; point < n_q_points; ++point)
1417  output_data.inverse_jacobians[point] =
1418  data.covariant[point].transpose();
1419  }
1420  }
1421 
1428  template <int dim,
1429  int spacedim,
1430  typename VectorType,
1431  typename DoFHandlerType>
1432  void
1433  do_fill_fe_face_values(
1434  const ::Mapping<dim, spacedim> &mapping,
1435  const typename ::Triangulation<dim, spacedim>::cell_iterator
1436  & cell,
1437  const unsigned int face_no,
1438  const unsigned int subface_no,
1439  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1440  const Quadrature<dim - 1> & quadrature,
1441  const typename ::
1442  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1443  InternalData & data,
1444  const FiniteElement<dim, spacedim> &fe,
1445  const ComponentMask & fe_mask,
1446  const std::vector<unsigned int> & fe_to_real,
1447  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1448  &output_data)
1449  {
1450  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1451  data_set,
1452  data,
1453  fe,
1454  fe_mask,
1455  fe_to_real,
1456  output_data.quadrature_points);
1457 
1458  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1459  data_set, data, fe, fe_mask, fe_to_real);
1460 
1461  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1462  data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1463 
1464  maybe_update_jacobian_pushed_forward_grads<dim,
1465  spacedim,
1466  VectorType,
1467  DoFHandlerType>(
1468  data_set,
1469  data,
1470  fe,
1471  fe_mask,
1472  fe_to_real,
1473  output_data.jacobian_pushed_forward_grads);
1474 
1475  maybe_update_jacobian_2nd_derivatives<dim,
1476  spacedim,
1477  VectorType,
1478  DoFHandlerType>(
1479  data_set,
1480  data,
1481  fe,
1482  fe_mask,
1483  fe_to_real,
1484  output_data.jacobian_2nd_derivatives);
1485 
1486  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1487  spacedim,
1488  VectorType,
1489  DoFHandlerType>(
1490  data_set,
1491  data,
1492  fe,
1493  fe_mask,
1494  fe_to_real,
1495  output_data.jacobian_pushed_forward_2nd_derivatives);
1496 
1497  maybe_update_jacobian_3rd_derivatives<dim,
1498  spacedim,
1499  VectorType,
1500  DoFHandlerType>(
1501  data_set,
1502  data,
1503  fe,
1504  fe_mask,
1505  fe_to_real,
1506  output_data.jacobian_3rd_derivatives);
1507 
1508  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1509  spacedim,
1510  VectorType,
1511  DoFHandlerType>(
1512  data_set,
1513  data,
1514  fe,
1515  fe_mask,
1516  fe_to_real,
1517  output_data.jacobian_pushed_forward_3rd_derivatives);
1518 
1519  maybe_compute_face_data<dim, spacedim, VectorType, DoFHandlerType>(
1520  mapping,
1521  cell,
1522  face_no,
1523  subface_no,
1524  quadrature.get_weights(),
1525  data,
1526  output_data);
1527  }
1528  } // namespace
1529  } // namespace MappingFEFieldImplementation
1530 } // namespace internal
1531 
1532 
1533 // Note that the CellSimilarity flag is modifiable, since MappingFEField can
1534 // need to recalculate data even when cells are similar.
1535 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1536 CellSimilarity::Similarity
1537 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::fill_fe_values(
1538  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1539  const CellSimilarity::Similarity,
1540  const Quadrature<dim> & quadrature,
1541  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1542  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1543  &output_data) const
1544 {
1545  // convert data object to internal data for this class. fails with an
1546  // exception if that is not possible
1547  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1548  ExcInternalError());
1549  const InternalData &data = static_cast<const InternalData &>(internal_data);
1550 
1551  const unsigned int n_q_points = quadrature.size();
1552 
1553  update_internal_dofs(cell, data);
1554 
1555  internal::MappingFEFieldImplementation::
1556  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1557  QProjector<dim>::DataSetDescriptor::cell(),
1558  data,
1559  euler_dof_handler->get_fe(),
1560  fe_mask,
1561  fe_to_real,
1562  output_data.quadrature_points);
1563 
1564  internal::MappingFEFieldImplementation::
1565  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1566  QProjector<dim>::DataSetDescriptor::cell(),
1567  data,
1568  euler_dof_handler->get_fe(),
1569  fe_mask,
1570  fe_to_real);
1571 
1572  const UpdateFlags update_flags = data.update_each;
1573  const std::vector<double> &weights = quadrature.get_weights();
1574 
1575  // Multiply quadrature weights by absolute value of Jacobian determinants or
1576  // the area element g=sqrt(DX^t DX) in case of codim > 0
1577 
1578  if (update_flags & (update_normal_vectors | update_JxW_values))
1579  {
1580  AssertDimension(output_data.JxW_values.size(), n_q_points);
1581 
1582  Assert(!(update_flags & update_normal_vectors) ||
1583  (output_data.normal_vectors.size() == n_q_points),
1584  ExcDimensionMismatch(output_data.normal_vectors.size(),
1585  n_q_points));
1586 
1587 
1588  for (unsigned int point = 0; point < n_q_points; ++point)
1589  {
1590  if (dim == spacedim)
1591  {
1592  const double det = data.contravariant[point].determinant();
1593 
1594  // check for distorted cells.
1595 
1596  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1597  // 1e12 in 2D. might want to find a finer
1598  // (dimension-independent) criterion
1599  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1600  cell->diameter() / std::sqrt(double(dim))),
1601  (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1602  cell->center(), det, point)));
1603  output_data.JxW_values[point] = weights[point] * det;
1604  }
1605  // if dim==spacedim, then there is no cell normal to
1606  // compute. since this is for FEValues (and not FEFaceValues),
1607  // there are also no face normals to compute
1608  else // codim>0 case
1609  {
1610  Tensor<1, spacedim> DX_t[dim];
1611  for (unsigned int i = 0; i < spacedim; ++i)
1612  for (unsigned int j = 0; j < dim; ++j)
1613  DX_t[j][i] = data.contravariant[point][i][j];
1614 
1615  Tensor<2, dim> G; // First fundamental form
1616  for (unsigned int i = 0; i < dim; ++i)
1617  for (unsigned int j = 0; j < dim; ++j)
1618  G[i][j] = DX_t[i] * DX_t[j];
1619 
1620  output_data.JxW_values[point] =
1621  std::sqrt(determinant(G)) * weights[point];
1622 
1623  if (update_flags & update_normal_vectors)
1624  {
1625  Assert(spacedim - dim == 1,
1626  ExcMessage("There is no cell normal in codim 2."));
1627 
1628  if (dim == 1)
1629  output_data.normal_vectors[point] =
1630  cross_product_2d(-DX_t[0]);
1631  else // dim == 2
1632  output_data.normal_vectors[point] =
1633  cross_product_3d(DX_t[0], DX_t[1]);
1634 
1635  output_data.normal_vectors[point] /=
1636  output_data.normal_vectors[point].norm();
1637 
1638  if (cell->direction_flag() == false)
1639  output_data.normal_vectors[point] *= -1.;
1640  }
1641  } // codim>0 case
1642  }
1643  }
1644 
1645  // copy values from InternalData to vector given by reference
1646  if (update_flags & update_jacobians)
1647  {
1648  AssertDimension(output_data.jacobians.size(), n_q_points);
1649  for (unsigned int point = 0; point < n_q_points; ++point)
1650  output_data.jacobians[point] = data.contravariant[point];
1651  }
1652 
1653  // copy values from InternalData to vector given by reference
1654  if (update_flags & update_inverse_jacobians)
1655  {
1656  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1657  for (unsigned int point = 0; point < n_q_points; ++point)
1658  output_data.inverse_jacobians[point] =
1659  data.covariant[point].transpose();
1660  }
1661 
1662  // calculate derivatives of the Jacobians
1663  internal::MappingFEFieldImplementation::
1664  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1665  QProjector<dim>::DataSetDescriptor::cell(),
1666  data,
1667  euler_dof_handler->get_fe(),
1668  fe_mask,
1669  fe_to_real,
1670  output_data.jacobian_grads);
1671 
1672  // calculate derivatives of the Jacobians pushed forward to real cell
1673  // coordinates
1674  internal::MappingFEFieldImplementation::
1675  maybe_update_jacobian_pushed_forward_grads<dim,
1676  spacedim,
1677  VectorType,
1678  DoFHandlerType>(
1679  QProjector<dim>::DataSetDescriptor::cell(),
1680  data,
1681  euler_dof_handler->get_fe(),
1682  fe_mask,
1683  fe_to_real,
1684  output_data.jacobian_pushed_forward_grads);
1685 
1686  // calculate hessians of the Jacobians
1687  internal::MappingFEFieldImplementation::maybe_update_jacobian_2nd_derivatives<
1688  dim,
1689  spacedim,
1690  VectorType,
1691  DoFHandlerType>(QProjector<dim>::DataSetDescriptor::cell(),
1692  data,
1693  euler_dof_handler->get_fe(),
1694  fe_mask,
1695  fe_to_real,
1696  output_data.jacobian_2nd_derivatives);
1697 
1698  // calculate hessians of the Jacobians pushed forward to real cell coordinates
1699  internal::MappingFEFieldImplementation::
1700  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1701  spacedim,
1702  VectorType,
1703  DoFHandlerType>(
1704  QProjector<dim>::DataSetDescriptor::cell(),
1705  data,
1706  euler_dof_handler->get_fe(),
1707  fe_mask,
1708  fe_to_real,
1709  output_data.jacobian_pushed_forward_2nd_derivatives);
1710 
1711  // calculate gradients of the hessians of the Jacobians
1712  internal::MappingFEFieldImplementation::maybe_update_jacobian_3rd_derivatives<
1713  dim,
1714  spacedim,
1715  VectorType,
1716  DoFHandlerType>(QProjector<dim>::DataSetDescriptor::cell(),
1717  data,
1718  euler_dof_handler->get_fe(),
1719  fe_mask,
1720  fe_to_real,
1721  output_data.jacobian_3rd_derivatives);
1722 
1723  // calculate gradients of the hessians of the Jacobians pushed forward to real
1724  // cell coordinates
1725  internal::MappingFEFieldImplementation::
1726  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1727  spacedim,
1728  VectorType,
1729  DoFHandlerType>(
1730  QProjector<dim>::DataSetDescriptor::cell(),
1731  data,
1732  euler_dof_handler->get_fe(),
1733  fe_mask,
1734  fe_to_real,
1735  output_data.jacobian_pushed_forward_3rd_derivatives);
1736 
1737  return CellSimilarity::invalid_next_cell;
1738 }
1739 
1740 
1741 
1742 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1743 void
1744 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::fill_fe_face_values(
1745  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1746  const unsigned int face_no,
1747  const Quadrature<dim - 1> & quadrature,
1748  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1749  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1750  &output_data) const
1751 {
1752  // convert data object to internal data for this class. fails with an
1753  // exception if that is not possible
1754  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1755  ExcInternalError());
1756  const InternalData &data = static_cast<const InternalData &>(internal_data);
1757 
1758  update_internal_dofs(cell, data);
1759 
1760  internal::MappingFEFieldImplementation::
1761  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1762  *this,
1763  cell,
1764  face_no,
1765  numbers::invalid_unsigned_int,
1766  QProjector<dim>::DataSetDescriptor::face(face_no,
1767  cell->face_orientation(face_no),
1768  cell->face_flip(face_no),
1769  cell->face_rotation(face_no),
1770  quadrature.size()),
1771  quadrature,
1772  data,
1773  euler_dof_handler->get_fe(),
1774  fe_mask,
1775  fe_to_real,
1776  output_data);
1777 }
1778 
1779 
1780 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1781 void
1782 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1783  fill_fe_subface_values(
1784  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1785  const unsigned int face_no,
1786  const unsigned int subface_no,
1787  const Quadrature<dim - 1> & quadrature,
1788  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1789  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1790  &output_data) const
1791 {
1792  // convert data object to internal data for this class. fails with an
1793  // exception if that is not possible
1794  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1795  ExcInternalError());
1796  const InternalData &data = static_cast<const InternalData &>(internal_data);
1797 
1798  update_internal_dofs(cell, data);
1799 
1800  internal::MappingFEFieldImplementation::
1801  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1802  *this,
1803  cell,
1804  face_no,
1805  numbers::invalid_unsigned_int,
1806  QProjector<dim>::DataSetDescriptor::subface(face_no,
1807  subface_no,
1808  cell->face_orientation(
1809  face_no),
1810  cell->face_flip(face_no),
1811  cell->face_rotation(face_no),
1812  quadrature.size(),
1813  cell->subface_case(face_no)),
1814  quadrature,
1815  data,
1816  euler_dof_handler->get_fe(),
1817  fe_mask,
1818  fe_to_real,
1819  output_data);
1820 }
1821 
1822 
1823 namespace internal
1824 {
1825  namespace MappingFEFieldImplementation
1826  {
1827  namespace
1828  {
1829  template <int dim,
1830  int spacedim,
1831  int rank,
1832  typename VectorType,
1833  typename DoFHandlerType>
1834  void
1835  transform_fields(
1836  const ArrayView<const Tensor<rank, dim>> & input,
1837  const MappingKind mapping_kind,
1838  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1839  const ArrayView<Tensor<rank, spacedim>> & output)
1840  {
1841  AssertDimension(input.size(), output.size());
1842  Assert((dynamic_cast<
1843  const typename ::
1844  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1845  InternalData *>(&mapping_data) != nullptr),
1846  ExcInternalError());
1847  const typename ::MappingFEField<dim,
1848  spacedim,
1849  VectorType,
1850  DoFHandlerType>::InternalData
1851  &data = static_cast<
1852  const typename ::
1853  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1854  InternalData &>(mapping_data);
1855 
1856  switch (mapping_kind)
1857  {
1858  case mapping_contravariant:
1859  {
1860  Assert(
1861  data.update_each & update_contravariant_transformation,
1862  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1863  "update_contravariant_transformation"));
1864 
1865  for (unsigned int i = 0; i < output.size(); ++i)
1866  output[i] =
1867  apply_transformation(data.contravariant[i], input[i]);
1868 
1869  return;
1870  }
1871 
1872  case mapping_piola:
1873  {
1874  Assert(
1875  data.update_each & update_contravariant_transformation,
1876  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1877  "update_contravariant_transformation"));
1878  Assert(
1879  data.update_each & update_volume_elements,
1880  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1881  "update_volume_elements"));
1882  Assert(rank == 1, ExcMessage("Only for rank 1"));
1883  for (unsigned int i = 0; i < output.size(); ++i)
1884  {
1885  output[i] =
1886  apply_transformation(data.contravariant[i], input[i]);
1887  output[i] /= data.volume_elements[i];
1888  }
1889  return;
1890  }
1891 
1892 
1893  // We still allow this operation as in the
1894  // reference cell Derivatives are Tensor
1895  // rather than DerivativeForm
1896  case mapping_covariant:
1897  {
1898  Assert(
1899  data.update_each & update_contravariant_transformation,
1900  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1901  "update_contravariant_transformation"));
1902 
1903  for (unsigned int i = 0; i < output.size(); ++i)
1904  output[i] = apply_transformation(data.covariant[i], input[i]);
1905 
1906  return;
1907  }
1908 
1909  default:
1910  Assert(false, ExcNotImplemented());
1911  }
1912  }
1913 
1914 
1915  template <int dim,
1916  int spacedim,
1917  int rank,
1918  typename VectorType,
1919  typename DoFHandlerType>
1920  void
1921  transform_differential_forms(
1922  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1923  const MappingKind mapping_kind,
1924  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1925  const ArrayView<Tensor<rank + 1, spacedim>> & output)
1926  {
1927  AssertDimension(input.size(), output.size());
1928  Assert((dynamic_cast<
1929  const typename ::
1930  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1931  InternalData *>(&mapping_data) != nullptr),
1932  ExcInternalError());
1933  const typename ::MappingFEField<dim,
1934  spacedim,
1935  VectorType,
1936  DoFHandlerType>::InternalData
1937  &data = static_cast<
1938  const typename ::
1939  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1940  InternalData &>(mapping_data);
1941 
1942  switch (mapping_kind)
1943  {
1944  case mapping_covariant:
1945  {
1946  Assert(
1947  data.update_each & update_contravariant_transformation,
1948  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1949  "update_contravariant_transformation"));
1950 
1951  for (unsigned int i = 0; i < output.size(); ++i)
1952  output[i] = apply_transformation(data.covariant[i], input[i]);
1953 
1954  return;
1955  }
1956  default:
1957  Assert(false, ExcNotImplemented());
1958  }
1959  }
1960  } // namespace
1961  } // namespace MappingFEFieldImplementation
1962 } // namespace internal
1963 
1964 
1965 
1966 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1967 void
1968 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::transform(
1969  const ArrayView<const Tensor<1, dim>> & input,
1970  const MappingKind mapping_kind,
1971  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1972  const ArrayView<Tensor<1, spacedim>> & output) const
1973 {
1974  AssertDimension(input.size(), output.size());
1975 
1976  internal::MappingFEFieldImplementation::
1977  transform_fields<dim, spacedim, 1, VectorType, DoFHandlerType>(input,
1978  mapping_kind,
1979  mapping_data,
1980  output);
1981 }
1982 
1983 
1984 
1985 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1986 void
1987 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::transform(
1988  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1989  const MappingKind mapping_kind,
1990  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1991  const ArrayView<Tensor<2, spacedim>> & output) const
1992 {
1993  AssertDimension(input.size(), output.size());
1994 
1995  internal::MappingFEFieldImplementation::
1996  transform_differential_forms<dim, spacedim, 1, VectorType, DoFHandlerType>(
1997  input, mapping_kind, mapping_data, output);
1998 }
1999 
2000 
2001 
2002 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2003 void
2004 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::transform(
2005  const ArrayView<const Tensor<2, dim>> &input,
2006  const MappingKind,
2007  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2008  const ArrayView<Tensor<2, spacedim>> & output) const
2009 {
2010  (void)input;
2011  (void)output;
2012  (void)mapping_data;
2013  AssertDimension(input.size(), output.size());
2014 
2015  AssertThrow(false, ExcNotImplemented());
2016 }
2017 
2018 
2019 
2020 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2021 void
2022 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::transform(
2023  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2024  const MappingKind mapping_kind,
2025  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2026  const ArrayView<Tensor<3, spacedim>> & output) const
2027 {
2028  AssertDimension(input.size(), output.size());
2029  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2030  ExcInternalError());
2031  const InternalData &data = static_cast<const InternalData &>(mapping_data);
2032 
2033  switch (mapping_kind)
2034  {
2035  case mapping_covariant_gradient:
2036  {
2037  Assert(data.update_each & update_contravariant_transformation,
2038  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2039  "update_covariant_transformation"));
2040 
2041  for (unsigned int q = 0; q < output.size(); ++q)
2042  for (unsigned int i = 0; i < spacedim; ++i)
2043  for (unsigned int j = 0; j < spacedim; ++j)
2044  for (unsigned int k = 0; k < spacedim; ++k)
2045  {
2046  output[q][i][j][k] = data.covariant[q][j][0] *
2047  data.covariant[q][k][0] *
2048  input[q][i][0][0];
2049  for (unsigned int J = 0; J < dim; ++J)
2050  {
2051  const unsigned int K0 = (0 == J) ? 1 : 0;
2052  for (unsigned int K = K0; K < dim; ++K)
2053  output[q][i][j][k] += data.covariant[q][j][J] *
2054  data.covariant[q][k][K] *
2055  input[q][i][J][K];
2056  }
2057  }
2058  return;
2059  }
2060 
2061  default:
2062  Assert(false, ExcNotImplemented());
2063  }
2064 }
2065 
2066 
2067 
2068 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2069 void
2070 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::transform(
2071  const ArrayView<const Tensor<3, dim>> &input,
2072  const MappingKind /*mapping_kind*/,
2073  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2074  const ArrayView<Tensor<3, spacedim>> & output) const
2075 {
2076  (void)input;
2077  (void)output;
2078  (void)mapping_data;
2079  AssertDimension(input.size(), output.size());
2080 
2081  AssertThrow(false, ExcNotImplemented());
2082 }
2083 
2084 
2085 
2086 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2087 Point<spacedim>
2088 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2089  transform_unit_to_real_cell(
2090  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2091  const Point<dim> & p) const
2092 {
2093  // Use the get_data function to create an InternalData with data vectors of
2094  // the right size and transformation shape values already computed at point
2095  // p.
2096  const Quadrature<dim> point_quadrature(p);
2097  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2098  get_data(update_quadrature_points | update_jacobians, point_quadrature));
2099  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2100  ExcInternalError());
2101 
2102  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2103 
2104  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2105 }
2106 
2107 
2108 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2109 Point<spacedim>
2110 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2111  do_transform_unit_to_real_cell(const InternalData &data) const
2112 {
2113  Point<spacedim> p_real;
2114 
2115  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2116  {
2117  unsigned int comp_i =
2118  euler_dof_handler->get_fe().system_to_component_index(i).first;
2119  if (fe_mask[comp_i])
2120  p_real[fe_to_real[comp_i]] +=
2121  data.local_dof_values[i] * data.shape(0, i);
2122  }
2123 
2124  return p_real;
2125 }
2126 
2127 
2128 
2129 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2130 Point<dim>
2131 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2132  transform_real_to_unit_cell(
2133  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2134  const Point<spacedim> & p) const
2135 {
2136  // first a Newton iteration based on the real mapping. It uses the center
2137  // point of the cell as a starting point
2138  Point<dim> initial_p_unit;
2139  try
2140  {
2141  initial_p_unit =
2142  StaticMappingQ1<dim, spacedim>::mapping.transform_real_to_unit_cell(
2143  cell, p);
2144  }
2145  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2146  {
2147  // mirror the conditions of the code below to determine if we need to
2148  // use an arbitrary starting point or if we just need to rethrow the
2149  // exception
2150  for (unsigned int d = 0; d < dim; ++d)
2151  initial_p_unit[d] = 0.5;
2152  }
2153 
2154  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2155 
2156  // for (unsigned int d=0; d<dim; ++d)
2157  // initial_p_unit[d] = 0.;
2158 
2159  const Quadrature<dim> point_quadrature(initial_p_unit);
2160 
2161  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
2162  if (spacedim > dim)
2163  update_flags |= update_jacobian_grads;
2164  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2165  get_data(update_flags, point_quadrature));
2166  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2167  ExcInternalError());
2168 
2169  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2170 
2171  return do_transform_real_to_unit_cell(cell,
2172  p,
2173  initial_p_unit,
2174  dynamic_cast<InternalData &>(*mdata));
2175 }
2176 
2177 
2178 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2179 Point<dim>
2180 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2181  do_transform_real_to_unit_cell(
2182  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2183  const Point<spacedim> & p,
2184  const Point<dim> & initial_p_unit,
2185  InternalData & mdata) const
2186 {
2187  const unsigned int n_shapes = mdata.shape_values.size();
2188  (void)n_shapes;
2189  Assert(n_shapes != 0, ExcInternalError());
2190  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2191 
2192 
2193  // Newton iteration to solve
2194  // f(x)=p(x)-p=0
2195  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2196  // The start value was set to be the
2197  // linear approximation to the cell
2198  // The shape values and derivatives
2199  // of the mapping at this point are
2200  // previously computed.
2201  // f(x)
2202  Point<dim> p_unit = initial_p_unit;
2203  Point<dim> f;
2204  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2205  Point<spacedim> p_real(do_transform_unit_to_real_cell(mdata));
2206  Tensor<1, spacedim> p_minus_F = p - p_real;
2207  const double eps = 1.e-12 * cell->diameter();
2208  const unsigned int newton_iteration_limit = 20;
2209  unsigned int newton_iteration = 0;
2210  while (p_minus_F.norm_square() > eps * eps)
2211  {
2212  // f'(x)
2213  Point<spacedim> DF[dim];
2214  Tensor<2, dim> df;
2215  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2216  {
2217  const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2218  const unsigned int comp_k =
2219  euler_dof_handler->get_fe().system_to_component_index(k).first;
2220  if (fe_mask[comp_k])
2221  for (unsigned int j = 0; j < dim; ++j)
2222  DF[j][fe_to_real[comp_k]] +=
2223  mdata.local_dof_values[k] * grad_k[j];
2224  }
2225  for (unsigned int j = 0; j < dim; ++j)
2226  {
2227  f[j] = DF[j] * p_minus_F;
2228  for (unsigned int l = 0; l < dim; ++l)
2229  df[j][l] = -DF[j] * DF[l];
2230  }
2231  // Solve [f'(x)]d=f(x)
2232  const Tensor<1, dim> delta =
2233  invert(df) * static_cast<const Tensor<1, dim> &>(f);
2234  // do a line search
2235  double step_length = 1;
2236  do
2237  {
2238  // update of p_unit. The
2239  // spacedimth component of
2240  // transformed point is simply
2241  // ignored in codimension one
2242  // case. When this component is
2243  // not zero, then we are
2244  // projecting the point to the
2245  // surface or curve identified
2246  // by the cell.
2247  Point<dim> p_unit_trial = p_unit;
2248  for (unsigned int i = 0; i < dim; ++i)
2249  p_unit_trial[i] -= step_length * delta[i];
2250  // shape values and derivatives
2251  // at new p_unit point
2252  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2253  mdata);
2254  // f(x)
2255  Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2256  const Tensor<1, spacedim> f_trial = p - p_real_trial;
2257  // see if we are making progress with the current step length
2258  // and if not, reduce it by a factor of two and try again
2259  if (f_trial.norm() < p_minus_F.norm())
2260  {
2261  p_real = p_real_trial;
2262  p_unit = p_unit_trial;
2263  p_minus_F = f_trial;
2264  break;
2265  }
2266  else if (step_length > 0.05)
2267  step_length /= 2;
2268  else
2269  goto failure;
2270  }
2271  while (true);
2272  ++newton_iteration;
2273  if (newton_iteration > newton_iteration_limit)
2274  goto failure;
2275  }
2276  return p_unit;
2277  // if we get to the following label, then we have either run out
2278  // of Newton iterations, or the line search has not converged.
2279  // in either case, we need to give up, so throw an exception that
2280  // can then be caught
2281 failure:
2282  AssertThrow(false,
2283  (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
2284  // ...the compiler wants us to return something, though we can
2285  // of course never get here...
2286  return Point<dim>();
2287 }
2288 
2289 
2290 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2291 unsigned int
2292 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_degree() const
2293 {
2294  return euler_dof_handler->get_fe().degree;
2295 }
2296 
2297 
2298 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2299 ComponentMask
2300 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::get_component_mask()
2301  const
2302 {
2303  return this->fe_mask;
2304 }
2305 
2306 
2307 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2308 std::unique_ptr<Mapping<dim, spacedim>>
2309 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::clone() const
2310 {
2311  return std_cxx14::make_unique<
2312  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>>(*this);
2313 }
2314 
2315 
2316 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2317 void
2318 MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::update_internal_dofs(
2319  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2320  const typename MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2321  InternalData &data) const
2322 {
2323  Assert(euler_dof_handler != nullptr,
2324  ExcMessage("euler_dof_handler is empty"));
2325 
2326  typename DoFHandlerType::cell_iterator dof_cell(*cell, euler_dof_handler);
2327  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2328  if (uses_level_dofs)
2329  {
2330  AssertIndexRange(cell->level(), euler_vector.size());
2331  AssertDimension(euler_vector[cell->level()]->size(),
2332  euler_dof_handler->n_dofs(cell->level()));
2333  }
2334  else
2335  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2336 
2337  if (uses_level_dofs)
2338  dof_cell->get_mg_dof_indices(data.local_dof_indices);
2339  else
2340  dof_cell->get_dof_indices(data.local_dof_indices);
2341 
2342  const VectorType &vector =
2343  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2344 
2345  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2346  data.local_dof_values[i] =
2347  internal::ElementAccess<VectorType>::get(vector,
2348  data.local_dof_indices[i]);
2349 }
2350 
2351 // explicit instantiations
2352 #define SPLIT_INSTANTIATIONS_COUNT 2
2353 #ifndef SPLIT_INSTANTIATIONS_INDEX
2354 # define SPLIT_INSTANTIATIONS_INDEX 0
2355 #endif
2356 #include "mapping_fe_field.inst"
2357 
2358 
2359 DEAL_II_NAMESPACE_CLOSE
update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
la_vector.h
MGLevelObject::max_level
unsigned int max_level() const
Definition: mg_level_object.h:246
update_values
Shape function values.
Definition: fe_update_flags.h:78
full_matrix.h
MappingFEField::get_data
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_fe_field.cc:608
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2645
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
dof_accessor.h
array_view.h
update_jacobian_2nd_derivatives
Definition: fe_update_flags.h:201
apply_transformation
Tensor< 1, spacedim, Number > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number > &grad_F, const Tensor< 1, dim, Number > &d_x)
Definition: derivative_form.h:399
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
petsc_block_vector.h
MappingFEField::compute_face_data
void compute_face_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
Definition: mapping_fe_field.cc:567
update_contravariant_transformation
Contravariant transformation.
Definition: fe_update_flags.h:175
update_jacobian_3rd_derivatives
Definition: fe_update_flags.h:210
MappingFEField::fe_values_mutex
std::mutex fe_values_mutex
Definition: mapping_fe_field.h:637
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
trilinos_tpetra_vector.h
MappingFEField::clone
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_fe_field.cc:2309
fe_values.h
Quadrature::get_weights
const std::vector< double > & get_weights() const
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
block_vector.h
MappingFEField::compute_shapes_virtual
virtual void compute_shapes_virtual(const std::vector< Point< dim >> &unit_points, typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
Definition: mapping_fe_field.cc:423
mapping_covariant_gradient
Definition: mapping.h:84
internal::FEValuesImplementation::MappingRelatedData::boundary_forms
std::vector< Tensor< 1, spacedim > > boundary_forms
Definition: fe_update_flags.h:511
update_jacobians
Volume element.
Definition: fe_update_flags.h:150
trilinos_parallel_block_vector.h
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
internal::FEValuesImplementation::MappingRelatedData::jacobian_2nd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Definition: fe_update_flags.h:476
update_boundary_forms
Outer normal vector, not normalized.
Definition: fe_update_flags.h:103
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=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:442
MappingFEField::fe_to_real
std::vector< unsigned int > fe_to_real
Definition: mapping_fe_field.h:626
update_volume_elements
Determinant of the Jacobian.
Definition: fe_update_flags.h:192
internal::FEValuesImplementation::MappingRelatedData::jacobian_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Definition: fe_update_flags.h:488
update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
MappingFEField::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_field.cc:1537
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
MappingFEField::fe_mask
const ComponentMask fe_mask
Definition: mapping_fe_field.h:615
second
Point< 2 > second
Definition: grid_out.cc:4353
MappingKind
MappingKind
Definition: mapping.h:62
memory_consumption.h
QProjector::DataSetDescriptor::cell
static DataSetDescriptor cell()
Definition: qprojector.h:346
quadrature_lib.h
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition: fe_update_flags.h:482
MappingFEField::do_transform_unit_to_real_cell
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
Definition: mapping_fe_field.cc:2111
MGLevelObject::min_level
unsigned int min_level() const
Definition: mg_level_object.h:238
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:2407
MappingFEField::get_degree
unsigned int get_degree() const
Definition: mapping_fe_field.cc:2292
MappingFEField::fill_fe_face_values
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_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_field.cc:1744
update_jacobian_pushed_forward_3rd_derivatives
Definition: fe_update_flags.h:215
update_jacobian_pushed_forward_grads
Definition: fe_update_flags.h:197
MappingFEField::fe_values
FEValues< dim, spacedim > fe_values
Definition: mapping_fe_field.h:632
petsc_vector.h
MappingFEField::InternalData::shape_third_derivatives
std::vector< Tensor< 3, dim > > shape_third_derivatives
Definition: mapping_fe_field.h:391
tensor_product_polynomials.h
MappingFEField::InternalData::shape_derivatives
std::vector< Tensor< 1, dim > > shape_derivatives
Definition: mapping_fe_field.h:375
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
MappingFEField::InternalData::shape
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: mapping_fe_field.cc:340
fe_system.h
MappingFEField::InternalData::shape_second_derivatives
std::vector< Tensor< 2, dim > > shape_second_derivatives
Definition: mapping_fe_field.h:383
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
MappingFEField::InternalData::shape_fourth_derivatives
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
Definition: mapping_fe_field.h:399
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Mapping
Abstract base class for mapping classes.
Definition: mapping.h:302
la_parallel_block_vector.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:607
MappingFEField::InternalData::volume_elements
std::vector< double > volume_elements
Definition: mapping_fe_field.h:464
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
trilinos_epetra_vector.h
ComponentMask::size
unsigned int size() const
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
update_jacobian_pushed_forward_2nd_derivatives
Definition: fe_update_flags.h:206
mapping_contravariant
Definition: mapping.h:78
update_inverse_jacobians
Volume element.
Definition: fe_update_flags.h:161
Tensor::norm_square
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
MappingFEField::get_face_data
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition: mapping_fe_field.cc:624
MappingFEField::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_field.cc:2132
MappingFEField::InternalData::covariant
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Definition: mapping_fe_field.h:449
mapping_covariant
Definition: mapping.h:73
update_jacobian_grads
Gradient of volume element.
Definition: fe_update_flags.h:155
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
MappingFEField::InternalData::local_dof_values
std::vector< double > local_dof_values
Definition: mapping_fe_field.h:479
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
MappingFEField::InternalData::third_derivative
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: mapping_fe_field.cc:147
fe_q.h
MappingFEField::get_vertices
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_fe_field.cc:362
Quadrature::size
unsigned int size() const
MappingFEField::InternalData::contravariant
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition: mapping_fe_field.h:458
MappingFEField::InternalData::shape_values
std::vector< double > shape_values
Definition: mapping_fe_field.h:368
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_grads
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Definition: fe_update_flags.h:470
MappingFEField::InternalData::aux
std::vector< std::vector< Tensor< 1, spacedim > > > aux
Definition: mapping_fe_field.h:469
MappingFEField::MappingFEField
friend class MappingFEField
Definition: mapping_fe_field.h:654
internal::FEValuesImplementation::MappingRelatedData::quadrature_points
std::vector< Point< spacedim > > quadrature_points
Definition: fe_update_flags.h:501
GeometryInfo::project_to_unit_cell
static Point< dim > project_to_unit_cell(const Point< dim > &p)
vector_tools.h
internal::ElementAccess::get
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
Definition: vector_element_access.h:73
MappingFEField::compute_data
void compute_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
Definition: mapping_fe_field.cc:516
tria_iterator.h
FiniteElement::system_to_component_index
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3082
MappingFEField::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_field.cc:1783
MappingFEField::InternalData::second_derivative
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: mapping_fe_field.cc:123
MappingFEField::uses_level_dofs
const bool uses_level_dofs
Definition: mapping_fe_field.h:536
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
mapping_fe_field.h
update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
MappingFEField::InternalData::fourth_derivative
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: mapping_fe_field.cc:170
internal::FEValuesImplementation::MappingRelatedData::jacobians
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
Definition: fe_update_flags.h:453
la_parallel_vector.h
MappingFEField::euler_dof_handler
SmartPointer< const DoFHandlerType, MappingFEField< dim, spacedim, VectorType, DoFHandlerType > > euler_dof_handler
Definition: mapping_fe_field.h:551
MappingFEField::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_field.cc:2089
fe_tools.h
MappingFEField::euler_vector
std::vector< SmartPointer< const VectorType, MappingFEField< dim, spacedim, VectorType, DoFHandlerType > > > euler_vector
Definition: mapping_fe_field.h:544
QProjector::DataSetDescriptor::subface
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
MappingFEField::InternalData::derivative
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: mapping_fe_field.cc:100
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
MappingFEField::get_component_mask
ComponentMask get_component_mask() const
Definition: mapping_fe_field.cc:2300
MappingFEField::update_internal_dofs
void update_internal_dofs(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
Definition: mapping_fe_field.cc:2318
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
MappingFEField::preserves_vertex_locations
virtual bool preserves_vertex_locations() const override
Definition: mapping_fe_field.cc:353
invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3467
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_3rd_derivatives
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
Definition: fe_update_flags.h:494
MGLevelObject::resize
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
Definition: mg_level_object.h:199
internal::FEValuesImplementation::MappingRelatedData::jacobian_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
Definition: fe_update_flags.h:459
mapping_piola
Definition: mapping.h:98
MappingFEField::ExcInactiveCell
static ::ExceptionBase & ExcInactiveCell()
internal::FEValuesImplementation::MappingRelatedData::inverse_jacobians
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
Definition: fe_update_flags.h:464
cross_product_2d
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2381
value
static const bool value
Definition: dof_tools_constraints.cc:433
MappingFEField::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_field.cc:1968
MappingFEField::InternalData::unit_tangentials
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
Definition: mapping_fe_field.h:416
mapping_q1.h
MappingFEField::do_transform_real_to_unit_cell
Point< dim > do_transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit, InternalData &mdata) const
Definition: mapping_fe_field.cc:2181
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
MappingFEField::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_field.cc:640
trilinos_vector.h
update_covariant_transformation
Covariant transformation.
Definition: fe_update_flags.h:168
internal::FEValuesImplementation::MappingRelatedData::normal_vectors
std::vector< Tensor< 1, spacedim > > normal_vectors
Definition: fe_update_flags.h:506
MappingFEField::InternalData::memory_consumption
virtual std::size_t memory_consumption() const override
Definition: mapping_fe_field.cc:78
MappingFEField::InternalData::InternalData
InternalData(const FiniteElement< dim, spacedim > &fe, const ComponentMask &mask)
Definition: mapping_fe_field.cc:64
MappingFEField::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_fe_field.cc:462
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Mapping::InternalDataBase::update_each
UpdateFlags update_each
Definition: mapping.h:630
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
memory.h
QProjector::DataSetDescriptor::face
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1141
Generated by   doxygen 1.8.13