]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some code to a .templates.h file.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 26 May 2025 23:05:17 +0000 (17:05 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 26 May 2025 23:08:43 +0000 (17:08 -0600)
include/deal.II/fe/mapping_fe_field.templates.h [new file with mode: 0644]
source/fe/mapping_fe_field.cc

diff --git a/include/deal.II/fe/mapping_fe_field.templates.h b/include/deal.II/fe/mapping_fe_field.templates.h
new file mode 100644 (file)
index 0000000..ea7b935
--- /dev/null
@@ -0,0 +1,2464 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2015 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+#ifndef dealii_mapping_fe_field_templates_h
+#define dealii_mapping_fe_field_templates_h
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <memory>
+#include <numeric>
+
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim, int spacedim, typename VectorType>
+MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
+  const FiniteElement<dim, spacedim> &fe,
+  const ComponentMask                &mask)
+  : fe(&fe)
+  , unit_tangentials()
+  , n_shape_functions(fe.n_dofs_per_cell())
+  , mask(mask)
+  , local_dof_indices(fe.n_dofs_per_cell())
+  , local_dof_values(fe.n_dofs_per_cell())
+{}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::InternalData::reinit(
+  const UpdateFlags      update_flags,
+  const Quadrature<dim> &quadrature)
+{
+  // store the flags in the internal data object so we can access them
+  // in fill_fe_*_values(). use the transitive hull of the required
+  // flags
+  this->update_each = update_flags;
+
+  const unsigned int             n_q_points = quadrature.size();
+  const std::vector<Point<dim>> &points     = quadrature.get_points();
+
+  // see if we need the (transformation) shape function values
+  // and/or gradients and resize the necessary arrays
+  if (update_flags & update_quadrature_points)
+    {
+      shape_values.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          shape(point, i) = fe->shape_value(i, points[point]);
+    }
+
+  if (update_flags &
+      (update_covariant_transformation | update_contravariant_transformation |
+       update_JxW_values | update_boundary_forms | update_normal_vectors |
+       update_jacobians | update_jacobian_grads | update_inverse_jacobians))
+    {
+      shape_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          derivative(point, i) = fe->shape_grad(i, points[point]);
+    }
+
+  if (update_flags & update_covariant_transformation)
+    covariant.resize(n_q_points);
+
+  if (update_flags & update_contravariant_transformation)
+    contravariant.resize(n_q_points);
+
+  if (update_flags & update_volume_elements)
+    volume_elements.resize(n_q_points);
+
+  if (update_flags &
+      (update_jacobian_grads | update_jacobian_pushed_forward_grads))
+    {
+      shape_second_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          second_derivative(point, i) = fe->shape_grad_grad(i, points[point]);
+    }
+
+  if (update_flags & (update_jacobian_2nd_derivatives |
+                      update_jacobian_pushed_forward_2nd_derivatives))
+    {
+      shape_third_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          third_derivative(point, i) =
+            fe->shape_3rd_derivative(i, points[point]);
+    }
+
+  if (update_flags & (update_jacobian_3rd_derivatives |
+                      update_jacobian_pushed_forward_3rd_derivatives))
+    {
+      shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          fourth_derivative(point, i) =
+            fe->shape_4th_derivative(i, points[point]);
+    }
+
+  // This (for face values and simplices) can be different for different
+  // calls, so always copy
+  quadrature_weights = quadrature.get_weights();
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+std::size_t
+MappingFEField<dim, spacedim, VectorType>::InternalData::memory_consumption()
+  const
+{
+  DEAL_II_NOT_IMPLEMENTED();
+  return 0;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+double &
+MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
+  const unsigned int qpoint,
+  const unsigned int shape_nr)
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
+  return shape_values[qpoint * n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+const Tensor<1, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr) const
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_derivatives.size());
+  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Tensor<1, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr)
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_derivatives.size());
+  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+const Tensor<2, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr) const
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_second_derivatives.size());
+  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Tensor<2, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr)
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_second_derivatives.size());
+  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+const Tensor<3, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr) const
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_third_derivatives.size());
+  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Tensor<3, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr)
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_third_derivatives.size());
+  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+const Tensor<4, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr) const
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_fourth_derivatives.size());
+  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Tensor<4, dim> &
+MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
+  const unsigned int qpoint,
+  const unsigned int shape_nr)
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
+                   shape_fourth_derivatives.size());
+  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+MappingFEField<dim, spacedim, VectorType>::MappingFEField(
+  const DoFHandler<dim, spacedim> &euler_dof_handler,
+  const VectorType                &euler_vector,
+  const ComponentMask             &mask)
+  : reference_cell(euler_dof_handler.get_fe().reference_cell())
+  , uses_level_dofs(false)
+  , euler_vector({&euler_vector})
+  , euler_dof_handler(&euler_dof_handler)
+  , fe_mask(mask.size() != 0u ?
+              mask :
+              ComponentMask(
+                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
+                true))
+  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
+  , fe_values(this->euler_dof_handler->get_fe(),
+              reference_cell.template get_nodal_type_quadrature<dim>(),
+              update_values)
+{
+  AssertDimension(euler_dof_handler.n_dofs(), euler_vector.size());
+  unsigned int size = 0;
+  for (unsigned int i = 0; i < fe_mask.size(); ++i)
+    {
+      if (fe_mask[i])
+        fe_to_real[i] = size++;
+    }
+  AssertDimension(size, spacedim);
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+MappingFEField<dim, spacedim, VectorType>::MappingFEField(
+  const DoFHandler<dim, spacedim> &euler_dof_handler,
+  const std::vector<VectorType>   &euler_vector,
+  const ComponentMask             &mask)
+  : reference_cell(euler_dof_handler.get_fe().reference_cell())
+  , uses_level_dofs(true)
+  , euler_dof_handler(&euler_dof_handler)
+  , fe_mask(mask.size() != 0u ?
+              mask :
+              ComponentMask(
+                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
+                true))
+  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
+  , fe_values(this->euler_dof_handler->get_fe(),
+              reference_cell.template get_nodal_type_quadrature<dim>(),
+              update_values)
+{
+  unsigned int size = 0;
+  for (unsigned int i = 0; i < fe_mask.size(); ++i)
+    {
+      if (fe_mask[i])
+        fe_to_real[i] = size++;
+    }
+  AssertDimension(size, spacedim);
+
+  Assert(euler_dof_handler.has_level_dofs(),
+         ExcMessage("The underlying DoFHandler object did not call "
+                    "distribute_mg_dofs(). In this case, the construction via "
+                    "level vectors does not make sense."));
+  AssertDimension(euler_vector.size(),
+                  euler_dof_handler.get_triangulation().n_global_levels());
+  this->euler_vector.clear();
+  this->euler_vector.resize(euler_vector.size());
+  for (unsigned int i = 0; i < euler_vector.size(); ++i)
+    {
+      AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
+      this->euler_vector[i] = &euler_vector[i];
+    }
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+MappingFEField<dim, spacedim, VectorType>::MappingFEField(
+  const DoFHandler<dim, spacedim> &euler_dof_handler,
+  const MGLevelObject<VectorType> &euler_vector,
+  const ComponentMask             &mask)
+  : reference_cell(euler_dof_handler.get_fe().reference_cell())
+  , uses_level_dofs(true)
+  , euler_dof_handler(&euler_dof_handler)
+  , fe_mask(mask.size() != 0u ?
+              mask :
+              ComponentMask(
+                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
+                true))
+  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
+  , fe_values(this->euler_dof_handler->get_fe(),
+              reference_cell.template get_nodal_type_quadrature<dim>(),
+              update_values)
+{
+  unsigned int size = 0;
+  for (unsigned int i = 0; i < fe_mask.size(); ++i)
+    {
+      if (fe_mask[i])
+        fe_to_real[i] = size++;
+    }
+  AssertDimension(size, spacedim);
+
+  Assert(euler_dof_handler.has_level_dofs(),
+         ExcMessage("The underlying DoFHandler object did not call "
+                    "distribute_mg_dofs(). In this case, the construction via "
+                    "level vectors does not make sense."));
+  AssertDimension(euler_vector.max_level() + 1,
+                  euler_dof_handler.get_triangulation().n_global_levels());
+  this->euler_vector.clear();
+  this->euler_vector.resize(
+    euler_dof_handler.get_triangulation().n_global_levels());
+  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
+       ++i)
+    {
+      AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
+      this->euler_vector[i] = &euler_vector[i];
+    }
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+MappingFEField<dim, spacedim, VectorType>::MappingFEField(
+  const MappingFEField<dim, spacedim, VectorType> &mapping)
+  : reference_cell(mapping.reference_cell)
+  , uses_level_dofs(mapping.uses_level_dofs)
+  , euler_vector(mapping.euler_vector)
+  , euler_dof_handler(mapping.euler_dof_handler)
+  , fe_mask(mapping.fe_mask)
+  , fe_to_real(mapping.fe_to_real)
+  , fe_values(mapping.euler_dof_handler->get_fe(),
+              reference_cell.template get_nodal_type_quadrature<dim>(),
+              update_values)
+{}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+inline const double &
+MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
+  const unsigned int qpoint,
+  const unsigned int shape_nr) const
+{
+  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
+  return shape_values[qpoint * n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+bool
+MappingFEField<dim, spacedim, VectorType>::preserves_vertex_locations() const
+{
+  return false;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+bool
+MappingFEField<dim, spacedim, VectorType>::is_compatible_with(
+  const ReferenceCell &reference_cell) const
+{
+  Assert(dim == reference_cell.get_dimension(),
+         ExcMessage("The dimension of your mapping (" +
+                    Utilities::to_string(dim) +
+                    ") and the reference cell cell_type (" +
+                    Utilities::to_string(reference_cell.get_dimension()) +
+                    " ) do not agree."));
+
+  return this->reference_cell == reference_cell;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+boost::container::small_vector<Point<spacedim>,
+#ifndef _MSC_VER
+                               ReferenceCells::max_n_vertices<dim>()
+#else
+                               GeometryInfo<dim>::vertices_per_cell
+#endif
+                               >
+MappingFEField<dim, spacedim, VectorType>::get_vertices(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
+{
+  // we transform our tria iterator into a dof iterator so we can access
+  // data not associated with triangulations
+  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
+    *cell, euler_dof_handler);
+
+  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
+  AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
+  AssertDimension(fe_to_real.size(),
+                  euler_dof_handler->get_fe().n_components());
+  if (uses_level_dofs)
+    {
+      AssertIndexRange(cell->level(), euler_vector.size());
+      AssertDimension(euler_vector[cell->level()]->size(),
+                      euler_dof_handler->n_dofs(cell->level()));
+    }
+  else
+    AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
+
+  {
+    std::lock_guard<std::mutex> lock(fe_values_mutex);
+    fe_values.reinit(dof_cell);
+  }
+  const unsigned int dofs_per_cell =
+    euler_dof_handler->get_fe().n_dofs_per_cell();
+  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+  if (uses_level_dofs)
+    dof_cell->get_mg_dof_indices(dof_indices);
+  else
+    dof_cell->get_dof_indices(dof_indices);
+
+  const VectorType &vector =
+    uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
+
+  boost::container::small_vector<Point<spacedim>,
+#ifndef _MSC_VER
+                                 ReferenceCells::max_n_vertices<dim>()
+#else
+                                 GeometryInfo<dim>::vertices_per_cell
+#endif
+                                 >
+    vertices(cell->n_vertices());
+  for (unsigned int i = 0; i < dofs_per_cell; ++i)
+    {
+      const unsigned int comp = fe_to_real
+        [euler_dof_handler->get_fe().system_to_component_index(i).first];
+      if (comp != numbers::invalid_unsigned_int)
+        {
+          typename VectorType::value_type value =
+            internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
+          if (euler_dof_handler->get_fe().is_primitive(i))
+            for (const unsigned int v : cell->vertex_indices())
+              vertices[v][comp] += fe_values.shape_value(i, v) * value;
+          else
+            DEAL_II_NOT_IMPLEMENTED();
+        }
+    }
+
+  return vertices;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+UpdateFlags
+MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
+  const UpdateFlags in) const
+{
+  // add flags if the respective quantities are necessary to compute
+  // what we need. note that some flags appear in both conditions and
+  // in subsequent set operations. this leads to some circular
+  // logic. the only way to treat this is to iterate. since there are
+  // 5 if-clauses in the loop, it will take at most 4 iterations to
+  // converge. do them:
+  UpdateFlags out = in;
+  for (unsigned int i = 0; i < 5; ++i)
+    {
+      // The following is a little incorrect:
+      // If not applied on a face,
+      // update_boundary_forms does not
+      // make sense. On the other hand,
+      // it is necessary on a
+      // face. Currently,
+      // update_boundary_forms is simply
+      // ignored for the interior of a
+      // cell.
+      if (out & (update_JxW_values | update_normal_vectors))
+        out |= update_boundary_forms;
+
+      if (out &
+          (update_covariant_transformation | update_jacobian_grads |
+           update_jacobians | update_boundary_forms | update_normal_vectors))
+        out |= update_contravariant_transformation;
+
+      if (out &
+          (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
+           update_jacobian_pushed_forward_2nd_derivatives |
+           update_jacobian_pushed_forward_3rd_derivatives))
+        out |= update_covariant_transformation;
+
+      // The contravariant transformation is used in the Piola
+      // transformation, which requires the determinant of the Jacobi
+      // matrix of the transformation.  Because we have no way of
+      // knowing here whether the finite element wants to use the
+      // contravariant or the Piola transforms, we add the volume elements
+      // to the list of flags to be updated for each cell.
+      if (out & update_contravariant_transformation)
+        out |= update_volume_elements;
+
+      if (out & update_normal_vectors)
+        out |= update_volume_elements;
+    }
+
+  return out;
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::compute_face_data(
+  const unsigned int n_original_q_points,
+  InternalData      &data) const
+{
+  // Set to the size of a single quadrature object for faces, as the size set
+  // in in reinit() is for all points
+  if (data.update_each & update_covariant_transformation)
+    data.covariant.resize(n_original_q_points);
+
+  if (data.update_each & update_contravariant_transformation)
+    data.contravariant.resize(n_original_q_points);
+
+  if (data.update_each & update_volume_elements)
+    data.volume_elements.resize(n_original_q_points);
+
+  if (dim > 1)
+    {
+      if (data.update_each & update_boundary_forms)
+        {
+          data.aux.resize(
+            dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
+
+
+          // TODO: only a single reference cell type possible...
+          const auto n_faces = reference_cell.n_faces();
+
+          // Compute tangentials to the unit cell.
+          for (unsigned int i = 0; i < n_faces; ++i)
+            {
+              data.unit_tangentials[i].resize(n_original_q_points);
+              std::fill(data.unit_tangentials[i].begin(),
+                        data.unit_tangentials[i].end(),
+                        reference_cell.template face_tangent_vector<dim>(i, 0));
+              if (dim > 2)
+                {
+                  data.unit_tangentials[n_faces + i].resize(
+                    n_original_q_points);
+                  std::fill(
+                    data.unit_tangentials[n_faces + i].begin(),
+                    data.unit_tangentials[n_faces + i].end(),
+                    reference_cell.template face_tangent_vector<dim>(i, 1));
+                }
+            }
+        }
+    }
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFEField<dim, spacedim, VectorType>::get_data(
+  const UpdateFlags      update_flags,
+  const Quadrature<dim> &quadrature) const
+{
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
+    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
+  data_ptr->reinit(requires_update_flags(update_flags), quadrature);
+
+  return data_ptr;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFEField<dim, spacedim, VectorType>::get_face_data(
+  const UpdateFlags               update_flags,
+  const hp::QCollection<dim - 1> &quadrature) const
+{
+  AssertDimension(quadrature.size(), 1);
+
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
+    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
+  const Quadrature<dim> q(
+    QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
+  data.reinit(requires_update_flags(update_flags), q);
+  this->compute_face_data(quadrature[0].size(), data);
+
+  return data_ptr;
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFEField<dim, spacedim, VectorType>::get_subface_data(
+  const UpdateFlags          update_flags,
+  const Quadrature<dim - 1> &quadrature) const
+{
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
+    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
+  const Quadrature<dim> q(
+    QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
+  data.reinit(requires_update_flags(update_flags), q);
+  this->compute_face_data(quadrature.size(), data);
+
+  return data_ptr;
+}
+
+
+
+namespace internal
+{
+  namespace MappingFEFieldImplementation
+  {
+    namespace
+    {
+      /**
+       * Compute the locations of quadrature points on the object described by
+       * the first argument (and the cell for which the mapping support points
+       * have already been set), but only if the update_flags of the @p data
+       * argument indicate so.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_compute_q_points(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real,
+        std::vector<Point<spacedim>>       &quadrature_points)
+      {
+        const UpdateFlags update_flags = data.update_each;
+
+        if (update_flags & update_quadrature_points)
+          {
+            for (unsigned int point = 0; point < quadrature_points.size();
+                 ++point)
+              {
+                Point<spacedim> result;
+                const double   *shape = &data.shape(point + data_set, 0);
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      result[fe_to_real[comp_k]] +=
+                        data.local_dof_values[k] * shape[k];
+                  }
+
+                quadrature_points[point] = result;
+              }
+          }
+      }
+
+      /**
+       * Update the co- and contravariant matrices as well as their determinant,
+       * for the cell described stored in the data object, but only if the
+       * update_flags of the @p data argument indicate so.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_Jacobians(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real)
+      {
+        const UpdateFlags update_flags = data.update_each;
+
+        // then Jacobians
+        if (update_flags & update_contravariant_transformation)
+          {
+            const unsigned int n_q_points = data.contravariant.size();
+
+            Assert(data.n_shape_functions > 0, ExcInternalError());
+
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<1, dim> *data_derv =
+                  &data.derivative(point + data_set, 0);
+
+                Tensor<1, dim> result[spacedim];
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      result[fe_to_real[comp_k]] +=
+                        data.local_dof_values[k] * data_derv[k];
+                  }
+
+                // write result into contravariant data
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  {
+                    data.contravariant[point][i] = result[i];
+                  }
+              }
+          }
+
+        if (update_flags & update_covariant_transformation)
+          {
+            AssertDimension(data.covariant.size(), data.contravariant.size());
+            for (unsigned int point = 0; point < data.contravariant.size();
+                 ++point)
+              data.covariant[point] =
+                (data.contravariant[point]).covariant_form();
+          }
+
+        if (update_flags & update_volume_elements)
+          {
+            AssertDimension(data.contravariant.size(),
+                            data.volume_elements.size());
+            for (unsigned int point = 0; point < data.contravariant.size();
+                 ++point)
+              data.volume_elements[point] =
+                data.contravariant[point].determinant();
+          }
+      }
+
+      /**
+       * Update the Hessian of the transformation from unit to real cell, the
+       * Jacobian gradients.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_grads(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                                &data,
+        const FiniteElement<dim, spacedim>            &fe,
+        const ComponentMask                           &fe_mask,
+        const std::vector<unsigned int>               &fe_to_real,
+        std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_grads)
+          {
+            const unsigned int n_q_points = jacobian_grads.size();
+
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<2, dim> *second =
+                  &data.second_derivative(point + data_set, 0);
+
+                DerivativeForm<2, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          result[fe_to_real[comp_k]][j][l] +=
+                            (second[k][j][l] * data.local_dof_values[k]);
+                  }
+
+                // never touch any data for j=dim in case dim<spacedim, so
+                // it will always be zero as it was initialized
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < dim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      jacobian_grads[point][i][j][l] = result[i][j][l];
+              }
+          }
+      }
+
+      /**
+       * Update the Hessian of the transformation from unit to real cell, the
+       * Jacobian gradients, pushed forward to the real cell coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_pushed_forward_grads(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real,
+        std::vector<Tensor<3, spacedim>>   &jacobian_pushed_forward_grads)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_grads)
+          {
+            const unsigned int n_q_points =
+              jacobian_pushed_forward_grads.size();
+
+            double tmp[spacedim][spacedim][spacedim];
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<2, dim> *second =
+                  &data.second_derivative(point + data_set, 0);
+
+                DerivativeForm<2, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          result[fe_to_real[comp_k]][j][l] +=
+                            (second[k][j][l] * data.local_dof_values[k]);
+                  }
+
+                // first push forward the j-components
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      {
+                        tmp[i][j][l] =
+                          result[i][0][l] * data.covariant[point][j][0];
+                        for (unsigned int jr = 1; jr < dim; ++jr)
+                          {
+                            tmp[i][j][l] +=
+                              result[i][jr][l] * data.covariant[point][j][jr];
+                          }
+                      }
+
+                // now, pushing forward the l-components
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      {
+                        jacobian_pushed_forward_grads[point][i][j][l] =
+                          tmp[i][j][0] * data.covariant[point][l][0];
+                        for (unsigned int lr = 1; lr < dim; ++lr)
+                          {
+                            jacobian_pushed_forward_grads[point][i][j][l] +=
+                              tmp[i][j][lr] * data.covariant[point][l][lr];
+                          }
+                      }
+              }
+          }
+      }
+
+      /**
+       * Update the third derivative of the transformation from unit to real
+       * cell, the Jacobian hessians.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_2nd_derivatives(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                                &data,
+        const FiniteElement<dim, spacedim>            &fe,
+        const ComponentMask                           &fe_mask,
+        const std::vector<unsigned int>               &fe_to_real,
+        std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_2nd_derivatives)
+          {
+            const unsigned int n_q_points = jacobian_2nd_derivatives.size();
+
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<3, dim> *third =
+                  &data.third_derivative(point + data_set, 0);
+
+                DerivativeForm<3, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          for (unsigned int m = 0; m < dim; ++m)
+                            result[fe_to_real[comp_k]][j][l][m] +=
+                              (third[k][j][l][m] * data.local_dof_values[k]);
+                  }
+
+                // never touch any data for j=dim in case dim<spacedim, so
+                // it will always be zero as it was initialized
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < dim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        jacobian_2nd_derivatives[point][i][j][l][m] =
+                          result[i][j][l][m];
+              }
+          }
+      }
+
+      /**
+       * Update the third derivative of the transformation from unit to real
+       * cell, the Jacobian hessians, pushed forward to the real cell
+       * coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_pushed_forward_2nd_derivatives(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real,
+        std::vector<Tensor<4, spacedim>>
+          &jacobian_pushed_forward_2nd_derivatives)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
+          {
+            const unsigned int n_q_points =
+              jacobian_pushed_forward_2nd_derivatives.size();
+
+            double tmp[spacedim][spacedim][spacedim][spacedim];
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<3, dim> *third =
+                  &data.third_derivative(point + data_set, 0);
+
+                DerivativeForm<3, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          for (unsigned int m = 0; m < dim; ++m)
+                            result[fe_to_real[comp_k]][j][l][m] +=
+                              (third[k][j][l][m] * data.local_dof_values[k]);
+                  }
+
+                // push forward the j-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        {
+                          jacobian_pushed_forward_2nd_derivatives
+                            [point][i][j][l][m] =
+                              result[i][0][l][m] * data.covariant[point][j][0];
+                          for (unsigned int jr = 1; jr < dim; ++jr)
+                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
+                                                                   [l][m] +=
+                              result[i][jr][l][m] *
+                              data.covariant[point][j][jr];
+                        }
+
+                // push forward the l-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        {
+                          tmp[i][j][l][m] =
+                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
+                                                                   [0][m] *
+                            data.covariant[point][l][0];
+                          for (unsigned int lr = 1; lr < dim; ++lr)
+                            tmp[i][j][l][m] +=
+                              jacobian_pushed_forward_2nd_derivatives[point][i]
+                                                                     [j][lr]
+                                                                     [m] *
+                              data.covariant[point][l][lr];
+                        }
+
+                // push forward the m-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      for (unsigned int m = 0; m < spacedim; ++m)
+                        {
+                          jacobian_pushed_forward_2nd_derivatives
+                            [point][i][j][l][m] =
+                              tmp[i][j][l][0] * data.covariant[point][m][0];
+                          for (unsigned int mr = 1; mr < dim; ++mr)
+                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
+                                                                   [l][m] +=
+                              tmp[i][j][l][mr] * data.covariant[point][m][mr];
+                        }
+              }
+          }
+      }
+
+      /**
+       * Update the fourth derivative of the transformation from unit to real
+       * cell, the Jacobian hessian gradients.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_3rd_derivatives(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                                &data,
+        const FiniteElement<dim, spacedim>            &fe,
+        const ComponentMask                           &fe_mask,
+        const std::vector<unsigned int>               &fe_to_real,
+        std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_3rd_derivatives)
+          {
+            const unsigned int n_q_points = jacobian_3rd_derivatives.size();
+
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<4, dim> *fourth =
+                  &data.fourth_derivative(point + data_set, 0);
+
+                DerivativeForm<4, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          for (unsigned int m = 0; m < dim; ++m)
+                            for (unsigned int n = 0; n < dim; ++n)
+                              result[fe_to_real[comp_k]][j][l][m][n] +=
+                                (fourth[k][j][l][m][n] *
+                                 data.local_dof_values[k]);
+                  }
+
+                // never touch any data for j,l,m,n=dim in case
+                // dim<spacedim, so it will always be zero as it was
+                // initialized
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < dim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        for (unsigned int n = 0; n < dim; ++n)
+                          jacobian_3rd_derivatives[point][i][j][l][m][n] =
+                            result[i][j][l][m][n];
+              }
+          }
+      }
+
+      /**
+       * Update the fourth derivative of the transformation from unit to real
+       * cell, the Jacobian hessian gradients, pushed forward to the real cell
+       * coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_update_jacobian_pushed_forward_3rd_derivatives(
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real,
+        std::vector<Tensor<5, spacedim>>
+          &jacobian_pushed_forward_3rd_derivatives)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
+          {
+            const unsigned int n_q_points =
+              jacobian_pushed_forward_3rd_derivatives.size();
+
+            double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
+            for (unsigned int point = 0; point < n_q_points; ++point)
+              {
+                const Tensor<4, dim> *fourth =
+                  &data.fourth_derivative(point + data_set, 0);
+
+                DerivativeForm<4, dim, spacedim> result;
+
+                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                  {
+                    const unsigned int comp_k =
+                      fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          for (unsigned int m = 0; m < dim; ++m)
+                            for (unsigned int n = 0; n < dim; ++n)
+                              result[fe_to_real[comp_k]][j][l][m][n] +=
+                                (fourth[k][j][l][m][n] *
+                                 data.local_dof_values[k]);
+                  }
+
+                // push-forward the j-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < dim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        for (unsigned int n = 0; n < dim; ++n)
+                          {
+                            tmp[i][j][l][m][n] = result[i][0][l][m][n] *
+                                                 data.covariant[point][j][0];
+                            for (unsigned int jr = 1; jr < dim; ++jr)
+                              tmp[i][j][l][m][n] +=
+                                result[i][jr][l][m][n] *
+                                data.covariant[point][j][jr];
+                          }
+
+                // push-forward the l-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      for (unsigned int m = 0; m < dim; ++m)
+                        for (unsigned int n = 0; n < dim; ++n)
+                          {
+                            jacobian_pushed_forward_3rd_derivatives
+                              [point][i][j][l][m][n] =
+                                tmp[i][j][0][m][n] *
+                                data.covariant[point][l][0];
+                            for (unsigned int lr = 1; lr < dim; ++lr)
+                              jacobian_pushed_forward_3rd_derivatives[point][i]
+                                                                     [j][l][m]
+                                                                     [n] +=
+                                tmp[i][j][lr][m][n] *
+                                data.covariant[point][l][lr];
+                          }
+
+                // push-forward the m-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      for (unsigned int m = 0; m < spacedim; ++m)
+                        for (unsigned int n = 0; n < dim; ++n)
+                          {
+                            tmp[i][j][l][m][n] =
+                              jacobian_pushed_forward_3rd_derivatives[point][i]
+                                                                     [j][l][0]
+                                                                     [n] *
+                              data.covariant[point][m][0];
+                            for (unsigned int mr = 1; mr < dim; ++mr)
+                              tmp[i][j][l][m][n] +=
+                                jacobian_pushed_forward_3rd_derivatives[point]
+                                                                       [i][j][l]
+                                                                       [mr][n] *
+                                data.covariant[point][m][mr];
+                          }
+
+                // push-forward the n-coordinate
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  for (unsigned int j = 0; j < spacedim; ++j)
+                    for (unsigned int l = 0; l < spacedim; ++l)
+                      for (unsigned int m = 0; m < spacedim; ++m)
+                        for (unsigned int n = 0; n < spacedim; ++n)
+                          {
+                            jacobian_pushed_forward_3rd_derivatives
+                              [point][i][j][l][m][n] =
+                                tmp[i][j][l][m][0] *
+                                data.covariant[point][n][0];
+                            for (unsigned int nr = 1; nr < dim; ++nr)
+                              jacobian_pushed_forward_3rd_derivatives[point][i]
+                                                                     [j][l][m]
+                                                                     [n] +=
+                                tmp[i][j][l][m][nr] *
+                                data.covariant[point][n][nr];
+                          }
+              }
+          }
+      }
+
+
+      /**
+       * Depending on what information is called for in the update flags of the
+       * @p data object, compute the various pieces of information that is
+       * required by the fill_fe_face_values() and fill_fe_subface_values()
+       * functions.  This function simply unifies the work that would be done by
+       * those two functions.
+       *
+       * The resulting data is put into the @p output_data argument.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      maybe_compute_face_data(
+        const dealii::Mapping<dim, spacedim> &mapping,
+        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+                                                         &cell,
+        const unsigned int                                face_no,
+        const unsigned int                                subface_no,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData &data,
+        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+          &output_data)
+      {
+        const UpdateFlags update_flags = data.update_each;
+
+        if (update_flags & update_boundary_forms)
+          {
+            const unsigned int n_q_points = output_data.boundary_forms.size();
+            if (update_flags & update_normal_vectors)
+              AssertDimension(output_data.normal_vectors.size(), n_q_points);
+            if (update_flags & update_JxW_values)
+              AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+            // map the unit tangentials to the real cell. checking for d!=dim-1
+            // eliminates compiler warnings regarding unsigned int expressions <
+            // 0.
+            for (unsigned int d = 0; d != dim - 1; ++d)
+              {
+                Assert(face_no + cell->n_faces() * d <
+                         data.unit_tangentials.size(),
+                       ExcInternalError());
+                Assert(
+                  data.aux[d].size() <=
+                    data.unit_tangentials[face_no + cell->n_faces() * d].size(),
+                  ExcInternalError());
+
+                mapping.transform(
+                  make_array_view(
+                    data.unit_tangentials[face_no + cell->n_faces() * d]),
+                  mapping_contravariant,
+                  data,
+                  make_array_view(data.aux[d]));
+              }
+
+            // if dim==spacedim, we can use the unit tangentials to compute the
+            // boundary form by simply taking the cross product
+            if (dim == spacedim)
+              {
+                for (unsigned int i = 0; i < n_q_points; ++i)
+                  switch (dim)
+                    {
+                      case 1:
+                        // in 1d, we don't have access to any of the data.aux
+                        // fields (because it has only dim-1 components), but we
+                        // can still compute the boundary form by simply looking
+                        // at the number of the face
+                        output_data.boundary_forms[i][0] =
+                          (face_no == 0 ? -1 : +1);
+                        break;
+                      case 2:
+                        output_data.boundary_forms[i] =
+                          cross_product_2d(data.aux[0][i]);
+                        break;
+                      case 3:
+                        output_data.boundary_forms[i] =
+                          cross_product_3d(data.aux[0][i], data.aux[1][i]);
+                        break;
+                      default:
+                        DEAL_II_NOT_IMPLEMENTED();
+                    }
+              }
+            else //(dim < spacedim)
+              {
+                // in the codim-one case, the boundary form results from the
+                // cross product of all the face tangential vectors and the cell
+                // normal vector
+                //
+                // to compute the cell normal, use the same method used in
+                // fill_fe_values for cells above
+                AssertDimension(data.contravariant.size(), n_q_points);
+
+                for (unsigned int point = 0; point < n_q_points; ++point)
+                  {
+                    if (dim == 1)
+                      {
+                        // J is a tangent vector
+                        output_data.boundary_forms[point] =
+                          data.contravariant[point].transpose()[0];
+                        output_data.boundary_forms[point] /=
+                          (face_no == 0 ? -1. : +1.) *
+                          output_data.boundary_forms[point].norm();
+                      }
+
+                    if (dim == 2)
+                      {
+                        const DerivativeForm<1, spacedim, dim> DX_t =
+                          data.contravariant[point].transpose();
+
+                        Tensor<1, spacedim> cell_normal =
+                          cross_product_3d(DX_t[0], DX_t[1]);
+                        cell_normal /= cell_normal.norm();
+
+                        // then compute the face normal from the face tangent
+                        // and the cell normal:
+                        output_data.boundary_forms[point] =
+                          cross_product_3d(data.aux[0][point], cell_normal);
+                      }
+                  }
+              }
+
+            if (update_flags & (update_normal_vectors | update_JxW_values))
+              for (unsigned int i = 0; i < output_data.boundary_forms.size();
+                   ++i)
+                {
+                  if (update_flags & update_JxW_values)
+                    {
+                      output_data.JxW_values[i] =
+                        output_data.boundary_forms[i].norm() *
+                        data.quadrature_weights[i + data_set];
+
+                      if (subface_no != numbers::invalid_unsigned_int)
+                        {
+                          // TODO
+                          const double area_ratio =
+                            GeometryInfo<dim>::subface_ratio(
+                              cell->subface_case(face_no), subface_no);
+                          output_data.JxW_values[i] *= area_ratio;
+                        }
+                    }
+
+                  if (update_flags & update_normal_vectors)
+                    output_data.normal_vectors[i] =
+                      Point<spacedim>(output_data.boundary_forms[i] /
+                                      output_data.boundary_forms[i].norm());
+                }
+          }
+      }
+
+      /**
+       * Do the work of MappingFEField::fill_fe_face_values() and
+       * MappingFEField::fill_fe_subface_values() in a generic way, using the
+       * 'data_set' to differentiate whether we will work on a face (and if so,
+       * which one) or subface.
+       */
+      template <int dim, int spacedim, typename VectorType>
+      void
+      do_fill_fe_face_values(
+        const dealii::Mapping<dim, spacedim> &mapping,
+        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+                                                                 &cell,
+        const unsigned int                                        face_no,
+        const unsigned int                                        subface_no,
+        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData                     &data,
+        const FiniteElement<dim, spacedim> &fe,
+        const ComponentMask                &fe_mask,
+        const std::vector<unsigned int>    &fe_to_real,
+        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+          &output_data)
+      {
+        maybe_compute_q_points<dim, spacedim, VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.quadrature_points);
+
+        maybe_update_Jacobians<dim, spacedim, VectorType>(
+          data_set, data, fe, fe_mask, fe_to_real);
+
+        const UpdateFlags  update_flags = data.update_each;
+        const unsigned int n_q_points   = data.contravariant.size();
+
+        if (update_flags & update_jacobians)
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            output_data.jacobians[point] = data.contravariant[point];
+
+        if (update_flags & update_inverse_jacobians)
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            output_data.inverse_jacobians[point] =
+              data.covariant[point].transpose();
+
+        maybe_update_jacobian_grads<dim, spacedim, VectorType>(
+          data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
+
+        maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_grads);
+
+        maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_2nd_derivatives);
+
+        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
+                                                             spacedim,
+                                                             VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_2nd_derivatives);
+
+        maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_3rd_derivatives);
+
+        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
+                                                             spacedim,
+                                                             VectorType>(
+          data_set,
+          data,
+          fe,
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_3rd_derivatives);
+
+        maybe_compute_face_data<dim, spacedim, VectorType>(
+          mapping, cell, face_no, subface_no, data_set, data, output_data);
+      }
+    } // namespace
+  }   // namespace MappingFEFieldImplementation
+} // namespace internal
+
+
+// Note that the CellSimilarity flag is modifiable, since MappingFEField can
+// need to recalculate data even when cells are similar.
+template <int dim, int spacedim, typename VectorType>
+CellSimilarity::Similarity
+MappingFEField<dim, spacedim, VectorType>::fill_fe_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const CellSimilarity::Similarity,
+  const Quadrature<dim>                                   &quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  const unsigned int n_q_points = quadrature.size();
+
+  update_internal_dofs(cell, data);
+
+  internal::MappingFEFieldImplementation::
+    maybe_compute_q_points<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.quadrature_points);
+
+  internal::MappingFEFieldImplementation::
+    maybe_update_Jacobians<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real);
+
+  const UpdateFlags          update_flags = data.update_each;
+  const std::vector<double> &weights      = quadrature.get_weights();
+
+  // Multiply quadrature weights by absolute value of Jacobian determinants or
+  // the area element g=sqrt(DX^t DX) in case of codim > 0
+
+  if (update_flags & (update_normal_vectors | update_JxW_values))
+    {
+      AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+      Assert(!(update_flags & update_normal_vectors) ||
+               (output_data.normal_vectors.size() == n_q_points),
+             ExcDimensionMismatch(output_data.normal_vectors.size(),
+                                  n_q_points));
+
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          if (dim == spacedim)
+            {
+              const double det = data.volume_elements[point];
+
+              // check for distorted cells.
+
+              // TODO: this allows for anisotropies of up to 1e6 in 3d and
+              // 1e12 in 2d. might want to find a finer
+              // (dimension-independent) criterion
+              Assert(det > 1e-12 * Utilities::fixed_power<dim>(
+                                     cell->diameter() / std::sqrt(double(dim))),
+                     (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
+                       cell->center(), det, point)));
+              output_data.JxW_values[point] = weights[point] * det;
+            }
+          // if dim==spacedim, then there is no cell normal to
+          // compute. since this is for FEValues (and not FEFaceValues),
+          // there are also no face normals to compute
+          else // codim>0 case
+            {
+              Tensor<1, spacedim> DX_t[dim];
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < dim; ++j)
+                  DX_t[j][i] = data.contravariant[point][i][j];
+
+              Tensor<2, dim> G; // First fundamental form
+              for (unsigned int i = 0; i < dim; ++i)
+                for (unsigned int j = 0; j < dim; ++j)
+                  G[i][j] = DX_t[i] * DX_t[j];
+
+              output_data.JxW_values[point] =
+                std::sqrt(determinant(G)) * weights[point];
+
+              if (update_flags & update_normal_vectors)
+                {
+                  Assert(spacedim - dim == 1,
+                         ExcMessage("There is no cell normal in codim 2."));
+
+                  if (dim == 1)
+                    output_data.normal_vectors[point] =
+                      cross_product_2d(-DX_t[0]);
+                  else
+                    {
+                      Assert(dim == 2, ExcInternalError());
+
+                      // dim-1==1 for the second argument, but this
+                      // avoids a compiler warning about array bounds:
+                      output_data.normal_vectors[point] =
+                        cross_product_3d(DX_t[0], DX_t[dim - 1]);
+                    }
+
+                  output_data.normal_vectors[point] /=
+                    output_data.normal_vectors[point].norm();
+
+                  if (cell->direction_flag() == false)
+                    output_data.normal_vectors[point] *= -1.;
+                }
+            } // codim>0 case
+        }
+    }
+
+  // copy values from InternalData to vector given by reference
+  if (update_flags & update_jacobians)
+    {
+      AssertDimension(output_data.jacobians.size(), n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        output_data.jacobians[point] = data.contravariant[point];
+    }
+
+  // copy values from InternalData to vector given by reference
+  if (update_flags & update_inverse_jacobians)
+    {
+      AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        output_data.inverse_jacobians[point] =
+          data.covariant[point].transpose();
+    }
+
+  // calculate derivatives of the Jacobians
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_grads<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_grads);
+
+  // calculate derivatives of the Jacobians pushed forward to real cell
+  // coordinates
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_pushed_forward_grads);
+
+  // calculate hessians of the Jacobians
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_2nd_derivatives);
+
+  // calculate hessians of the Jacobians pushed forward to real cell coordinates
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
+                                                         spacedim,
+                                                         VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_pushed_forward_2nd_derivatives);
+
+  // calculate gradients of the hessians of the Jacobians
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_3rd_derivatives);
+
+  // calculate gradients of the hessians of the Jacobians pushed forward to real
+  // cell coordinates
+  internal::MappingFEFieldImplementation::
+    maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
+                                                         spacedim,
+                                                         VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.jacobian_pushed_forward_3rd_derivatives);
+
+  return CellSimilarity::invalid_next_cell;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::fill_fe_face_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const hp::QCollection<dim - 1>                             &quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  AssertDimension(quadrature.size(), 1);
+
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  update_internal_dofs(cell, data);
+
+  internal::MappingFEFieldImplementation::
+    do_fill_fe_face_values<dim, spacedim, VectorType>(
+      *this,
+      cell,
+      face_no,
+      numbers::invalid_unsigned_int,
+      QProjector<dim>::DataSetDescriptor::face(reference_cell,
+                                               face_no,
+                                               cell->combined_face_orientation(
+                                                 face_no),
+                                               quadrature[0].size()),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data);
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::fill_fe_subface_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const unsigned int                                          subface_no,
+  const Quadrature<dim - 1>                                  &quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  update_internal_dofs(cell, data);
+
+  internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
+                                                                 spacedim,
+                                                                 VectorType>(
+    *this,
+    cell,
+    face_no,
+    numbers::invalid_unsigned_int,
+    QProjector<dim>::DataSetDescriptor::subface(reference_cell,
+                                                face_no,
+                                                subface_no,
+                                                cell->combined_face_orientation(
+                                                  face_no),
+                                                quadrature.size(),
+                                                cell->subface_case(face_no)),
+    data,
+    euler_dof_handler->get_fe(),
+    fe_mask,
+    fe_to_real,
+    output_data);
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::fill_fe_immersed_surface_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const NonMatching::ImmersedSurfaceQuadrature<dim>          &quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  AssertDimension(dim, spacedim);
+  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  const unsigned int n_q_points = quadrature.size();
+
+  update_internal_dofs(cell, data);
+
+  internal::MappingFEFieldImplementation::
+    maybe_compute_q_points<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real,
+      output_data.quadrature_points);
+
+  internal::MappingFEFieldImplementation::
+    maybe_update_Jacobians<dim, spacedim, VectorType>(
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      euler_dof_handler->get_fe(),
+      fe_mask,
+      fe_to_real);
+
+  const UpdateFlags          update_flags = data.update_each;
+  const std::vector<double> &weights      = quadrature.get_weights();
+
+  if (update_flags & (update_normal_vectors | update_JxW_values))
+    {
+      AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+      Assert(!(update_flags & update_normal_vectors) ||
+               (output_data.normal_vectors.size() == n_q_points),
+             ExcDimensionMismatch(output_data.normal_vectors.size(),
+                                  n_q_points));
+
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          const double det = data.volume_elements[point];
+
+          // check for distorted cells.
+
+          // TODO: this allows for anisotropies of up to 1e6 in 3d and
+          // 1e12 in 2d. might want to find a finer
+          // (dimension-independent) criterion
+          Assert(det > 1e-12 * Utilities::fixed_power<dim>(
+                                 cell->diameter() / std::sqrt(double(dim))),
+                 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
+                   cell->center(), det, point)));
+
+          // The normals are n = J^{-T} * \hat{n} before normalizing.
+          Tensor<1, spacedim> normal;
+          for (unsigned int d = 0; d < spacedim; d++)
+            normal[d] =
+              data.covariant[point][d] * quadrature.normal_vector(point);
+
+          output_data.JxW_values[point] = weights[point] * det * normal.norm();
+
+          if ((update_flags & update_normal_vectors) != 0u)
+            {
+              normal /= normal.norm();
+              output_data.normal_vectors[point] = normal;
+            }
+        }
+
+      // copy values from InternalData to vector given by reference
+      if (update_flags & update_jacobians)
+        {
+          AssertDimension(output_data.jacobians.size(), n_q_points);
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            output_data.jacobians[point] = data.contravariant[point];
+        }
+
+      // copy values from InternalData to vector given by reference
+      if (update_flags & update_inverse_jacobians)
+        {
+          AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            output_data.inverse_jacobians[point] =
+              data.covariant[point].transpose();
+        }
+
+      // calculate derivatives of the Jacobians
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_grads<dim, spacedim, VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_grads);
+
+      // calculate derivatives of the Jacobians pushed forward to real cell
+      // coordinates
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_grads);
+
+      // calculate hessians of the Jacobians
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_2nd_derivatives);
+
+      // calculate hessians of the Jacobians pushed forward to real cell
+      // coordinates
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
+                                                             spacedim,
+                                                             VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_2nd_derivatives);
+
+      // calculate gradients of the hessians of the Jacobians
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_3rd_derivatives);
+
+      // calculate gradients of the hessians of the Jacobians pushed forward to
+      // real cell coordinates
+      internal::MappingFEFieldImplementation::
+        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
+                                                             spacedim,
+                                                             VectorType>(
+          QProjector<dim>::DataSetDescriptor::cell(),
+          data,
+          euler_dof_handler->get_fe(),
+          fe_mask,
+          fe_to_real,
+          output_data.jacobian_pushed_forward_3rd_derivatives);
+    }
+}
+
+namespace internal
+{
+  namespace MappingFEFieldImplementation
+  {
+    namespace
+    {
+      template <int dim, int spacedim, int rank, typename VectorType>
+      void
+      transform_fields(
+        const ArrayView<const Tensor<rank, dim>>                &input,
+        const MappingKind                                        mapping_kind,
+        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+        const ArrayView<Tensor<rank, spacedim>>                 &output)
+      {
+        AssertDimension(input.size(), output.size());
+        Assert((dynamic_cast<
+                  const typename dealii::
+                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
+                  &mapping_data) != nullptr),
+               ExcInternalError());
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData &data = static_cast<
+            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+              InternalData &>(mapping_data);
+
+        switch (mapping_kind)
+          {
+            case mapping_contravariant:
+              {
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  output[i] =
+                    apply_transformation(data.contravariant[i], input[i]);
+
+                return;
+              }
+
+            case mapping_piola:
+              {
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+                Assert(
+                  data.update_each & update_volume_elements,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_volume_elements"));
+                Assert(rank == 1, ExcMessage("Only for rank 1"));
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  {
+                    output[i] =
+                      apply_transformation(data.contravariant[i], input[i]);
+                    output[i] /= data.volume_elements[i];
+                  }
+                return;
+              }
+
+
+            // We still allow this operation as in the
+            // reference cell Derivatives are Tensor
+            // rather than DerivativeForm
+            case mapping_covariant:
+              {
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  output[i] = apply_transformation(data.covariant[i], input[i]);
+
+                return;
+              }
+
+            default:
+              DEAL_II_NOT_IMPLEMENTED();
+          }
+      }
+
+
+      template <int dim, int spacedim, int rank, typename VectorType>
+      void
+      transform_differential_forms(
+        const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
+        const MappingKind                                        mapping_kind,
+        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+        const ArrayView<Tensor<rank + 1, spacedim>>             &output)
+      {
+        AssertDimension(input.size(), output.size());
+        Assert((dynamic_cast<
+                  const typename dealii::
+                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
+                  &mapping_data) != nullptr),
+               ExcInternalError());
+        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+          InternalData &data = static_cast<
+            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+              InternalData &>(mapping_data);
+
+        switch (mapping_kind)
+          {
+            case mapping_covariant:
+              {
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  output[i] = apply_transformation(data.covariant[i], input[i]);
+
+                return;
+              }
+            default:
+              DEAL_II_NOT_IMPLEMENTED();
+          }
+      }
+    } // namespace
+  }   // namespace MappingFEFieldImplementation
+} // namespace internal
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::transform(
+  const ArrayView<const Tensor<1, dim>>                   &input,
+  const MappingKind                                        mapping_kind,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<1, spacedim>>                    &output) const
+{
+  AssertDimension(input.size(), output.size());
+
+  internal::MappingFEFieldImplementation::
+    transform_fields<dim, spacedim, 1, VectorType>(input,
+                                                   mapping_kind,
+                                                   mapping_data,
+                                                   output);
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::transform(
+  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
+  const MappingKind                                        mapping_kind,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<2, spacedim>>                    &output) const
+{
+  AssertDimension(input.size(), output.size());
+
+  internal::MappingFEFieldImplementation::
+    transform_differential_forms<dim, spacedim, 1, VectorType>(input,
+                                                               mapping_kind,
+                                                               mapping_data,
+                                                               output);
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::transform(
+  const ArrayView<const Tensor<2, dim>> &input,
+  const MappingKind,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<2, spacedim>>                    &output) const
+{
+  (void)input;
+  (void)output;
+  (void)mapping_data;
+  AssertDimension(input.size(), output.size());
+
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::transform(
+  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
+  const MappingKind                                        mapping_kind,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<3, spacedim>>                    &output) const
+{
+  AssertDimension(input.size(), output.size());
+  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_kind)
+    {
+      case mapping_covariant_gradient:
+        {
+          Assert(data.update_each & update_contravariant_transformation,
+                 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                   "update_covariant_transformation"));
+
+          for (unsigned int q = 0; q < output.size(); ++q)
+            for (unsigned int i = 0; i < spacedim; ++i)
+              for (unsigned int j = 0; j < spacedim; ++j)
+                for (unsigned int k = 0; k < spacedim; ++k)
+                  {
+                    output[q][i][j][k] = data.covariant[q][j][0] *
+                                         data.covariant[q][k][0] *
+                                         input[q][i][0][0];
+                    for (unsigned int J = 0; J < dim; ++J)
+                      {
+                        const unsigned int K0 = (0 == J) ? 1 : 0;
+                        for (unsigned int K = K0; K < dim; ++K)
+                          output[q][i][j][k] += data.covariant[q][j][J] *
+                                                data.covariant[q][k][K] *
+                                                input[q][i][J][K];
+                      }
+                  }
+          return;
+        }
+
+      default:
+        DEAL_II_NOT_IMPLEMENTED();
+    }
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::transform(
+  const ArrayView<const Tensor<3, dim>> &input,
+  const MappingKind /*mapping_kind*/,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<3, spacedim>>                    &output) const
+{
+  (void)input;
+  (void)output;
+  (void)mapping_data;
+  AssertDimension(input.size(), output.size());
+
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Point<spacedim>
+MappingFEField<dim, spacedim, VectorType>::transform_unit_to_real_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const Point<dim>                                           &p) const
+{
+  //  Use the get_data function to create an InternalData with data vectors of
+  //  the right size and transformation shape values already computed at point
+  //  p.
+  const Quadrature<dim> point_quadrature(p);
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
+    get_data(update_quadrature_points | update_jacobians, point_quadrature));
+  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
+         ExcInternalError());
+
+  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
+
+  return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+Point<spacedim>
+MappingFEField<dim, spacedim, VectorType>::do_transform_unit_to_real_cell(
+  const InternalData &data) const
+{
+  Point<spacedim> p_real;
+
+  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
+    {
+      unsigned int comp_i =
+        euler_dof_handler->get_fe().system_to_component_index(i).first;
+      if (fe_mask[comp_i])
+        p_real[fe_to_real[comp_i]] +=
+          data.local_dof_values[i] * data.shape(0, i);
+    }
+
+  return p_real;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+Point<dim>
+MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const Point<spacedim>                                      &p) const
+{
+  // first a Newton iteration based on the real mapping. It uses the center
+  // point of the cell as a starting point
+  Point<dim> initial_p_unit;
+  try
+    {
+      initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
+                         .transform_real_to_unit_cell(cell, p);
+    }
+  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
+    {
+      // mirror the conditions of the code below to determine if we need to
+      // use an arbitrary starting point or if we just need to rethrow the
+      // exception
+      for (unsigned int d = 0; d < dim; ++d)
+        initial_p_unit[d] = 0.5;
+    }
+
+  initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
+
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim > dim)
+    update_flags |= update_jacobian_grads;
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
+    get_data(update_flags, Quadrature<dim>(initial_p_unit)));
+  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
+         ExcInternalError());
+
+  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
+
+  return do_transform_real_to_unit_cell(cell,
+                                        p,
+                                        initial_p_unit,
+                                        static_cast<InternalData &>(*mdata));
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+Point<dim>
+MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const Point<spacedim>                                      &p,
+  const Point<dim>                                           &starting_guess,
+  InternalData                                               &mdata) const
+{
+  const unsigned int n_shapes = mdata.shape_values.size();
+  (void)n_shapes;
+  Assert(n_shapes != 0, ExcInternalError());
+  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
+
+
+  // Newton iteration to solve
+  // f(x)=p(x)-p=0
+  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+  // The start value was set to be the
+  // linear approximation to the cell
+  // The shape values and derivatives
+  // of the mapping at this point are
+  // previously computed.
+
+  Point<dim> p_unit = starting_guess;
+  Point<dim> f;
+  mdata.reinit(mdata.update_each, Quadrature<dim>(starting_guess));
+
+  Point<spacedim>     p_real(do_transform_unit_to_real_cell(mdata));
+  Tensor<1, spacedim> p_minus_F              = p - p_real;
+  const double        eps                    = 1.e-12 * cell->diameter();
+  const unsigned int  newton_iteration_limit = 20;
+  unsigned int        newton_iteration       = 0;
+  while (p_minus_F.norm_square() > eps * eps)
+    {
+      // f'(x)
+      Point<spacedim> DF[dim];
+      Tensor<2, dim>  df;
+      for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
+        {
+          const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
+          const unsigned int    comp_k =
+            euler_dof_handler->get_fe().system_to_component_index(k).first;
+          if (fe_mask[comp_k])
+            for (unsigned int j = 0; j < dim; ++j)
+              DF[j][fe_to_real[comp_k]] +=
+                mdata.local_dof_values[k] * grad_k[j];
+        }
+      for (unsigned int j = 0; j < dim; ++j)
+        {
+          f[j] = DF[j] * p_minus_F;
+          for (unsigned int l = 0; l < dim; ++l)
+            df[j][l] = -DF[j] * DF[l];
+        }
+      // Solve  [f'(x)]d=f(x)
+      const Tensor<1, dim> delta =
+        invert(df) * static_cast<const Tensor<1, dim> &>(f);
+      // do a line search
+      double step_length = 1;
+      do
+        {
+          // update of p_unit. The
+          // spacedimth component of
+          // transformed point is simply
+          // ignored in codimension one
+          // case. When this component is
+          // not zero, then we are
+          // projecting the point to the
+          // surface or curve identified
+          // by the cell.
+          Point<dim> p_unit_trial = p_unit;
+          for (unsigned int i = 0; i < dim; ++i)
+            p_unit_trial[i] -= step_length * delta[i];
+          // shape values and derivatives
+          // at new p_unit point
+          mdata.reinit(mdata.update_each, Quadrature<dim>(p_unit_trial));
+          // f(x)
+          const Point<spacedim> p_real_trial =
+            do_transform_unit_to_real_cell(mdata);
+          const Tensor<1, spacedim> f_trial = p - p_real_trial;
+          // see if we are making progress with the current step length
+          // and if not, reduce it by a factor of two and try again
+          if (f_trial.norm() < p_minus_F.norm())
+            {
+              p_real    = p_real_trial;
+              p_unit    = p_unit_trial;
+              p_minus_F = f_trial;
+              break;
+            }
+          else if (step_length > 0.05)
+            step_length /= 2;
+          else
+            goto failure;
+        }
+      while (true);
+      ++newton_iteration;
+      if (newton_iteration > newton_iteration_limit)
+        goto failure;
+    }
+  return p_unit;
+  // if we get to the following label, then we have either run out
+  // of Newton iterations, or the line search has not converged.
+  // in either case, we need to give up, so throw an exception that
+  // can then be caught
+failure:
+  AssertThrow(false,
+              (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
+  // ...the compiler wants us to return something, though we can
+  // of course never get here...
+  return {};
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+unsigned int
+MappingFEField<dim, spacedim, VectorType>::get_degree() const
+{
+  return euler_dof_handler->get_fe().degree;
+}
+
+
+
+template <int dim, int spacedim, typename VectorType>
+ComponentMask
+MappingFEField<dim, spacedim, VectorType>::get_component_mask() const
+{
+  return this->fe_mask;
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+std::unique_ptr<Mapping<dim, spacedim>>
+MappingFEField<dim, spacedim, VectorType>::clone() const
+{
+  return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
+}
+
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::update_internal_dofs(
+  const typename Triangulation<dim, spacedim>::cell_iterator             &cell,
+  const typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
+  const
+{
+  Assert(euler_dof_handler != nullptr,
+         ExcMessage("euler_dof_handler is empty"));
+
+  typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
+                                                             euler_dof_handler);
+  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
+  if (uses_level_dofs)
+    {
+      AssertIndexRange(cell->level(), euler_vector.size());
+      AssertDimension(euler_vector[cell->level()]->size(),
+                      euler_dof_handler->n_dofs(cell->level()));
+    }
+  else
+    AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
+
+  if (uses_level_dofs)
+    dof_cell->get_mg_dof_indices(data.local_dof_indices);
+  else
+    dof_cell->get_dof_indices(data.local_dof_indices);
+
+  const VectorType &vector =
+    uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
+
+  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
+    data.local_dof_values[i] =
+      internal::ElementAccess<VectorType>::get(vector,
+                                               data.local_dof_indices[i]);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 8028a37918d8c8864eea496f3752c1a95c9c5fb4..3457b5deecc0fe896d002ff5a56e199afdf805e1 100644 (file)
 //
 // ------------------------------------------------------------------------
 
-#include <deal.II/base/array_view.h>
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/polynomial.h>
-#include <deal.II/base/qprojector.h>
-#include <deal.II/base/quadrature.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/tensor_product_polynomials.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/dofs/dof_accessor.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/fe/mapping_fe_field.h>
-
-#include <deal.II/grid/tria_iterator.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/petsc_block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/lac/trilinos_epetra_vector.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_tpetra_block_vector.h>
-#include <deal.II/lac/trilinos_tpetra_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/numerics/vector_tools.h>
-
-#include <fstream>
-#include <memory>
-#include <numeric>
-
 
+#include <deal.II/fe/mapping_fe_field.templates.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-
-template <int dim, int spacedim, typename VectorType>
-MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
-  const FiniteElement<dim, spacedim> &fe,
-  const ComponentMask                &mask)
-  : fe(&fe)
-  , unit_tangentials()
-  , n_shape_functions(fe.n_dofs_per_cell())
-  , mask(mask)
-  , local_dof_indices(fe.n_dofs_per_cell())
-  , local_dof_values(fe.n_dofs_per_cell())
-{}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::InternalData::reinit(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &quadrature)
-{
-  // store the flags in the internal data object so we can access them
-  // in fill_fe_*_values(). use the transitive hull of the required
-  // flags
-  this->update_each = update_flags;
-
-  const unsigned int             n_q_points = quadrature.size();
-  const std::vector<Point<dim>> &points     = quadrature.get_points();
-
-  // see if we need the (transformation) shape function values
-  // and/or gradients and resize the necessary arrays
-  if (update_flags & update_quadrature_points)
-    {
-      shape_values.resize(n_shape_functions * n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        for (unsigned int i = 0; i < n_shape_functions; ++i)
-          shape(point, i) = fe->shape_value(i, points[point]);
-    }
-
-  if (update_flags &
-      (update_covariant_transformation | update_contravariant_transformation |
-       update_JxW_values | update_boundary_forms | update_normal_vectors |
-       update_jacobians | update_jacobian_grads | update_inverse_jacobians))
-    {
-      shape_derivatives.resize(n_shape_functions * n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        for (unsigned int i = 0; i < n_shape_functions; ++i)
-          derivative(point, i) = fe->shape_grad(i, points[point]);
-    }
-
-  if (update_flags & update_covariant_transformation)
-    covariant.resize(n_q_points);
-
-  if (update_flags & update_contravariant_transformation)
-    contravariant.resize(n_q_points);
-
-  if (update_flags & update_volume_elements)
-    volume_elements.resize(n_q_points);
-
-  if (update_flags &
-      (update_jacobian_grads | update_jacobian_pushed_forward_grads))
-    {
-      shape_second_derivatives.resize(n_shape_functions * n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        for (unsigned int i = 0; i < n_shape_functions; ++i)
-          second_derivative(point, i) = fe->shape_grad_grad(i, points[point]);
-    }
-
-  if (update_flags & (update_jacobian_2nd_derivatives |
-                      update_jacobian_pushed_forward_2nd_derivatives))
-    {
-      shape_third_derivatives.resize(n_shape_functions * n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        for (unsigned int i = 0; i < n_shape_functions; ++i)
-          third_derivative(point, i) =
-            fe->shape_3rd_derivative(i, points[point]);
-    }
-
-  if (update_flags & (update_jacobian_3rd_derivatives |
-                      update_jacobian_pushed_forward_3rd_derivatives))
-    {
-      shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        for (unsigned int i = 0; i < n_shape_functions; ++i)
-          fourth_derivative(point, i) =
-            fe->shape_4th_derivative(i, points[point]);
-    }
-
-  // This (for face values and simplices) can be different for different
-  // calls, so always copy
-  quadrature_weights = quadrature.get_weights();
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-std::size_t
-MappingFEField<dim, spacedim, VectorType>::InternalData::memory_consumption()
-  const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-  return 0;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-double &
-MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
-  return shape_values[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-const Tensor<1, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_derivatives.size());
-  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Tensor<1, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_derivatives.size());
-  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-const Tensor<2, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_second_derivatives.size());
-  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Tensor<2, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_second_derivatives.size());
-  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-const Tensor<3, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_third_derivatives.size());
-  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Tensor<3, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_third_derivatives.size());
-  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-const Tensor<4, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_fourth_derivatives.size());
-  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Tensor<4, dim> &
-MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_fourth_derivatives.size());
-  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-MappingFEField<dim, spacedim, VectorType>::MappingFEField(
-  const DoFHandler<dim, spacedim> &euler_dof_handler,
-  const VectorType                &euler_vector,
-  const ComponentMask             &mask)
-  : reference_cell(euler_dof_handler.get_fe().reference_cell())
-  , uses_level_dofs(false)
-  , euler_vector({&euler_vector})
-  , euler_dof_handler(&euler_dof_handler)
-  , fe_mask(mask.size() != 0u ?
-              mask :
-              ComponentMask(
-                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
-                true))
-  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
-  , fe_values(this->euler_dof_handler->get_fe(),
-              reference_cell.template get_nodal_type_quadrature<dim>(),
-              update_values)
-{
-  AssertDimension(euler_dof_handler.n_dofs(), euler_vector.size());
-  unsigned int size = 0;
-  for (unsigned int i = 0; i < fe_mask.size(); ++i)
-    {
-      if (fe_mask[i])
-        fe_to_real[i] = size++;
-    }
-  AssertDimension(size, spacedim);
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-MappingFEField<dim, spacedim, VectorType>::MappingFEField(
-  const DoFHandler<dim, spacedim> &euler_dof_handler,
-  const std::vector<VectorType>   &euler_vector,
-  const ComponentMask             &mask)
-  : reference_cell(euler_dof_handler.get_fe().reference_cell())
-  , uses_level_dofs(true)
-  , euler_dof_handler(&euler_dof_handler)
-  , fe_mask(mask.size() != 0u ?
-              mask :
-              ComponentMask(
-                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
-                true))
-  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
-  , fe_values(this->euler_dof_handler->get_fe(),
-              reference_cell.template get_nodal_type_quadrature<dim>(),
-              update_values)
-{
-  unsigned int size = 0;
-  for (unsigned int i = 0; i < fe_mask.size(); ++i)
-    {
-      if (fe_mask[i])
-        fe_to_real[i] = size++;
-    }
-  AssertDimension(size, spacedim);
-
-  Assert(euler_dof_handler.has_level_dofs(),
-         ExcMessage("The underlying DoFHandler object did not call "
-                    "distribute_mg_dofs(). In this case, the construction via "
-                    "level vectors does not make sense."));
-  AssertDimension(euler_vector.size(),
-                  euler_dof_handler.get_triangulation().n_global_levels());
-  this->euler_vector.clear();
-  this->euler_vector.resize(euler_vector.size());
-  for (unsigned int i = 0; i < euler_vector.size(); ++i)
-    {
-      AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
-      this->euler_vector[i] = &euler_vector[i];
-    }
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-MappingFEField<dim, spacedim, VectorType>::MappingFEField(
-  const DoFHandler<dim, spacedim> &euler_dof_handler,
-  const MGLevelObject<VectorType> &euler_vector,
-  const ComponentMask             &mask)
-  : reference_cell(euler_dof_handler.get_fe().reference_cell())
-  , uses_level_dofs(true)
-  , euler_dof_handler(&euler_dof_handler)
-  , fe_mask(mask.size() != 0u ?
-              mask :
-              ComponentMask(
-                euler_dof_handler.get_fe().get_nonzero_components(0).size(),
-                true))
-  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
-  , fe_values(this->euler_dof_handler->get_fe(),
-              reference_cell.template get_nodal_type_quadrature<dim>(),
-              update_values)
-{
-  unsigned int size = 0;
-  for (unsigned int i = 0; i < fe_mask.size(); ++i)
-    {
-      if (fe_mask[i])
-        fe_to_real[i] = size++;
-    }
-  AssertDimension(size, spacedim);
-
-  Assert(euler_dof_handler.has_level_dofs(),
-         ExcMessage("The underlying DoFHandler object did not call "
-                    "distribute_mg_dofs(). In this case, the construction via "
-                    "level vectors does not make sense."));
-  AssertDimension(euler_vector.max_level() + 1,
-                  euler_dof_handler.get_triangulation().n_global_levels());
-  this->euler_vector.clear();
-  this->euler_vector.resize(
-    euler_dof_handler.get_triangulation().n_global_levels());
-  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
-       ++i)
-    {
-      AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
-      this->euler_vector[i] = &euler_vector[i];
-    }
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-MappingFEField<dim, spacedim, VectorType>::MappingFEField(
-  const MappingFEField<dim, spacedim, VectorType> &mapping)
-  : reference_cell(mapping.reference_cell)
-  , uses_level_dofs(mapping.uses_level_dofs)
-  , euler_vector(mapping.euler_vector)
-  , euler_dof_handler(mapping.euler_dof_handler)
-  , fe_mask(mapping.fe_mask)
-  , fe_to_real(mapping.fe_to_real)
-  , fe_values(mapping.euler_dof_handler->get_fe(),
-              reference_cell.template get_nodal_type_quadrature<dim>(),
-              update_values)
-{}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-inline const double &
-MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
-  return shape_values[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-bool
-MappingFEField<dim, spacedim, VectorType>::preserves_vertex_locations() const
-{
-  return false;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-bool
-MappingFEField<dim, spacedim, VectorType>::is_compatible_with(
-  const ReferenceCell &reference_cell) const
-{
-  Assert(dim == reference_cell.get_dimension(),
-         ExcMessage("The dimension of your mapping (" +
-                    Utilities::to_string(dim) +
-                    ") and the reference cell cell_type (" +
-                    Utilities::to_string(reference_cell.get_dimension()) +
-                    " ) do not agree."));
-
-  return this->reference_cell == reference_cell;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-boost::container::small_vector<Point<spacedim>,
-#ifndef _MSC_VER
-                               ReferenceCells::max_n_vertices<dim>()
-#else
-                               GeometryInfo<dim>::vertices_per_cell
-#endif
-                               >
-MappingFEField<dim, spacedim, VectorType>::get_vertices(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
-{
-  // we transform our tria iterator into a dof iterator so we can access
-  // data not associated with triangulations
-  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
-    *cell, euler_dof_handler);
-
-  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
-  AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
-  AssertDimension(fe_to_real.size(),
-                  euler_dof_handler->get_fe().n_components());
-  if (uses_level_dofs)
-    {
-      AssertIndexRange(cell->level(), euler_vector.size());
-      AssertDimension(euler_vector[cell->level()]->size(),
-                      euler_dof_handler->n_dofs(cell->level()));
-    }
-  else
-    AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
-
-  {
-    std::lock_guard<std::mutex> lock(fe_values_mutex);
-    fe_values.reinit(dof_cell);
-  }
-  const unsigned int dofs_per_cell =
-    euler_dof_handler->get_fe().n_dofs_per_cell();
-  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
-  if (uses_level_dofs)
-    dof_cell->get_mg_dof_indices(dof_indices);
-  else
-    dof_cell->get_dof_indices(dof_indices);
-
-  const VectorType &vector =
-    uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
-
-  boost::container::small_vector<Point<spacedim>,
-#ifndef _MSC_VER
-                                 ReferenceCells::max_n_vertices<dim>()
-#else
-                                 GeometryInfo<dim>::vertices_per_cell
-#endif
-                                 >
-    vertices(cell->n_vertices());
-  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-    {
-      const unsigned int comp = fe_to_real
-        [euler_dof_handler->get_fe().system_to_component_index(i).first];
-      if (comp != numbers::invalid_unsigned_int)
-        {
-          typename VectorType::value_type value =
-            internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
-          if (euler_dof_handler->get_fe().is_primitive(i))
-            for (const unsigned int v : cell->vertex_indices())
-              vertices[v][comp] += fe_values.shape_value(i, v) * value;
-          else
-            DEAL_II_NOT_IMPLEMENTED();
-        }
-    }
-
-  return vertices;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-UpdateFlags
-MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
-  const UpdateFlags in) const
-{
-  // add flags if the respective quantities are necessary to compute
-  // what we need. note that some flags appear in both conditions and
-  // in subsequent set operations. this leads to some circular
-  // logic. the only way to treat this is to iterate. since there are
-  // 5 if-clauses in the loop, it will take at most 4 iterations to
-  // converge. do them:
-  UpdateFlags out = in;
-  for (unsigned int i = 0; i < 5; ++i)
-    {
-      // The following is a little incorrect:
-      // If not applied on a face,
-      // update_boundary_forms does not
-      // make sense. On the other hand,
-      // it is necessary on a
-      // face. Currently,
-      // update_boundary_forms is simply
-      // ignored for the interior of a
-      // cell.
-      if (out & (update_JxW_values | update_normal_vectors))
-        out |= update_boundary_forms;
-
-      if (out &
-          (update_covariant_transformation | update_jacobian_grads |
-           update_jacobians | update_boundary_forms | update_normal_vectors))
-        out |= update_contravariant_transformation;
-
-      if (out &
-          (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
-           update_jacobian_pushed_forward_2nd_derivatives |
-           update_jacobian_pushed_forward_3rd_derivatives))
-        out |= update_covariant_transformation;
-
-      // The contravariant transformation is used in the Piola
-      // transformation, which requires the determinant of the Jacobi
-      // matrix of the transformation.  Because we have no way of
-      // knowing here whether the finite element wants to use the
-      // contravariant or the Piola transforms, we add the volume elements
-      // to the list of flags to be updated for each cell.
-      if (out & update_contravariant_transformation)
-        out |= update_volume_elements;
-
-      if (out & update_normal_vectors)
-        out |= update_volume_elements;
-    }
-
-  return out;
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::compute_face_data(
-  const unsigned int n_original_q_points,
-  InternalData      &data) const
-{
-  // Set to the size of a single quadrature object for faces, as the size set
-  // in in reinit() is for all points
-  if (data.update_each & update_covariant_transformation)
-    data.covariant.resize(n_original_q_points);
-
-  if (data.update_each & update_contravariant_transformation)
-    data.contravariant.resize(n_original_q_points);
-
-  if (data.update_each & update_volume_elements)
-    data.volume_elements.resize(n_original_q_points);
-
-  if (dim > 1)
-    {
-      if (data.update_each & update_boundary_forms)
-        {
-          data.aux.resize(
-            dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
-
-
-          // TODO: only a single reference cell type possible...
-          const auto n_faces = reference_cell.n_faces();
-
-          // Compute tangentials to the unit cell.
-          for (unsigned int i = 0; i < n_faces; ++i)
-            {
-              data.unit_tangentials[i].resize(n_original_q_points);
-              std::fill(data.unit_tangentials[i].begin(),
-                        data.unit_tangentials[i].end(),
-                        reference_cell.template face_tangent_vector<dim>(i, 0));
-              if (dim > 2)
-                {
-                  data.unit_tangentials[n_faces + i].resize(
-                    n_original_q_points);
-                  std::fill(
-                    data.unit_tangentials[n_faces + i].begin(),
-                    data.unit_tangentials[n_faces + i].end(),
-                    reference_cell.template face_tangent_vector<dim>(i, 1));
-                }
-            }
-        }
-    }
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
-MappingFEField<dim, spacedim, VectorType>::get_data(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &quadrature) const
-{
-  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
-    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  data_ptr->reinit(requires_update_flags(update_flags), quadrature);
-
-  return data_ptr;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
-MappingFEField<dim, spacedim, VectorType>::get_face_data(
-  const UpdateFlags               update_flags,
-  const hp::QCollection<dim - 1> &quadrature) const
-{
-  AssertDimension(quadrature.size(), 1);
-
-  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
-    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-
-  const Quadrature<dim> q(
-    QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
-  data.reinit(requires_update_flags(update_flags), q);
-  this->compute_face_data(quadrature[0].size(), data);
-
-  return data_ptr;
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
-MappingFEField<dim, spacedim, VectorType>::get_subface_data(
-  const UpdateFlags          update_flags,
-  const Quadrature<dim - 1> &quadrature) const
-{
-  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
-    std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-
-  const Quadrature<dim> q(
-    QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
-  data.reinit(requires_update_flags(update_flags), q);
-  this->compute_face_data(quadrature.size(), data);
-
-  return data_ptr;
-}
-
-
-
-namespace internal
-{
-  namespace MappingFEFieldImplementation
-  {
-    namespace
-    {
-      /**
-       * Compute the locations of quadrature points on the object described by
-       * the first argument (and the cell for which the mapping support points
-       * have already been set), but only if the update_flags of the @p data
-       * argument indicate so.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_compute_q_points(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Point<spacedim>>       &quadrature_points)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        if (update_flags & update_quadrature_points)
-          {
-            for (unsigned int point = 0; point < quadrature_points.size();
-                 ++point)
-              {
-                Point<spacedim> result;
-                const double   *shape = &data.shape(point + data_set, 0);
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      result[fe_to_real[comp_k]] +=
-                        data.local_dof_values[k] * shape[k];
-                  }
-
-                quadrature_points[point] = result;
-              }
-          }
-      }
-
-      /**
-       * Update the co- and contravariant matrices as well as their determinant,
-       * for the cell described stored in the data object, but only if the
-       * update_flags of the @p data argument indicate so.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_Jacobians(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        // then Jacobians
-        if (update_flags & update_contravariant_transformation)
-          {
-            const unsigned int n_q_points = data.contravariant.size();
-
-            Assert(data.n_shape_functions > 0, ExcInternalError());
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<1, dim> *data_derv =
-                  &data.derivative(point + data_set, 0);
-
-                Tensor<1, dim> result[spacedim];
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      result[fe_to_real[comp_k]] +=
-                        data.local_dof_values[k] * data_derv[k];
-                  }
-
-                // write result into contravariant data
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  {
-                    data.contravariant[point][i] = result[i];
-                  }
-              }
-          }
-
-        if (update_flags & update_covariant_transformation)
-          {
-            AssertDimension(data.covariant.size(), data.contravariant.size());
-            for (unsigned int point = 0; point < data.contravariant.size();
-                 ++point)
-              data.covariant[point] =
-                (data.contravariant[point]).covariant_form();
-          }
-
-        if (update_flags & update_volume_elements)
-          {
-            AssertDimension(data.contravariant.size(),
-                            data.volume_elements.size());
-            for (unsigned int point = 0; point < data.contravariant.size();
-                 ++point)
-              data.volume_elements[point] =
-                data.contravariant[point].determinant();
-          }
-      }
-
-      /**
-       * Update the Hessian of the transformation from unit to real cell, the
-       * Jacobian gradients.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_grads(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_grads)
-          {
-            const unsigned int n_q_points = jacobian_grads.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<2, dim> *second =
-                  &data.second_derivative(point + data_set, 0);
-
-                DerivativeForm<2, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          result[fe_to_real[comp_k]][j][l] +=
-                            (second[k][j][l] * data.local_dof_values[k]);
-                  }
-
-                // never touch any data for j=dim in case dim<spacedim, so
-                // it will always be zero as it was initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      jacobian_grads[point][i][j][l] = result[i][j][l];
-              }
-          }
-      }
-
-      /**
-       * Update the Hessian of the transformation from unit to real cell, the
-       * Jacobian gradients, pushed forward to the real cell coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_grads(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<3, spacedim>>   &jacobian_pushed_forward_grads)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_grads)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_grads.size();
-
-            double tmp[spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<2, dim> *second =
-                  &data.second_derivative(point + data_set, 0);
-
-                DerivativeForm<2, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          result[fe_to_real[comp_k]][j][l] +=
-                            (second[k][j][l] * data.local_dof_values[k]);
-                  }
-
-                // first push forward the j-components
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      {
-                        tmp[i][j][l] =
-                          result[i][0][l] * data.covariant[point][j][0];
-                        for (unsigned int jr = 1; jr < dim; ++jr)
-                          {
-                            tmp[i][j][l] +=
-                              result[i][jr][l] * data.covariant[point][j][jr];
-                          }
-                      }
-
-                // now, pushing forward the l-components
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      {
-                        jacobian_pushed_forward_grads[point][i][j][l] =
-                          tmp[i][j][0] * data.covariant[point][l][0];
-                        for (unsigned int lr = 1; lr < dim; ++lr)
-                          {
-                            jacobian_pushed_forward_grads[point][i][j][l] +=
-                              tmp[i][j][lr] * data.covariant[point][l][lr];
-                          }
-                      }
-              }
-          }
-      }
-
-      /**
-       * Update the third derivative of the transformation from unit to real
-       * cell, the Jacobian hessians.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_2nd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_2nd_derivatives)
-          {
-            const unsigned int n_q_points = jacobian_2nd_derivatives.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<3, dim> *third =
-                  &data.third_derivative(point + data_set, 0);
-
-                DerivativeForm<3, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            result[fe_to_real[comp_k]][j][l][m] +=
-                              (third[k][j][l][m] * data.local_dof_values[k]);
-                  }
-
-                // never touch any data for j=dim in case dim<spacedim, so
-                // it will always be zero as it was initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        jacobian_2nd_derivatives[point][i][j][l][m] =
-                          result[i][j][l][m];
-              }
-          }
-      }
-
-      /**
-       * Update the third derivative of the transformation from unit to real
-       * cell, the Jacobian hessians, pushed forward to the real cell
-       * coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_2nd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<4, spacedim>>
-          &jacobian_pushed_forward_2nd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_2nd_derivatives.size();
-
-            double tmp[spacedim][spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<3, dim> *third =
-                  &data.third_derivative(point + data_set, 0);
-
-                DerivativeForm<3, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            result[fe_to_real[comp_k]][j][l][m] +=
-                              (third[k][j][l][m] * data.local_dof_values[k]);
-                  }
-
-                // push forward the j-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        {
-                          jacobian_pushed_forward_2nd_derivatives
-                            [point][i][j][l][m] =
-                              result[i][0][l][m] * data.covariant[point][j][0];
-                          for (unsigned int jr = 1; jr < dim; ++jr)
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [l][m] +=
-                              result[i][jr][l][m] *
-                              data.covariant[point][j][jr];
-                        }
-
-                // push forward the l-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        {
-                          tmp[i][j][l][m] =
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [0][m] *
-                            data.covariant[point][l][0];
-                          for (unsigned int lr = 1; lr < dim; ++lr)
-                            tmp[i][j][l][m] +=
-                              jacobian_pushed_forward_2nd_derivatives[point][i]
-                                                                     [j][lr]
-                                                                     [m] *
-                              data.covariant[point][l][lr];
-                        }
-
-                // push forward the m-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < spacedim; ++m)
-                        {
-                          jacobian_pushed_forward_2nd_derivatives
-                            [point][i][j][l][m] =
-                              tmp[i][j][l][0] * data.covariant[point][m][0];
-                          for (unsigned int mr = 1; mr < dim; ++mr)
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [l][m] +=
-                              tmp[i][j][l][mr] * data.covariant[point][m][mr];
-                        }
-              }
-          }
-      }
-
-      /**
-       * Update the fourth derivative of the transformation from unit to real
-       * cell, the Jacobian hessian gradients.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_3rd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_3rd_derivatives)
-          {
-            const unsigned int n_q_points = jacobian_3rd_derivatives.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<4, dim> *fourth =
-                  &data.fourth_derivative(point + data_set, 0);
-
-                DerivativeForm<4, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.local_dof_values[k]);
-                  }
-
-                // never touch any data for j,l,m,n=dim in case
-                // dim<spacedim, so it will always be zero as it was
-                // initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          jacobian_3rd_derivatives[point][i][j][l][m][n] =
-                            result[i][j][l][m][n];
-              }
-          }
-      }
-
-      /**
-       * Update the fourth derivative of the transformation from unit to real
-       * cell, the Jacobian hessian gradients, pushed forward to the real cell
-       * coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_3rd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<5, spacedim>>
-          &jacobian_pushed_forward_3rd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_3rd_derivatives.size();
-
-            double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<4, dim> *fourth =
-                  &data.fourth_derivative(point + data_set, 0);
-
-                DerivativeForm<4, dim, spacedim> result;
-
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.local_dof_values[k]);
-                  }
-
-                // push-forward the j-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            tmp[i][j][l][m][n] = result[i][0][l][m][n] *
-                                                 data.covariant[point][j][0];
-                            for (unsigned int jr = 1; jr < dim; ++jr)
-                              tmp[i][j][l][m][n] +=
-                                result[i][jr][l][m][n] *
-                                data.covariant[point][j][jr];
-                          }
-
-                // push-forward the l-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives
-                              [point][i][j][l][m][n] =
-                                tmp[i][j][0][m][n] *
-                                data.covariant[point][l][0];
-                            for (unsigned int lr = 1; lr < dim; ++lr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][m]
-                                                                     [n] +=
-                                tmp[i][j][lr][m][n] *
-                                data.covariant[point][l][lr];
-                          }
-
-                // push-forward the m-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < spacedim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            tmp[i][j][l][m][n] =
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][0]
-                                                                     [n] *
-                              data.covariant[point][m][0];
-                            for (unsigned int mr = 1; mr < dim; ++mr)
-                              tmp[i][j][l][m][n] +=
-                                jacobian_pushed_forward_3rd_derivatives[point]
-                                                                       [i][j][l]
-                                                                       [mr][n] *
-                                data.covariant[point][m][mr];
-                          }
-
-                // push-forward the n-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < spacedim; ++m)
-                        for (unsigned int n = 0; n < spacedim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives
-                              [point][i][j][l][m][n] =
-                                tmp[i][j][l][m][0] *
-                                data.covariant[point][n][0];
-                            for (unsigned int nr = 1; nr < dim; ++nr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][m]
-                                                                     [n] +=
-                                tmp[i][j][l][m][nr] *
-                                data.covariant[point][n][nr];
-                          }
-              }
-          }
-      }
-
-
-      /**
-       * Depending on what information is called for in the update flags of the
-       * @p data object, compute the various pieces of information that is
-       * required by the fill_fe_face_values() and fill_fe_subface_values()
-       * functions.  This function simply unifies the work that would be done by
-       * those two functions.
-       *
-       * The resulting data is put into the @p output_data argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_compute_face_data(
-        const dealii::Mapping<dim, spacedim> &mapping,
-        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
-                                                         &cell,
-        const unsigned int                                face_no,
-        const unsigned int                                subface_no,
-        const typename QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data,
-        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-          &output_data)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        if (update_flags & update_boundary_forms)
-          {
-            const unsigned int n_q_points = output_data.boundary_forms.size();
-            if (update_flags & update_normal_vectors)
-              AssertDimension(output_data.normal_vectors.size(), n_q_points);
-            if (update_flags & update_JxW_values)
-              AssertDimension(output_data.JxW_values.size(), n_q_points);
-
-            // map the unit tangentials to the real cell. checking for d!=dim-1
-            // eliminates compiler warnings regarding unsigned int expressions <
-            // 0.
-            for (unsigned int d = 0; d != dim - 1; ++d)
-              {
-                Assert(face_no + cell->n_faces() * d <
-                         data.unit_tangentials.size(),
-                       ExcInternalError());
-                Assert(
-                  data.aux[d].size() <=
-                    data.unit_tangentials[face_no + cell->n_faces() * d].size(),
-                  ExcInternalError());
-
-                mapping.transform(
-                  make_array_view(
-                    data.unit_tangentials[face_no + cell->n_faces() * d]),
-                  mapping_contravariant,
-                  data,
-                  make_array_view(data.aux[d]));
-              }
-
-            // if dim==spacedim, we can use the unit tangentials to compute the
-            // boundary form by simply taking the cross product
-            if (dim == spacedim)
-              {
-                for (unsigned int i = 0; i < n_q_points; ++i)
-                  switch (dim)
-                    {
-                      case 1:
-                        // in 1d, we don't have access to any of the data.aux
-                        // fields (because it has only dim-1 components), but we
-                        // can still compute the boundary form by simply looking
-                        // at the number of the face
-                        output_data.boundary_forms[i][0] =
-                          (face_no == 0 ? -1 : +1);
-                        break;
-                      case 2:
-                        output_data.boundary_forms[i] =
-                          cross_product_2d(data.aux[0][i]);
-                        break;
-                      case 3:
-                        output_data.boundary_forms[i] =
-                          cross_product_3d(data.aux[0][i], data.aux[1][i]);
-                        break;
-                      default:
-                        DEAL_II_NOT_IMPLEMENTED();
-                    }
-              }
-            else //(dim < spacedim)
-              {
-                // in the codim-one case, the boundary form results from the
-                // cross product of all the face tangential vectors and the cell
-                // normal vector
-                //
-                // to compute the cell normal, use the same method used in
-                // fill_fe_values for cells above
-                AssertDimension(data.contravariant.size(), n_q_points);
-
-                for (unsigned int point = 0; point < n_q_points; ++point)
-                  {
-                    if (dim == 1)
-                      {
-                        // J is a tangent vector
-                        output_data.boundary_forms[point] =
-                          data.contravariant[point].transpose()[0];
-                        output_data.boundary_forms[point] /=
-                          (face_no == 0 ? -1. : +1.) *
-                          output_data.boundary_forms[point].norm();
-                      }
-
-                    if (dim == 2)
-                      {
-                        const DerivativeForm<1, spacedim, dim> DX_t =
-                          data.contravariant[point].transpose();
-
-                        Tensor<1, spacedim> cell_normal =
-                          cross_product_3d(DX_t[0], DX_t[1]);
-                        cell_normal /= cell_normal.norm();
-
-                        // then compute the face normal from the face tangent
-                        // and the cell normal:
-                        output_data.boundary_forms[point] =
-                          cross_product_3d(data.aux[0][point], cell_normal);
-                      }
-                  }
-              }
-
-            if (update_flags & (update_normal_vectors | update_JxW_values))
-              for (unsigned int i = 0; i < output_data.boundary_forms.size();
-                   ++i)
-                {
-                  if (update_flags & update_JxW_values)
-                    {
-                      output_data.JxW_values[i] =
-                        output_data.boundary_forms[i].norm() *
-                        data.quadrature_weights[i + data_set];
-
-                      if (subface_no != numbers::invalid_unsigned_int)
-                        {
-                          // TODO
-                          const double area_ratio =
-                            GeometryInfo<dim>::subface_ratio(
-                              cell->subface_case(face_no), subface_no);
-                          output_data.JxW_values[i] *= area_ratio;
-                        }
-                    }
-
-                  if (update_flags & update_normal_vectors)
-                    output_data.normal_vectors[i] =
-                      Point<spacedim>(output_data.boundary_forms[i] /
-                                      output_data.boundary_forms[i].norm());
-                }
-          }
-      }
-
-      /**
-       * Do the work of MappingFEField::fill_fe_face_values() and
-       * MappingFEField::fill_fe_subface_values() in a generic way, using the
-       * 'data_set' to differentiate whether we will work on a face (and if so,
-       * which one) or subface.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      do_fill_fe_face_values(
-        const dealii::Mapping<dim, spacedim> &mapping,
-        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
-                                                                 &cell,
-        const unsigned int                                        face_no,
-        const unsigned int                                        subface_no,
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-          &output_data)
-      {
-        maybe_compute_q_points<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.quadrature_points);
-
-        maybe_update_Jacobians<dim, spacedim, VectorType>(
-          data_set, data, fe, fe_mask, fe_to_real);
-
-        const UpdateFlags  update_flags = data.update_each;
-        const unsigned int n_q_points   = data.contravariant.size();
-
-        if (update_flags & update_jacobians)
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.jacobians[point] = data.contravariant[point];
-
-        if (update_flags & update_inverse_jacobians)
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.inverse_jacobians[point] =
-              data.covariant[point].transpose();
-
-        maybe_update_jacobian_grads<dim, spacedim, VectorType>(
-          data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
-
-        maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_grads);
-
-        maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_2nd_derivatives);
-
-        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_2nd_derivatives);
-
-        maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_3rd_derivatives);
-
-        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_3rd_derivatives);
-
-        maybe_compute_face_data<dim, spacedim, VectorType>(
-          mapping, cell, face_no, subface_no, data_set, data, output_data);
-      }
-    } // namespace
-  }   // namespace MappingFEFieldImplementation
-} // namespace internal
-
-
-// Note that the CellSimilarity flag is modifiable, since MappingFEField can
-// need to recalculate data even when cells are similar.
-template <int dim, int spacedim, typename VectorType>
-CellSimilarity::Similarity
-MappingFEField<dim, spacedim, VectorType>::fill_fe_values(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const CellSimilarity::Similarity,
-  const Quadrature<dim>                                   &quadrature,
-  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
-  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-    &output_data) const
-{
-  // convert data object to internal data for this class. fails with an
-  // exception if that is not possible
-  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
-         ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(internal_data);
-
-  const unsigned int n_q_points = quadrature.size();
-
-  update_internal_dofs(cell, data);
-
-  internal::MappingFEFieldImplementation::
-    maybe_compute_q_points<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.quadrature_points);
-
-  internal::MappingFEFieldImplementation::
-    maybe_update_Jacobians<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real);
-
-  const UpdateFlags          update_flags = data.update_each;
-  const std::vector<double> &weights      = quadrature.get_weights();
-
-  // Multiply quadrature weights by absolute value of Jacobian determinants or
-  // the area element g=sqrt(DX^t DX) in case of codim > 0
-
-  if (update_flags & (update_normal_vectors | update_JxW_values))
-    {
-      AssertDimension(output_data.JxW_values.size(), n_q_points);
-
-      Assert(!(update_flags & update_normal_vectors) ||
-               (output_data.normal_vectors.size() == n_q_points),
-             ExcDimensionMismatch(output_data.normal_vectors.size(),
-                                  n_q_points));
-
-
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        {
-          if (dim == spacedim)
-            {
-              const double det = data.volume_elements[point];
-
-              // check for distorted cells.
-
-              // TODO: this allows for anisotropies of up to 1e6 in 3d and
-              // 1e12 in 2d. might want to find a finer
-              // (dimension-independent) criterion
-              Assert(det > 1e-12 * Utilities::fixed_power<dim>(
-                                     cell->diameter() / std::sqrt(double(dim))),
-                     (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
-                       cell->center(), det, point)));
-              output_data.JxW_values[point] = weights[point] * det;
-            }
-          // if dim==spacedim, then there is no cell normal to
-          // compute. since this is for FEValues (and not FEFaceValues),
-          // there are also no face normals to compute
-          else // codim>0 case
-            {
-              Tensor<1, spacedim> DX_t[dim];
-              for (unsigned int i = 0; i < spacedim; ++i)
-                for (unsigned int j = 0; j < dim; ++j)
-                  DX_t[j][i] = data.contravariant[point][i][j];
-
-              Tensor<2, dim> G; // First fundamental form
-              for (unsigned int i = 0; i < dim; ++i)
-                for (unsigned int j = 0; j < dim; ++j)
-                  G[i][j] = DX_t[i] * DX_t[j];
-
-              output_data.JxW_values[point] =
-                std::sqrt(determinant(G)) * weights[point];
-
-              if (update_flags & update_normal_vectors)
-                {
-                  Assert(spacedim - dim == 1,
-                         ExcMessage("There is no cell normal in codim 2."));
-
-                  if (dim == 1)
-                    output_data.normal_vectors[point] =
-                      cross_product_2d(-DX_t[0]);
-                  else
-                    {
-                      Assert(dim == 2, ExcInternalError());
-
-                      // dim-1==1 for the second argument, but this
-                      // avoids a compiler warning about array bounds:
-                      output_data.normal_vectors[point] =
-                        cross_product_3d(DX_t[0], DX_t[dim - 1]);
-                    }
-
-                  output_data.normal_vectors[point] /=
-                    output_data.normal_vectors[point].norm();
-
-                  if (cell->direction_flag() == false)
-                    output_data.normal_vectors[point] *= -1.;
-                }
-            } // codim>0 case
-        }
-    }
-
-  // copy values from InternalData to vector given by reference
-  if (update_flags & update_jacobians)
-    {
-      AssertDimension(output_data.jacobians.size(), n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        output_data.jacobians[point] = data.contravariant[point];
-    }
-
-  // copy values from InternalData to vector given by reference
-  if (update_flags & update_inverse_jacobians)
-    {
-      AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        output_data.inverse_jacobians[point] =
-          data.covariant[point].transpose();
-    }
-
-  // calculate derivatives of the Jacobians
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_grads<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_grads);
-
-  // calculate derivatives of the Jacobians pushed forward to real cell
-  // coordinates
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_pushed_forward_grads);
-
-  // calculate hessians of the Jacobians
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_2nd_derivatives);
-
-  // calculate hessians of the Jacobians pushed forward to real cell coordinates
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
-                                                         spacedim,
-                                                         VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_pushed_forward_2nd_derivatives);
-
-  // calculate gradients of the hessians of the Jacobians
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_3rd_derivatives);
-
-  // calculate gradients of the hessians of the Jacobians pushed forward to real
-  // cell coordinates
-  internal::MappingFEFieldImplementation::
-    maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
-                                                         spacedim,
-                                                         VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.jacobian_pushed_forward_3rd_derivatives);
-
-  return CellSimilarity::invalid_next_cell;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::fill_fe_face_values(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const unsigned int                                          face_no,
-  const hp::QCollection<dim - 1>                             &quadrature,
-  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
-  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-    &output_data) const
-{
-  AssertDimension(quadrature.size(), 1);
-
-  // convert data object to internal data for this class. fails with an
-  // exception if that is not possible
-  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
-         ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(internal_data);
-
-  update_internal_dofs(cell, data);
-
-  internal::MappingFEFieldImplementation::
-    do_fill_fe_face_values<dim, spacedim, VectorType>(
-      *this,
-      cell,
-      face_no,
-      numbers::invalid_unsigned_int,
-      QProjector<dim>::DataSetDescriptor::face(reference_cell,
-                                               face_no,
-                                               cell->combined_face_orientation(
-                                                 face_no),
-                                               quadrature[0].size()),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data);
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::fill_fe_subface_values(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const unsigned int                                          face_no,
-  const unsigned int                                          subface_no,
-  const Quadrature<dim - 1>                                  &quadrature,
-  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
-  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-    &output_data) const
-{
-  // convert data object to internal data for this class. fails with an
-  // exception if that is not possible
-  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
-         ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(internal_data);
-
-  update_internal_dofs(cell, data);
-
-  internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
-                                                                 spacedim,
-                                                                 VectorType>(
-    *this,
-    cell,
-    face_no,
-    numbers::invalid_unsigned_int,
-    QProjector<dim>::DataSetDescriptor::subface(reference_cell,
-                                                face_no,
-                                                subface_no,
-                                                cell->combined_face_orientation(
-                                                  face_no),
-                                                quadrature.size(),
-                                                cell->subface_case(face_no)),
-    data,
-    euler_dof_handler->get_fe(),
-    fe_mask,
-    fe_to_real,
-    output_data);
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::fill_fe_immersed_surface_values(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const NonMatching::ImmersedSurfaceQuadrature<dim>          &quadrature,
-  const typename Mapping<dim, spacedim>::InternalDataBase    &internal_data,
-  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-    &output_data) const
-{
-  AssertDimension(dim, spacedim);
-  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
-         ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(internal_data);
-
-  const unsigned int n_q_points = quadrature.size();
-
-  update_internal_dofs(cell, data);
-
-  internal::MappingFEFieldImplementation::
-    maybe_compute_q_points<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real,
-      output_data.quadrature_points);
-
-  internal::MappingFEFieldImplementation::
-    maybe_update_Jacobians<dim, spacedim, VectorType>(
-      QProjector<dim>::DataSetDescriptor::cell(),
-      data,
-      euler_dof_handler->get_fe(),
-      fe_mask,
-      fe_to_real);
-
-  const UpdateFlags          update_flags = data.update_each;
-  const std::vector<double> &weights      = quadrature.get_weights();
-
-  if (update_flags & (update_normal_vectors | update_JxW_values))
-    {
-      AssertDimension(output_data.JxW_values.size(), n_q_points);
-
-      Assert(!(update_flags & update_normal_vectors) ||
-               (output_data.normal_vectors.size() == n_q_points),
-             ExcDimensionMismatch(output_data.normal_vectors.size(),
-                                  n_q_points));
-
-
-      for (unsigned int point = 0; point < n_q_points; ++point)
-        {
-          const double det = data.volume_elements[point];
-
-          // check for distorted cells.
-
-          // TODO: this allows for anisotropies of up to 1e6 in 3d and
-          // 1e12 in 2d. might want to find a finer
-          // (dimension-independent) criterion
-          Assert(det > 1e-12 * Utilities::fixed_power<dim>(
-                                 cell->diameter() / std::sqrt(double(dim))),
-                 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
-                   cell->center(), det, point)));
-
-          // The normals are n = J^{-T} * \hat{n} before normalizing.
-          Tensor<1, spacedim> normal;
-          for (unsigned int d = 0; d < spacedim; d++)
-            normal[d] =
-              data.covariant[point][d] * quadrature.normal_vector(point);
-
-          output_data.JxW_values[point] = weights[point] * det * normal.norm();
-
-          if ((update_flags & update_normal_vectors) != 0u)
-            {
-              normal /= normal.norm();
-              output_data.normal_vectors[point] = normal;
-            }
-        }
-
-      // copy values from InternalData to vector given by reference
-      if (update_flags & update_jacobians)
-        {
-          AssertDimension(output_data.jacobians.size(), n_q_points);
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.jacobians[point] = data.contravariant[point];
-        }
-
-      // copy values from InternalData to vector given by reference
-      if (update_flags & update_inverse_jacobians)
-        {
-          AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.inverse_jacobians[point] =
-              data.covariant[point].transpose();
-        }
-
-      // calculate derivatives of the Jacobians
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_grads<dim, spacedim, VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_grads);
-
-      // calculate derivatives of the Jacobians pushed forward to real cell
-      // coordinates
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_grads);
-
-      // calculate hessians of the Jacobians
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_2nd_derivatives);
-
-      // calculate hessians of the Jacobians pushed forward to real cell
-      // coordinates
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_2nd_derivatives);
-
-      // calculate gradients of the hessians of the Jacobians
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_3rd_derivatives);
-
-      // calculate gradients of the hessians of the Jacobians pushed forward to
-      // real cell coordinates
-      internal::MappingFEFieldImplementation::
-        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          QProjector<dim>::DataSetDescriptor::cell(),
-          data,
-          euler_dof_handler->get_fe(),
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_3rd_derivatives);
-    }
-}
-
-namespace internal
-{
-  namespace MappingFEFieldImplementation
-  {
-    namespace
-    {
-      template <int dim, int spacedim, int rank, typename VectorType>
-      void
-      transform_fields(
-        const ArrayView<const Tensor<rank, dim>>                &input,
-        const MappingKind                                        mapping_kind,
-        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-        const ArrayView<Tensor<rank, spacedim>>                 &output)
-      {
-        AssertDimension(input.size(), output.size());
-        Assert((dynamic_cast<
-                  const typename dealii::
-                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
-                  &mapping_data) != nullptr),
-               ExcInternalError());
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data = static_cast<
-            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-              InternalData &>(mapping_data);
-
-        switch (mapping_kind)
-          {
-            case mapping_contravariant:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
-
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  output[i] =
-                    apply_transformation(data.contravariant[i], input[i]);
-
-                return;
-              }
-
-            case mapping_piola:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
-                Assert(
-                  data.update_each & update_volume_elements,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_volume_elements"));
-                Assert(rank == 1, ExcMessage("Only for rank 1"));
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  {
-                    output[i] =
-                      apply_transformation(data.contravariant[i], input[i]);
-                    output[i] /= data.volume_elements[i];
-                  }
-                return;
-              }
-
-
-            // We still allow this operation as in the
-            // reference cell Derivatives are Tensor
-            // rather than DerivativeForm
-            case mapping_covariant:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
-
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  output[i] = apply_transformation(data.covariant[i], input[i]);
-
-                return;
-              }
-
-            default:
-              DEAL_II_NOT_IMPLEMENTED();
-          }
-      }
-
-
-      template <int dim, int spacedim, int rank, typename VectorType>
-      void
-      transform_differential_forms(
-        const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
-        const MappingKind                                        mapping_kind,
-        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-        const ArrayView<Tensor<rank + 1, spacedim>>             &output)
-      {
-        AssertDimension(input.size(), output.size());
-        Assert((dynamic_cast<
-                  const typename dealii::
-                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
-                  &mapping_data) != nullptr),
-               ExcInternalError());
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data = static_cast<
-            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-              InternalData &>(mapping_data);
-
-        switch (mapping_kind)
-          {
-            case mapping_covariant:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
-
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  output[i] = apply_transformation(data.covariant[i], input[i]);
-
-                return;
-              }
-            default:
-              DEAL_II_NOT_IMPLEMENTED();
-          }
-      }
-    } // namespace
-  }   // namespace MappingFEFieldImplementation
-} // namespace internal
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::transform(
-  const ArrayView<const Tensor<1, dim>>                   &input,
-  const MappingKind                                        mapping_kind,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-  const ArrayView<Tensor<1, spacedim>>                    &output) const
-{
-  AssertDimension(input.size(), output.size());
-
-  internal::MappingFEFieldImplementation::
-    transform_fields<dim, spacedim, 1, VectorType>(input,
-                                                   mapping_kind,
-                                                   mapping_data,
-                                                   output);
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::transform(
-  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
-  const MappingKind                                        mapping_kind,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-  const ArrayView<Tensor<2, spacedim>>                    &output) const
-{
-  AssertDimension(input.size(), output.size());
-
-  internal::MappingFEFieldImplementation::
-    transform_differential_forms<dim, spacedim, 1, VectorType>(input,
-                                                               mapping_kind,
-                                                               mapping_data,
-                                                               output);
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::transform(
-  const ArrayView<const Tensor<2, dim>> &input,
-  const MappingKind,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-  const ArrayView<Tensor<2, spacedim>>                    &output) const
-{
-  (void)input;
-  (void)output;
-  (void)mapping_data;
-  AssertDimension(input.size(), output.size());
-
-  AssertThrow(false, ExcNotImplemented());
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::transform(
-  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
-  const MappingKind                                        mapping_kind,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-  const ArrayView<Tensor<3, spacedim>>                    &output) const
-{
-  AssertDimension(input.size(), output.size());
-  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
-         ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(mapping_data);
-
-  switch (mapping_kind)
-    {
-      case mapping_covariant_gradient:
-        {
-          Assert(data.update_each & update_contravariant_transformation,
-                 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                   "update_covariant_transformation"));
-
-          for (unsigned int q = 0; q < output.size(); ++q)
-            for (unsigned int i = 0; i < spacedim; ++i)
-              for (unsigned int j = 0; j < spacedim; ++j)
-                for (unsigned int k = 0; k < spacedim; ++k)
-                  {
-                    output[q][i][j][k] = data.covariant[q][j][0] *
-                                         data.covariant[q][k][0] *
-                                         input[q][i][0][0];
-                    for (unsigned int J = 0; J < dim; ++J)
-                      {
-                        const unsigned int K0 = (0 == J) ? 1 : 0;
-                        for (unsigned int K = K0; K < dim; ++K)
-                          output[q][i][j][k] += data.covariant[q][j][J] *
-                                                data.covariant[q][k][K] *
-                                                input[q][i][J][K];
-                      }
-                  }
-          return;
-        }
-
-      default:
-        DEAL_II_NOT_IMPLEMENTED();
-    }
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::transform(
-  const ArrayView<const Tensor<3, dim>> &input,
-  const MappingKind /*mapping_kind*/,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-  const ArrayView<Tensor<3, spacedim>>                    &output) const
-{
-  (void)input;
-  (void)output;
-  (void)mapping_data;
-  AssertDimension(input.size(), output.size());
-
-  AssertThrow(false, ExcNotImplemented());
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Point<spacedim>
-MappingFEField<dim, spacedim, VectorType>::transform_unit_to_real_cell(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const Point<dim>                                           &p) const
-{
-  //  Use the get_data function to create an InternalData with data vectors of
-  //  the right size and transformation shape values already computed at point
-  //  p.
-  const Quadrature<dim> point_quadrature(p);
-  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
-    get_data(update_quadrature_points | update_jacobians, point_quadrature));
-  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
-         ExcInternalError());
-
-  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
-
-  return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-Point<spacedim>
-MappingFEField<dim, spacedim, VectorType>::do_transform_unit_to_real_cell(
-  const InternalData &data) const
-{
-  Point<spacedim> p_real;
-
-  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-    {
-      unsigned int comp_i =
-        euler_dof_handler->get_fe().system_to_component_index(i).first;
-      if (fe_mask[comp_i])
-        p_real[fe_to_real[comp_i]] +=
-          data.local_dof_values[i] * data.shape(0, i);
-    }
-
-  return p_real;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-Point<dim>
-MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const Point<spacedim>                                      &p) const
-{
-  // first a Newton iteration based on the real mapping. It uses the center
-  // point of the cell as a starting point
-  Point<dim> initial_p_unit;
-  try
-    {
-      initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
-                         .transform_real_to_unit_cell(cell, p);
-    }
-  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
-    {
-      // mirror the conditions of the code below to determine if we need to
-      // use an arbitrary starting point or if we just need to rethrow the
-      // exception
-      for (unsigned int d = 0; d < dim; ++d)
-        initial_p_unit[d] = 0.5;
-    }
-
-  initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
-
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim > dim)
-    update_flags |= update_jacobian_grads;
-  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
-    get_data(update_flags, Quadrature<dim>(initial_p_unit)));
-  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
-         ExcInternalError());
-
-  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
-
-  return do_transform_real_to_unit_cell(cell,
-                                        p,
-                                        initial_p_unit,
-                                        static_cast<InternalData &>(*mdata));
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-Point<dim>
-MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const Point<spacedim>                                      &p,
-  const Point<dim>                                           &starting_guess,
-  InternalData                                               &mdata) const
-{
-  const unsigned int n_shapes = mdata.shape_values.size();
-  (void)n_shapes;
-  Assert(n_shapes != 0, ExcInternalError());
-  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
-
-
-  // Newton iteration to solve
-  // f(x)=p(x)-p=0
-  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
-  // The start value was set to be the
-  // linear approximation to the cell
-  // The shape values and derivatives
-  // of the mapping at this point are
-  // previously computed.
-
-  Point<dim> p_unit = starting_guess;
-  Point<dim> f;
-  mdata.reinit(mdata.update_each, Quadrature<dim>(starting_guess));
-
-  Point<spacedim>     p_real(do_transform_unit_to_real_cell(mdata));
-  Tensor<1, spacedim> p_minus_F              = p - p_real;
-  const double        eps                    = 1.e-12 * cell->diameter();
-  const unsigned int  newton_iteration_limit = 20;
-  unsigned int        newton_iteration       = 0;
-  while (p_minus_F.norm_square() > eps * eps)
-    {
-      // f'(x)
-      Point<spacedim> DF[dim];
-      Tensor<2, dim>  df;
-      for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
-        {
-          const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
-          const unsigned int    comp_k =
-            euler_dof_handler->get_fe().system_to_component_index(k).first;
-          if (fe_mask[comp_k])
-            for (unsigned int j = 0; j < dim; ++j)
-              DF[j][fe_to_real[comp_k]] +=
-                mdata.local_dof_values[k] * grad_k[j];
-        }
-      for (unsigned int j = 0; j < dim; ++j)
-        {
-          f[j] = DF[j] * p_minus_F;
-          for (unsigned int l = 0; l < dim; ++l)
-            df[j][l] = -DF[j] * DF[l];
-        }
-      // Solve  [f'(x)]d=f(x)
-      const Tensor<1, dim> delta =
-        invert(df) * static_cast<const Tensor<1, dim> &>(f);
-      // do a line search
-      double step_length = 1;
-      do
-        {
-          // update of p_unit. The
-          // spacedimth component of
-          // transformed point is simply
-          // ignored in codimension one
-          // case. When this component is
-          // not zero, then we are
-          // projecting the point to the
-          // surface or curve identified
-          // by the cell.
-          Point<dim> p_unit_trial = p_unit;
-          for (unsigned int i = 0; i < dim; ++i)
-            p_unit_trial[i] -= step_length * delta[i];
-          // shape values and derivatives
-          // at new p_unit point
-          mdata.reinit(mdata.update_each, Quadrature<dim>(p_unit_trial));
-          // f(x)
-          const Point<spacedim> p_real_trial =
-            do_transform_unit_to_real_cell(mdata);
-          const Tensor<1, spacedim> f_trial = p - p_real_trial;
-          // see if we are making progress with the current step length
-          // and if not, reduce it by a factor of two and try again
-          if (f_trial.norm() < p_minus_F.norm())
-            {
-              p_real    = p_real_trial;
-              p_unit    = p_unit_trial;
-              p_minus_F = f_trial;
-              break;
-            }
-          else if (step_length > 0.05)
-            step_length /= 2;
-          else
-            goto failure;
-        }
-      while (true);
-      ++newton_iteration;
-      if (newton_iteration > newton_iteration_limit)
-        goto failure;
-    }
-  return p_unit;
-  // if we get to the following label, then we have either run out
-  // of Newton iterations, or the line search has not converged.
-  // in either case, we need to give up, so throw an exception that
-  // can then be caught
-failure:
-  AssertThrow(false,
-              (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
-  // ...the compiler wants us to return something, though we can
-  // of course never get here...
-  return {};
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-unsigned int
-MappingFEField<dim, spacedim, VectorType>::get_degree() const
-{
-  return euler_dof_handler->get_fe().degree;
-}
-
-
-
-template <int dim, int spacedim, typename VectorType>
-ComponentMask
-MappingFEField<dim, spacedim, VectorType>::get_component_mask() const
-{
-  return this->fe_mask;
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-std::unique_ptr<Mapping<dim, spacedim>>
-MappingFEField<dim, spacedim, VectorType>::clone() const
-{
-  return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::update_internal_dofs(
-  const typename Triangulation<dim, spacedim>::cell_iterator             &cell,
-  const typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
-  const
-{
-  Assert(euler_dof_handler != nullptr,
-         ExcMessage("euler_dof_handler is empty"));
-
-  typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
-                                                             euler_dof_handler);
-  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
-  if (uses_level_dofs)
-    {
-      AssertIndexRange(cell->level(), euler_vector.size());
-      AssertDimension(euler_vector[cell->level()]->size(),
-                      euler_dof_handler->n_dofs(cell->level()));
-    }
-  else
-    AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
-
-  if (uses_level_dofs)
-    dof_cell->get_mg_dof_indices(data.local_dof_indices);
-  else
-    dof_cell->get_dof_indices(data.local_dof_indices);
-
-  const VectorType &vector =
-    uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
-
-  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
-    data.local_dof_values[i] =
-      internal::ElementAccess<VectorType>::get(vector,
-                                               data.local_dof_indices[i]);
-}
-
 // explicit instantiations
 #define SPLIT_INSTANTIATIONS_COUNT 2
 #ifndef SPLIT_INSTANTIATIONS_INDEX

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.