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