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