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