]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Intoduce MappingFE 10659/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 3 Jul 2020 14:44:55 +0000 (16:44 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 27 Jul 2020 13:03:45 +0000 (15:03 +0200)
include/deal.II/fe/mapping_fe.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/mapping_fe.cc [new file with mode: 0644]
source/fe/mapping_fe.inst.in [new file with mode: 0644]

diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h
new file mode 100644 (file)
index 0000000..78dee1d
--- /dev/null
@@ -0,0 +1,604 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_mapping_fe_h
+#define dealii_mapping_fe_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/tria_iterator.h>
+
+#include <array>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*!@addtogroup mapping */
+/*@{*/
+
+
+/**
+ * This class consistently uses a user-provided finite element on all cells of a
+ * triangulation to implement a polynomial mapping.
+ *
+ * If one initializes this class with the same FiniteElement as the
+ * discretization, one obtains an iso-parametric mapping.
+ *
+ * If one initializes this class with an FE_Q(degree) object, then this class is
+ * equivalent to MappingQGeneric(degree). Please note that no optimizations
+ * exploiting tensor-product structures of finite elements have been added here.
+ *
+ * @node Currently, only implemented for degree==1.
+ *
+ * @ingroupalso simplex
+ */
+template <int dim, int spacedim = dim>
+class MappingFE : public Mapping<dim, spacedim>
+{
+public:
+  /**
+   * Constructor.
+   */
+  explicit MappingFE(const FiniteElement<dim, spacedim> &fe);
+
+  /**
+   * Copy constructor.
+   */
+  MappingFE(const MappingFE<dim, spacedim> &mapping);
+
+  // for documentation, see the Mapping base class
+  virtual std::unique_ptr<Mapping<dim, spacedim>>
+  clone() const override;
+
+  /**
+   * Return the degree of the mapping, i.e., the degree of the finite element
+   * which was passed to the constructor.
+   */
+  unsigned int
+  get_degree() const;
+
+  /**
+   * Always returns @p true because the default implementation of functions in
+   * this class preserves vertex locations.
+   */
+  virtual bool
+  preserves_vertex_locations() const override;
+
+  /**
+   * @name Mapping points between reference and real cells
+   * @{
+   */
+
+  // for documentation, see the Mapping base class
+  virtual Point<spacedim>
+  transform_unit_to_real_cell(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const Point<dim> &p) const override;
+
+  /**
+   * for documentation, see the Mapping base class
+   *
+   * note: not implemented yet
+   */
+  virtual Point<dim>
+  transform_real_to_unit_cell(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const Point<spacedim> &p) const override;
+
+  /**
+   * @}
+   */
+
+  /**
+   * @name Functions to transform tensors from reference to real coordinates
+   * @{
+   */
+
+  // for documentation, see the Mapping base class
+  virtual void
+  transform(const ArrayView<const Tensor<1, dim>> &                  input,
+            const MappingKind                                        kind,
+            const typename Mapping<dim, spacedim>::InternalDataBase &internal,
+            const ArrayView<Tensor<1, spacedim>> &output) const override;
+
+  // for documentation, see the Mapping base class
+  virtual void
+  transform(const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
+            const MappingKind                                        kind,
+            const typename Mapping<dim, spacedim>::InternalDataBase &internal,
+            const ArrayView<Tensor<2, spacedim>> &output) const override;
+
+  // for documentation, see the Mapping base class
+  virtual void
+  transform(const ArrayView<const Tensor<2, dim>> &                  input,
+            const MappingKind                                        kind,
+            const typename Mapping<dim, spacedim>::InternalDataBase &internal,
+            const ArrayView<Tensor<2, spacedim>> &output) const override;
+
+  // for documentation, see the Mapping base class
+  virtual void
+  transform(const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
+            const MappingKind                                        kind,
+            const typename Mapping<dim, spacedim>::InternalDataBase &internal,
+            const ArrayView<Tensor<3, spacedim>> &output) const override;
+
+  // for documentation, see the Mapping base class
+  virtual void
+  transform(const ArrayView<const Tensor<3, dim>> &                  input,
+            const MappingKind                                        kind,
+            const typename Mapping<dim, spacedim>::InternalDataBase &internal,
+            const ArrayView<Tensor<3, spacedim>> &output) const override;
+
+  /**
+   * @}
+   */
+
+  /**
+   * @name Interface with FEValues
+   * @{
+   */
+
+  /**
+   * Storage for internal data of polynomial mappings. See
+   * Mapping::InternalDataBase for an extensive description.
+   *
+   * For the current class, the InternalData class stores data that is
+   * computed once when the object is created (in get_data()) as well as data
+   * the class wants to store from between the call to fill_fe_values(),
+   * fill_fe_face_values(), or fill_fe_subface_values() until possible later
+   * calls from the finite element to functions such as transform(). The
+   * latter class of member variables are marked as 'mutable'.
+   */
+  class InternalData : public Mapping<dim, spacedim>::InternalDataBase
+  {
+  public:
+    /**
+     * Constructor.
+     */
+    InternalData(const FiniteElement<dim, spacedim> &fe);
+
+    /**
+     * Initialize the object's member variables related to cell data based on
+     * the given arguments.
+     *
+     * The function also calls compute_shape_function_values() to actually set
+     * the member variables related to the values and derivatives of the
+     * mapping shape functions.
+     */
+    void
+    initialize(const UpdateFlags      update_flags,
+               const Quadrature<dim> &quadrature,
+               const unsigned int     n_original_q_points);
+
+    /**
+     * Initialize the object's member variables related to cell and face data
+     * based on the given arguments. In order to initialize cell data, this
+     * function calls initialize().
+     */
+    void
+    initialize_face(const UpdateFlags      update_flags,
+                    const Quadrature<dim> &quadrature,
+                    const unsigned int     n_original_q_points);
+
+    /**
+     * Compute the values and/or derivatives of the shape functions used for
+     * the mapping.
+     */
+    void
+    compute_shape_function_values(const std::vector<Point<dim>> &unit_points);
+
+
+    /**
+     * Shape function at quadrature point. Shape functions are in tensor
+     * product order, so vertices must be reordered to obtain transformation.
+     */
+    const double &
+    shape(const unsigned int qpoint, const unsigned int shape_nr) const;
+
+    /**
+     * Shape function at quadrature point. See above.
+     */
+    double &
+    shape(const unsigned int qpoint, const unsigned int shape_nr);
+
+    /**
+     * Gradient of shape function in quadrature point. See above.
+     */
+    const Tensor<1, dim> &
+    derivative(const unsigned int qpoint, const unsigned int shape_nr) const;
+
+    /**
+     * Gradient of shape function in quadrature point. See above.
+     */
+    Tensor<1, dim> &
+    derivative(const unsigned int qpoint, const unsigned int shape_nr);
+
+    /**
+     * Second derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<2, dim> &
+    second_derivative(const unsigned int qpoint,
+                      const unsigned int shape_nr) const;
+
+    /**
+     * Second derivative of shape function in quadrature point. See above.
+     */
+    Tensor<2, dim> &
+    second_derivative(const unsigned int qpoint, const unsigned int shape_nr);
+
+    /**
+     * third derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<3, dim> &
+    third_derivative(const unsigned int qpoint,
+                     const unsigned int shape_nr) const;
+
+    /**
+     * third derivative of shape function in quadrature point. See above.
+     */
+    Tensor<3, dim> &
+    third_derivative(const unsigned int qpoint, const unsigned int shape_nr);
+
+    /**
+     * fourth derivative of shape function in quadrature point. See above.
+     */
+    const Tensor<4, dim> &
+    fourth_derivative(const unsigned int qpoint,
+                      const unsigned int shape_nr) const;
+
+    /**
+     * fourth derivative of shape function in quadrature point. See above.
+     */
+    Tensor<4, dim> &
+    fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr);
+
+    /**
+     * Return an estimate (in bytes) for the memory consumption of this object.
+     */
+    virtual std::size_t
+    memory_consumption() const override;
+
+    /**
+     * Values of shape functions. Access by function @p shape.
+     *
+     * Computed once.
+     */
+    std::vector<double> shape_values;
+
+    /**
+     * Values of shape function derivatives. Access by function @p derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<1, dim>> shape_derivatives;
+
+    /**
+     * Values of shape function second derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<2, dim>> shape_second_derivatives;
+
+    /**
+     * Values of shape function third derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<3, dim>> shape_third_derivatives;
+
+    /**
+     * Values of shape function fourth derivatives. Access by function @p
+     * second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<4, dim>> shape_fourth_derivatives;
+
+    /**
+     * Unit tangential vectors. Used for the computation of boundary forms and
+     * normal vectors.
+     *
+     * Filled once.
+     */
+    std::array<std::vector<Tensor<1, dim>>,
+               GeometryInfo<dim>::faces_per_cell *(dim - 1)>
+      unit_tangentials;
+
+    /**
+     * Underlying finite element.
+     */
+    const FiniteElement<dim, spacedim> &fe;
+
+    /**
+     * The polynomial degree of the mapping.
+     */
+    const unsigned int polynomial_degree;
+
+    /**
+     * Number of shape functions.
+     */
+    const unsigned int n_shape_functions;
+
+    /**
+     * Tensors of covariant transformation at each of the quadrature points.
+     * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
+     * Jacobian, is the first fundamental form of the map; if dim=spacedim
+     * then it reduces to the transpose of the inverse of the Jacobian matrix,
+     * which itself is stored in the @p contravariant field of this structure.
+     *
+     * Computed on each cell.
+     */
+    mutable std::vector<DerivativeForm<1, dim, spacedim>> covariant;
+
+    /**
+     * Tensors of contravariant transformation at each of the quadrature
+     * points. The contravariant matrix is the Jacobian of the transformation,
+     * i.e. $J_{ij}=dx_i/d\hat x_j$.
+     *
+     * Computed on each cell.
+     */
+    mutable std::vector<DerivativeForm<1, dim, spacedim>> contravariant;
+
+    /**
+     * Auxiliary vectors for internal use.
+     */
+    mutable std::vector<std::vector<Tensor<1, spacedim>>> aux;
+
+    /**
+     * Stores the support points of the mapping shape functions on the @p
+     * cell_of_current_support_points.
+     */
+    mutable std::vector<Point<spacedim>> mapping_support_points;
+
+    /**
+     * Stores the cell of which the @p mapping_support_points are stored.
+     */
+    mutable typename Triangulation<dim, spacedim>::cell_iterator
+      cell_of_current_support_points;
+
+    /**
+     * The determinant of the Jacobian in each quadrature point. Filled if
+     * #update_volume_elements.
+     */
+    mutable std::vector<double> volume_elements;
+
+    /**
+     * Projected quadrature weights.
+     */
+    mutable std::vector<double> quadrature_weights;
+  };
+
+
+  // documentation can be found in Mapping::requires_update_flags()
+  virtual UpdateFlags
+  requires_update_flags(const UpdateFlags update_flags) const override;
+
+  // documentation can be found in Mapping::get_data()
+  virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+  get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+
+  // documentation can be found in Mapping::get_face_data()
+  virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+  get_face_data(const UpdateFlags          flags,
+                const Quadrature<dim - 1> &quadrature) const override;
+
+  // documentation can be found in Mapping::get_subface_data()
+  virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+  get_subface_data(const UpdateFlags          flags,
+                   const Quadrature<dim - 1> &quadrature) const override;
+
+  // documentation can be found in Mapping::fill_fe_values()
+  virtual CellSimilarity::Similarity
+  fill_fe_values(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const CellSimilarity::Similarity                            cell_similarity,
+    const Quadrature<dim> &                                     quadrature,
+    const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      &output_data) const override;
+
+  // documentation can be found in Mapping::fill_fe_face_values()
+  virtual void
+  fill_fe_face_values(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_no,
+    const Quadrature<dim - 1> &                                 quadrature,
+    const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      &output_data) const override;
+
+  // documentation can be found in Mapping::fill_fe_subface_values()
+  virtual void
+  fill_fe_subface_values(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_no,
+    const unsigned int                                          subface_no,
+    const Quadrature<dim - 1> &                                 quadrature,
+    const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      &output_data) const override;
+
+  /**
+   * @}
+   */
+
+protected:
+  const std::unique_ptr<FiniteElement<dim, spacedim>> fe;
+
+  /**
+   * The degree of the polynomials used as shape functions for the mapping of
+   * cells.
+   */
+  const unsigned int polynomial_degree;
+
+  /**
+   * Return the locations of support points for the mapping.
+   */
+  virtual std::vector<Point<spacedim>>
+  compute_mapping_support_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
+};
+
+
+
+/*@}*/
+
+/*----------------------------------------------------------------------*/
+
+#ifndef DOXYGEN
+
+template <int dim, int spacedim>
+inline const double &
+MappingFE<dim, spacedim>::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>
+inline double &
+MappingFE<dim, spacedim>::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>
+inline const Tensor<1, dim> &
+MappingFE<dim, spacedim>::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>
+inline Tensor<1, dim> &
+MappingFE<dim, spacedim>::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>
+inline const Tensor<2, dim> &
+MappingFE<dim, spacedim>::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>
+inline Tensor<2, dim> &
+MappingFE<dim, spacedim>::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>
+inline const Tensor<3, dim> &
+MappingFE<dim, spacedim>::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>
+inline Tensor<3, dim> &
+MappingFE<dim, spacedim>::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>
+inline const Tensor<4, dim> &
+MappingFE<dim, spacedim>::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>
+inline Tensor<4, dim> &
+MappingFE<dim, spacedim>::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>
+inline bool
+MappingFE<dim, spacedim>::preserves_vertex_locations() const
+{
+  return true;
+}
+
+#endif // DOXYGEN
+
+/* -------------- declaration of explicit specializations ------------- */
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 4eacd86c356c5993968396891876f74948388a18..666999b0a5b0985de05fe1526267904af725f143 100644 (file)
@@ -74,6 +74,7 @@ SET(_separate_src
   fe_tools.cc
   fe_tools_interpolate.cc
   fe_tools_extrapolate.cc
+  mapping_fe.cc
   mapping_q_generic.cc
   mapping_q1_eulerian.cc
   mapping_q_eulerian.cc
@@ -130,6 +131,7 @@ SET(_inst
   mapping_c1.inst.in
   mapping_cartesian.inst.in
   mapping.inst.in
+  mapping_fe.inst.in
   mapping_fe_field.inst.in
   mapping_q_generic.inst.in
   mapping_q1_eulerian.inst.in
diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc
new file mode 100644 (file)
index 0000000..e47f04d
--- /dev/null
@@ -0,0 +1,2210 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2000 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+
+#include <deal.II/fe/fe_poly.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/tensor_product_matrix.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <memory>
+#include <numeric>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim, int spacedim>
+MappingFE<dim, spacedim>::InternalData::InternalData(
+  const FiniteElement<dim, spacedim> &fe)
+  : fe(fe)
+  , polynomial_degree(fe.tensor_degree())
+  , n_shape_functions(fe.n_dofs_per_cell())
+{}
+
+
+
+template <int dim, int spacedim>
+std::size_t
+MappingFE<dim, spacedim>::InternalData::memory_consumption() const
+{
+  return (
+    Mapping<dim, spacedim>::InternalDataBase::memory_consumption() +
+    MemoryConsumption::memory_consumption(shape_values) +
+    MemoryConsumption::memory_consumption(shape_derivatives) +
+    MemoryConsumption::memory_consumption(covariant) +
+    MemoryConsumption::memory_consumption(contravariant) +
+    MemoryConsumption::memory_consumption(unit_tangentials) +
+    MemoryConsumption::memory_consumption(aux) +
+    MemoryConsumption::memory_consumption(mapping_support_points) +
+    MemoryConsumption::memory_consumption(cell_of_current_support_points) +
+    MemoryConsumption::memory_consumption(volume_elements) +
+    MemoryConsumption::memory_consumption(polynomial_degree) +
+    MemoryConsumption::memory_consumption(n_shape_functions));
+}
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::InternalData::initialize(
+  const UpdateFlags      update_flags,
+  const Quadrature<dim> &q,
+  const unsigned int     n_original_q_points)
+{
+  // store the flags in the internal data object so we can access them
+  // in fill_fe_*_values()
+  this->update_each = update_flags;
+
+  const unsigned int n_q_points = q.size();
+
+  if (this->update_each & update_covariant_transformation)
+    covariant.resize(n_original_q_points);
+
+  if (this->update_each & update_contravariant_transformation)
+    contravariant.resize(n_original_q_points);
+
+  if (this->update_each & update_volume_elements)
+    volume_elements.resize(n_original_q_points);
+
+  // see if we need the (transformation) shape function values
+  // and/or gradients and resize the necessary arrays
+  if (this->update_each & update_quadrature_points)
+    shape_values.resize(n_shape_functions * n_q_points);
+
+  if (this->update_each &
+      (update_covariant_transformation | update_contravariant_transformation |
+       update_JxW_values | update_boundary_forms | update_normal_vectors |
+       update_jacobians | update_jacobian_grads | update_inverse_jacobians |
+       update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
+       update_jacobian_pushed_forward_2nd_derivatives |
+       update_jacobian_3rd_derivatives |
+       update_jacobian_pushed_forward_3rd_derivatives))
+    shape_derivatives.resize(n_shape_functions * n_q_points);
+
+  if (this->update_each &
+      (update_jacobian_grads | update_jacobian_pushed_forward_grads))
+    shape_second_derivatives.resize(n_shape_functions * n_q_points);
+
+  if (this->update_each & (update_jacobian_2nd_derivatives |
+                           update_jacobian_pushed_forward_2nd_derivatives))
+    shape_third_derivatives.resize(n_shape_functions * n_q_points);
+
+  if (this->update_each & (update_jacobian_3rd_derivatives |
+                           update_jacobian_pushed_forward_3rd_derivatives))
+    shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
+
+  // now also fill the various fields with their correct values
+  compute_shape_function_values(q.get_points());
+
+  // copy (projected) quadrature weights
+  quadrature_weights = q.get_weights();
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::InternalData::initialize_face(
+  const UpdateFlags      update_flags,
+  const Quadrature<dim> &q,
+  const unsigned int     n_original_q_points)
+{
+  initialize(update_flags, q, n_original_q_points);
+
+  if (this->update_each &
+      (update_boundary_forms | update_normal_vectors | update_jacobians |
+       update_JxW_values | update_inverse_jacobians))
+    {
+      aux.resize(dim - 1,
+                 std::vector<Tensor<1, spacedim>>(n_original_q_points));
+
+      // Compute tangentials to the unit cell.
+      if (this->fe.reference_cell_type() == ReferenceCell::get_hypercube(dim))
+        {
+          for (const unsigned int i : GeometryInfo<dim>::face_indices())
+            {
+              unit_tangentials[i].resize(n_original_q_points);
+              std::fill(unit_tangentials[i].begin(),
+                        unit_tangentials[i].end(),
+                        GeometryInfo<dim>::unit_tangential_vectors[i][0]);
+              if (dim > 2)
+                {
+                  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
+                    .resize(n_original_q_points);
+                  std::fill(
+                    unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
+                      .begin(),
+                    unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
+                      .end(),
+                    GeometryInfo<dim>::unit_tangential_vectors[i][1]);
+                }
+            }
+        }
+      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tri)
+        {
+          Tensor<1, dim> t1;
+          t1[0] = 1;
+          t1[1] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[0].emplace_back(t1);
+          t1[0] = -std::sqrt(0.5);
+          t1[1] = +std::sqrt(0.5);
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[1].emplace_back(t1);
+          t1[0] = 0;
+          t1[1] = -1;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[2].emplace_back(t1);
+        }
+      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tet)
+        {
+          Tensor<1, dim> t1;
+
+          t1[0] = 0;
+          t1[1] = 1;
+          t1[2] = 0; // face 0
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[0].emplace_back(t1);
+
+          t1[0] = 1;
+          t1[1] = 0;
+          t1[2] = 0; // face 0
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[4].emplace_back(t1);
+
+          t1[0] = 1;
+          t1[1] = 0;
+          t1[2] = 0; // face 1
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[1].emplace_back(t1);
+
+          t1[0] = 0;
+          t1[1] = 0;
+          t1[2] = 1; // face 1
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[5].emplace_back(t1);
+
+          t1[0] = 0;
+          t1[1] = 0;
+          t1[2] = 1; // face 2
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[2].emplace_back(t1);
+
+          t1[0] = 0;
+          t1[1] = 1;
+          t1[2] = 0; // face 2
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[6].emplace_back(t1);
+
+          t1[0] = -std::pow(1.0 / 3.0, 1.0 / 4.0);
+          t1[1] = +std::pow(1.0 / 3.0, 1.0 / 4.0);
+          t1[2] = +0; // face 3
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[3].emplace_back(t1);
+
+          t1[0] = -std::pow(1.0 / 3.0, 1.0 / 4.0);
+          t1[1] = +0;
+          t1[2] = +std::pow(1.0 / 3.0, 1.0 / 4.0); // face 3
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[7].emplace_back(t1);
+        }
+      else
+        {
+          Assert(false, ExcNotImplemented());
+        }
+    }
+}
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::InternalData::compute_shape_function_values(
+  const std::vector<Point<dim>> &unit_points)
+{
+  const auto fe_poly = dynamic_cast<const FE_Poly<dim, spacedim> *>(&this->fe);
+
+  Assert(fe_poly != nullptr, ExcNotImplemented());
+
+  const auto &tensor_pols = fe_poly->get_poly_space();
+
+  const unsigned int n_shape_functions = fe.n_dofs_per_cell();
+  const unsigned int n_points          = unit_points.size();
+
+  std::vector<double>         values;
+  std::vector<Tensor<1, dim>> grads;
+  if (shape_values.size() != 0)
+    {
+      Assert(shape_values.size() == n_shape_functions * n_points,
+             ExcInternalError());
+      values.resize(n_shape_functions);
+    }
+  if (shape_derivatives.size() != 0)
+    {
+      Assert(shape_derivatives.size() == n_shape_functions * n_points,
+             ExcInternalError());
+      grads.resize(n_shape_functions);
+    }
+
+  std::vector<Tensor<2, dim>> grad2;
+  if (shape_second_derivatives.size() != 0)
+    {
+      Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
+             ExcInternalError());
+      grad2.resize(n_shape_functions);
+    }
+
+  std::vector<Tensor<3, dim>> grad3;
+  if (shape_third_derivatives.size() != 0)
+    {
+      Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
+             ExcInternalError());
+      grad3.resize(n_shape_functions);
+    }
+
+  std::vector<Tensor<4, dim>> grad4;
+  if (shape_fourth_derivatives.size() != 0)
+    {
+      Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
+             ExcInternalError());
+      grad4.resize(n_shape_functions);
+    }
+
+
+  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
+      shape_second_derivatives.size() != 0 ||
+      shape_third_derivatives.size() != 0 ||
+      shape_fourth_derivatives.size() != 0)
+    for (unsigned int point = 0; point < n_points; ++point)
+      {
+        tensor_pols.evaluate(
+          unit_points[point], values, grads, grad2, grad3, grad4);
+
+        if (shape_values.size() != 0)
+          for (unsigned int i = 0; i < n_shape_functions; ++i)
+            shape(point, i) = values[i];
+
+        if (shape_derivatives.size() != 0)
+          for (unsigned int i = 0; i < n_shape_functions; ++i)
+            derivative(point, i) = grads[i];
+
+        if (shape_second_derivatives.size() != 0)
+          for (unsigned int i = 0; i < n_shape_functions; ++i)
+            second_derivative(point, i) = grad2[i];
+
+        if (shape_third_derivatives.size() != 0)
+          for (unsigned int i = 0; i < n_shape_functions; ++i)
+            third_derivative(point, i) = grad3[i];
+
+        if (shape_fourth_derivatives.size() != 0)
+          for (unsigned int i = 0; i < n_shape_functions; ++i)
+            fourth_derivative(point, i) = grad4[i];
+      }
+}
+
+
+namespace internal
+{
+  namespace MappingFEImplementation
+  {
+    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>
+      void
+      maybe_compute_q_points(
+        const typename QProjector<dim>::DataSetDescriptor              data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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)
+            {
+              const double *  shape = &data.shape(point + data_set, 0);
+              Point<spacedim> result =
+                (shape[0] * data.mapping_support_points[0]);
+              for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  result[i] += shape[k] * data.mapping_support_points[k][i];
+              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>
+      void
+      maybe_update_Jacobians(
+        const CellSimilarity::Similarity cell_similarity,
+        const typename dealii::QProjector<dim>::DataSetDescriptor      data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data)
+      {
+        const UpdateFlags update_flags = data.update_each;
+
+        if (update_flags & update_contravariant_transformation)
+          // if the current cell is just a
+          // translation of the previous one, no
+          // need to recompute jacobians...
+          if (cell_similarity != CellSimilarity::translation)
+            {
+              const unsigned int n_q_points = data.contravariant.size();
+
+              std::fill(data.contravariant.begin(),
+                        data.contravariant.end(),
+                        DerivativeForm<1, dim, spacedim>());
+
+              Assert(data.n_shape_functions > 0, ExcInternalError());
+
+              const Tensor<1, spacedim> *supp_pts =
+                data.mapping_support_points.data();
+
+              for (unsigned int point = 0; point < n_q_points; ++point)
+                {
+                  const Tensor<1, dim> *data_derv =
+                    &data.derivative(point + data_set, 0);
+
+                  double result[spacedim][dim];
+
+                  // peel away part of sum to avoid zeroing the
+                  // entries and adding for the first time
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    for (unsigned int j = 0; j < dim; ++j)
+                      result[i][j] = data_derv[0][j] * supp_pts[0][i];
+                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                    for (unsigned int i = 0; i < spacedim; ++i)
+                      for (unsigned int j = 0; j < dim; ++j)
+                        result[i][j] += data_derv[k][j] * supp_pts[k][i];
+
+                  // write result into contravariant data. for
+                  // j=dim in the case dim<spacedim, there will
+                  // never be any nonzero data that arrives in
+                  // here, so it is ok anyway because it was
+                  // initialized to zero at the initialization
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    for (unsigned int j = 0; j < dim; ++j)
+                      data.contravariant[point][i][j] = result[i][j];
+                }
+            }
+
+        if (update_flags & update_covariant_transformation)
+          if (cell_similarity != CellSimilarity::translation)
+            {
+              const unsigned int n_q_points = data.contravariant.size();
+              for (unsigned int point = 0; point < n_q_points; ++point)
+                {
+                  data.covariant[point] =
+                    (data.contravariant[point]).covariant_form();
+                }
+            }
+
+        if (update_flags & update_volume_elements)
+          if (cell_similarity != CellSimilarity::translation)
+            {
+              const unsigned int n_q_points = data.contravariant.size();
+              for (unsigned int point = 0; point < n_q_points; ++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>
+      void
+      maybe_update_jacobian_grads(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              for (unsigned int point = 0; point < n_q_points; ++point)
+                {
+                  const Tensor<2, dim> *second =
+                    &data.second_derivative(point + data_set, 0);
+                  double result[spacedim][dim][dim];
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        result[i][j][l] =
+                          (second[0][j][l] * data.mapping_support_points[0][i]);
+                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                    for (unsigned int i = 0; i < spacedim; ++i)
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          result[i][j][l] +=
+                            (second[k][j][l] *
+                             data.mapping_support_points[k][i]);
+
+                  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>
+      void
+      maybe_update_jacobian_pushed_forward_grads(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              {
+                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);
+                    double result[spacedim][dim][dim];
+                    for (unsigned int i = 0; i < spacedim; ++i)
+                      for (unsigned int j = 0; j < dim; ++j)
+                        for (unsigned int l = 0; l < dim; ++l)
+                          result[i][j][l] = (second[0][j][l] *
+                                             data.mapping_support_points[0][i]);
+                    for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                      for (unsigned int i = 0; i < spacedim; ++i)
+                        for (unsigned int j = 0; j < dim; ++j)
+                          for (unsigned int l = 0; l < dim; ++l)
+                            result[i][j][l] +=
+                              (second[k][j][l] *
+                               data.mapping_support_points[k][i]);
+
+                    // 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 derivatives 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>
+      void
+      maybe_update_jacobian_2nd_derivatives(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              {
+                for (unsigned int point = 0; point < n_q_points; ++point)
+                  {
+                    const Tensor<3, dim> *third =
+                      &data.third_derivative(point + data_set, 0);
+                    double result[spacedim][dim][dim][dim];
+                    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)
+                            result[i][j][l][m] =
+                              (third[0][j][l][m] *
+                               data.mapping_support_points[0][i]);
+                    for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                      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)
+                              result[i][j][l][m] +=
+                                (third[k][j][l][m] *
+                                 data.mapping_support_points[k][i]);
+
+                    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 Hessian of the Hessian 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>
+      void
+      maybe_update_jacobian_pushed_forward_2nd_derivatives(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              {
+                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);
+                    double result[spacedim][dim][dim][dim];
+                    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)
+                            result[i][j][l][m] =
+                              (third[0][j][l][m] *
+                               data.mapping_support_points[0][i]);
+                    for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                      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)
+                              result[i][j][l][m] +=
+                                (third[k][j][l][m] *
+                                 data.mapping_support_points[k][i]);
+
+                    // 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 derivatives 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>
+      void
+      maybe_update_jacobian_3rd_derivatives(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              {
+                for (unsigned int point = 0; point < n_q_points; ++point)
+                  {
+                    const Tensor<4, dim> *fourth =
+                      &data.fourth_derivative(point + data_set, 0);
+                    double result[spacedim][dim][dim][dim][dim];
+                    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)
+                              result[i][j][l][m][n] =
+                                (fourth[0][j][l][m][n] *
+                                 data.mapping_support_points[0][i]);
+                    for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                      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)
+                                result[i][j][l][m][n] +=
+                                  (fourth[k][j][l][m][n] *
+                                   data.mapping_support_points[k][i]);
+
+                    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 Hessian gradient 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>
+      void
+      maybe_update_jacobian_pushed_forward_3rd_derivatives(
+        const CellSimilarity::Similarity                  cell_similarity,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        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();
+
+            if (cell_similarity != CellSimilarity::translation)
+              {
+                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);
+                    double result[spacedim][dim][dim][dim][dim];
+                    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)
+                              result[i][j][l][m][n] =
+                                (fourth[0][j][l][m][n] *
+                                 data.mapping_support_points[0][i]);
+                    for (unsigned int k = 1; k < data.n_shape_functions; ++k)
+                      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)
+                                result[i][j][l][m][n] +=
+                                  (fourth[k][j][l][m][n] *
+                                   data.mapping_support_points[k][i]);
+
+                    // 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];
+                              }
+                  }
+              }
+          }
+      }
+    } // namespace
+  }   // namespace MappingFEImplementation
+} // namespace internal
+
+
+
+template <int dim, int spacedim>
+MappingFE<dim, spacedim>::MappingFE(const FiniteElement<dim, spacedim> &fe)
+  : fe(fe.clone())
+  , polynomial_degree(fe.tensor_degree())
+{
+  Assert(polynomial_degree >= 1,
+         ExcMessage("It only makes sense to create polynomial mappings "
+                    "with a polynomial degree greater or equal to one."));
+}
+
+
+
+template <int dim, int spacedim>
+MappingFE<dim, spacedim>::MappingFE(const MappingFE<dim, spacedim> &mapping)
+  : fe(mapping.fe->clone())
+  , polynomial_degree(mapping.polynomial_degree)
+{}
+
+
+
+template <int dim, int spacedim>
+std::unique_ptr<Mapping<dim, spacedim>>
+MappingFE<dim, spacedim>::clone() const
+{
+  return std::make_unique<MappingFE<dim, spacedim>>(*this);
+}
+
+
+
+template <int dim, int spacedim>
+unsigned int
+MappingFE<dim, spacedim>::get_degree() const
+{
+  return polynomial_degree;
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+MappingFE<dim, spacedim>::transform_unit_to_real_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const Point<dim> &                                          p) const
+{
+  const std::vector<Point<spacedim>> support_points =
+    this->compute_mapping_support_points(cell);
+
+  Point<spacedim> mapped_point;
+
+  for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
+    mapped_point += support_points[i] * this->fe->shape_value(i, p);
+
+  return mapped_point;
+}
+
+
+
+template <int dim, int spacedim>
+Point<dim>
+MappingFE<dim, spacedim>::transform_real_to_unit_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const Point<spacedim> &                                     p) const
+{
+  Assert(false, StandardExceptions::ExcNotImplemented());
+
+  (void)cell;
+  (void)p;
+
+  return Point<dim>();
+}
+
+
+
+template <int dim, int spacedim>
+UpdateFlags
+MappingFE<dim, spacedim>::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 the 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 5 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_JxW_values |
+                 update_jacobians | update_jacobian_grads |
+                 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 JxW values
+      // to the list of flags to be updated for each cell.
+      if (out & update_contravariant_transformation)
+        out |= update_volume_elements;
+
+      // the same is true when computing normal vectors: they require
+      // the determinant of the Jacobian
+      if (out & update_normal_vectors)
+        out |= update_volume_elements;
+    }
+
+  return out;
+}
+
+
+
+template <int dim, int spacedim>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFE<dim, spacedim>::get_data(const UpdateFlags      update_flags,
+                                   const Quadrature<dim> &q) const
+{
+  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
+    std::make_unique<InternalData>(*this->fe);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+  data.initialize(this->requires_update_flags(update_flags), q, q.size());
+
+  return data_ptr;
+}
+
+
+
+template <int dim, int spacedim>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFE<dim, spacedim>::get_face_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>(*this->fe);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+  data.initialize_face(this->requires_update_flags(update_flags),
+                       QProjector<dim>::project_to_all_faces(
+                         this->fe->reference_cell_type(), quadrature),
+                       quadrature.size());
+
+  return data_ptr;
+}
+
+
+
+template <int dim, int spacedim>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+MappingFE<dim, spacedim>::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>(*this->fe);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+  data.initialize_face(this->requires_update_flags(update_flags),
+                       QProjector<dim>::project_to_all_subfaces(
+                         this->fe->reference_cell_type(), quadrature),
+                       quadrature.size());
+
+  return data_ptr;
+}
+
+
+
+template <int dim, int spacedim>
+CellSimilarity::Similarity
+MappingFE<dim, spacedim>::fill_fe_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const CellSimilarity::Similarity                            cell_similarity,
+  const Quadrature<dim> &                                     quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  // ensure that the following static_cast is really correct:
+  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();
+
+  // recompute the support points of the transformation of this
+  // cell. we tried to be clever here in an earlier version of the
+  // library by checking whether the cell is the same as the one we
+  // had visited last, but it turns out to be difficult to determine
+  // that because a cell for the purposes of a mapping is
+  // characterized not just by its (triangulation, level, index)
+  // triple, but also by the locations of its vertices, the manifold
+  // object attached to the cell and all of its bounding faces/edges,
+  // etc. to reliably test that the "cell" we are on is, therefore,
+  // not easily done
+  data.mapping_support_points = this->compute_mapping_support_points(cell);
+  data.cell_of_current_support_points = cell;
+
+  // if the order of the mapping is greater than 1, then do not reuse any cell
+  // similarity information. This is necessary because the cell similarity
+  // value is computed with just cell vertices and does not take into account
+  // cell curvature.
+  const CellSimilarity::Similarity computed_cell_similarity =
+    (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
+
+  internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
+    QProjector<dim>::DataSetDescriptor::cell(),
+    data,
+    output_data.quadrature_points);
+
+  internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
+    computed_cell_similarity, QProjector<dim>::DataSetDescriptor::cell(), data);
+
+  internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
+    computed_cell_similarity,
+    QProjector<dim>::DataSetDescriptor::cell(),
+    data,
+    output_data.jacobian_grads);
+
+  internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads<
+    dim,
+    spacedim>(computed_cell_similarity,
+              QProjector<dim>::DataSetDescriptor::cell(),
+              data,
+              output_data.jacobian_pushed_forward_grads);
+
+  internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives<
+    dim,
+    spacedim>(computed_cell_similarity,
+              QProjector<dim>::DataSetDescriptor::cell(),
+              data,
+              output_data.jacobian_2nd_derivatives);
+
+  internal::MappingFEImplementation::
+    maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
+      computed_cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      output_data.jacobian_pushed_forward_2nd_derivatives);
+
+  internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives<
+    dim,
+    spacedim>(computed_cell_similarity,
+              QProjector<dim>::DataSetDescriptor::cell(),
+              data,
+              output_data.jacobian_3rd_derivatives);
+
+  internal::MappingFEImplementation::
+    maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
+      computed_cell_similarity,
+      QProjector<dim>::DataSetDescriptor::cell(),
+      data,
+      output_data.jacobian_pushed_forward_3rd_derivatives);
+
+  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));
+
+
+      if (computed_cell_similarity != CellSimilarity::translation)
+        for (unsigned int point = 0; point < n_q_points; ++point)
+          {
+            if (dim == spacedim)
+              {
+                const double det = data.contravariant[point].determinant();
+
+                // 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 (computed_cell_similarity ==
+                    CellSimilarity::inverted_translation)
+                  {
+                    // we only need to flip the normal
+                    if (update_flags & update_normal_vectors)
+                      output_data.normal_vectors[point] *= -1.;
+                  }
+                else
+                  {
+                    if (update_flags & update_normal_vectors)
+                      {
+                        Assert(spacedim == dim + 1,
+                               ExcMessage(
+                                 "There is no (unique) cell normal for " +
+                                 Utilities::int_to_string(dim) +
+                                 "-dimensional cells in " +
+                                 Utilities::int_to_string(spacedim) +
+                                 "-dimensional space. This only works if the "
+                                 "space dimension is one greater than the "
+                                 "dimensionality of the mesh cells."));
+
+                        if (dim == 1)
+                          output_data.normal_vectors[point] =
+                            cross_product_2d(-DX_t[0]);
+                        else // dim == 2
+                          output_data.normal_vectors[point] =
+                            cross_product_3d(DX_t[0], DX_t[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);
+      if (computed_cell_similarity != CellSimilarity::translation)
+        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);
+      if (computed_cell_similarity != CellSimilarity::translation)
+        for (unsigned int point = 0; point < n_q_points; ++point)
+          output_data.inverse_jacobians[point] =
+            data.covariant[point].transpose();
+    }
+
+  return computed_cell_similarity;
+}
+
+
+
+namespace internal
+{
+  namespace MappingFEImplementation
+  {
+    namespace
+    {
+      /**
+       * 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>
+      void
+      maybe_compute_face_data(
+        const dealii::MappingFE<dim, spacedim> &mapping,
+        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+          &                                               cell,
+        const unsigned int                                face_no,
+        const unsigned int                                subface_no,
+        const unsigned int                                n_q_points,
+        const typename QProjector<dim>::DataSetDescriptor data_set,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+          &output_data)
+      {
+        const UpdateFlags update_flags = data.update_each;
+
+        if (update_flags &
+            (update_boundary_forms | update_normal_vectors | update_jacobians |
+             update_JxW_values | update_inverse_jacobians))
+          {
+            if (update_flags & update_boundary_forms)
+              AssertDimension(output_data.boundary_forms.size(), n_q_points);
+            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);
+
+            Assert(data.aux.size() + 1 >= dim, ExcInternalError());
+
+            // first compute some common data that is used for evaluating
+            // all of the flags below
+
+            // 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 (update_flags & update_boundary_forms)
+              {
+                // 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:
+                            Assert(false, ExcNotImplemented());
+                        }
+                  }
+                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_JxW_values)
+              for (unsigned int i = 0; i < output_data.boundary_forms.size();
+                   ++i)
+                {
+                  output_data.JxW_values[i] =
+                    output_data.boundary_forms[i].norm() *
+                    data.quadrature_weights[i + data_set];
+
+                  if (subface_no != numbers::invalid_unsigned_int)
+                    {
+#if false
+                       const double area_ratio =
+                        GeometryInfo<dim>::subface_ratio(
+                          cell->subface_case(face_no), subface_no);
+                       output_data.JxW_values[i] *= area_ratio;
+#else
+                      Assert(false, ExcNotImplemented());
+#endif
+                    }
+                }
+
+            if (update_flags & update_normal_vectors)
+              for (unsigned int i = 0; i < output_data.normal_vectors.size();
+                   ++i)
+                output_data.normal_vectors[i] =
+                  Point<spacedim>(output_data.boundary_forms[i] /
+                                  output_data.boundary_forms[i].norm());
+
+            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();
+          }
+      }
+
+
+      /**
+       * Do the work of MappingFE::fill_fe_face_values() and
+       * MappingFE::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>
+      void
+      do_fill_fe_face_values(
+        const dealii::MappingFE<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 Quadrature<dim - 1> &                       quadrature,
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+          &output_data)
+      {
+        maybe_compute_q_points<dim, spacedim>(data_set,
+                                              data,
+                                              output_data.quadrature_points);
+        maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
+                                              data_set,
+                                              data);
+        maybe_update_jacobian_grads<dim, spacedim>(CellSimilarity::none,
+                                                   data_set,
+                                                   data,
+                                                   output_data.jacobian_grads);
+        maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
+          CellSimilarity::none,
+          data_set,
+          data,
+          output_data.jacobian_pushed_forward_grads);
+        maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
+          CellSimilarity::none,
+          data_set,
+          data,
+          output_data.jacobian_2nd_derivatives);
+        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
+          CellSimilarity::none,
+          data_set,
+          data,
+          output_data.jacobian_pushed_forward_2nd_derivatives);
+        maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
+          CellSimilarity::none,
+          data_set,
+          data,
+          output_data.jacobian_3rd_derivatives);
+        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
+          CellSimilarity::none,
+          data_set,
+          data,
+          output_data.jacobian_pushed_forward_3rd_derivatives);
+
+        maybe_compute_face_data(mapping,
+                                cell,
+                                face_no,
+                                subface_no,
+                                quadrature.size(),
+                                data_set,
+                                data,
+                                output_data);
+      }
+    } // namespace
+  }   // namespace MappingFEImplementation
+} // namespace internal
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::fill_fe_face_values(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const Quadrature<dim - 1> &                                 quadrature,
+  const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
+  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &output_data) const
+{
+  // ensure that the following cast is really correct:
+  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  // if necessary, recompute the support points of the transformation of this
+  // cell (note that we need to first check the triangulation pointer, since
+  // otherwise the second test might trigger an exception if the triangulations
+  // are not the same)
+  if ((data.mapping_support_points.size() == 0) ||
+      (&cell->get_triangulation() !=
+       &data.cell_of_current_support_points->get_triangulation()) ||
+      (cell != data.cell_of_current_support_points))
+    {
+      data.mapping_support_points = this->compute_mapping_support_points(cell);
+      data.cell_of_current_support_points = cell;
+    }
+
+  internal::MappingFEImplementation::do_fill_fe_face_values(
+    *this,
+    cell,
+    face_no,
+    numbers::invalid_unsigned_int,
+    QProjector<dim>::DataSetDescriptor::face(this->fe->reference_cell_type(),
+                                             face_no,
+                                             cell->face_orientation(face_no),
+                                             cell->face_flip(face_no),
+                                             cell->face_rotation(face_no),
+                                             quadrature.size()),
+    quadrature,
+    data,
+    output_data);
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::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
+{
+  // ensure that the following cast is really correct:
+  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
+         ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+  // if necessary, recompute the support points of the transformation of this
+  // cell (note that we need to first check the triangulation pointer, since
+  // otherwise the second test might trigger an exception if the triangulations
+  // are not the same)
+  if ((data.mapping_support_points.size() == 0) ||
+      (&cell->get_triangulation() !=
+       &data.cell_of_current_support_points->get_triangulation()) ||
+      (cell != data.cell_of_current_support_points))
+    {
+      data.mapping_support_points = this->compute_mapping_support_points(cell);
+      data.cell_of_current_support_points = cell;
+    }
+
+  internal::MappingFEImplementation::do_fill_fe_face_values(
+    *this,
+    cell,
+    face_no,
+    subface_no,
+    QProjector<dim>::DataSetDescriptor::subface(this->fe->reference_cell_type(),
+                                                face_no,
+                                                subface_no,
+                                                cell->face_orientation(face_no),
+                                                cell->face_flip(face_no),
+                                                cell->face_rotation(face_no),
+                                                quadrature.size(),
+                                                cell->subface_case(face_no)),
+    quadrature,
+    data,
+    output_data);
+}
+
+
+
+namespace internal
+{
+  namespace MappingFEImplementation
+  {
+    namespace
+    {
+      template <int dim, int spacedim, int rank>
+      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::MappingFE<dim, spacedim>::InternalData *>(
+             &mapping_data) != nullptr),
+          ExcInternalError());
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data =
+          static_cast<
+            const typename dealii::MappingFE<dim, spacedim>::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"));
+                if (rank != 1)
+                  return;
+
+                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_covariant_transformation"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  output[i] = apply_transformation(data.covariant[i], input[i]);
+
+                return;
+              }
+
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+      }
+
+
+      template <int dim, int spacedim, int rank>
+      void
+      transform_gradients(
+        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::MappingFE<dim, spacedim>::InternalData *>(
+             &mapping_data) != nullptr),
+          ExcInternalError());
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data =
+          static_cast<
+            const typename dealii::MappingFE<dim, spacedim>::InternalData &>(
+            mapping_data);
+
+        switch (mapping_kind)
+          {
+            case mapping_contravariant_gradient:
+              {
+                Assert(
+                  data.update_each & update_covariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+                Assert(rank == 2, ExcMessage("Only for rank 2"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  {
+                    const DerivativeForm<1, spacedim, dim> A =
+                      apply_transformation(data.contravariant[i],
+                                           transpose(input[i]));
+                    output[i] =
+                      apply_transformation(data.covariant[i], A.transpose());
+                  }
+
+                return;
+              }
+
+            case mapping_covariant_gradient:
+              {
+                Assert(
+                  data.update_each & update_covariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+                Assert(rank == 2, ExcMessage("Only for rank 2"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  {
+                    const DerivativeForm<1, spacedim, dim> A =
+                      apply_transformation(data.covariant[i],
+                                           transpose(input[i]));
+                    output[i] =
+                      apply_transformation(data.covariant[i], A.transpose());
+                  }
+
+                return;
+              }
+
+            case mapping_piola_gradient:
+              {
+                Assert(
+                  data.update_each & update_covariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+                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 == 2, ExcMessage("Only for rank 2"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  {
+                    const DerivativeForm<1, spacedim, dim> A =
+                      apply_transformation(data.covariant[i], input[i]);
+                    const Tensor<2, spacedim> T =
+                      apply_transformation(data.contravariant[i],
+                                           A.transpose());
+
+                    output[i] = transpose(T);
+                    output[i] /= data.volume_elements[i];
+                  }
+
+                return;
+              }
+
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+      }
+
+
+
+      template <int dim, int spacedim>
+      void
+      transform_hessians(
+        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)
+      {
+        AssertDimension(input.size(), output.size());
+        Assert(
+          (dynamic_cast<
+             const typename dealii::MappingFE<dim, spacedim>::InternalData *>(
+             &mapping_data) != nullptr),
+          ExcInternalError());
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data =
+          static_cast<
+            const typename dealii::MappingFE<dim, spacedim>::InternalData &>(
+            mapping_data);
+
+        switch (mapping_kind)
+          {
+            case mapping_contravariant_hessian:
+              {
+                Assert(
+                  data.update_each & update_covariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_contravariant_transformation"));
+
+                for (unsigned int q = 0; q < output.size(); ++q)
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    {
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
+                        for (unsigned int K = 0; K < dim; ++K)
+                          {
+                            tmp1[J][K] =
+                              data.contravariant[q][i][0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] +=
+                                data.contravariant[q][i][I] * input[q][I][J][K];
+                          }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] =
+                                data.covariant[q][k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] +=
+                                  data.covariant[q][k][K] * tmp2[K];
+                            }
+                        }
+                    }
+                return;
+              }
+
+            case mapping_covariant_hessian:
+              {
+                Assert(
+                  data.update_each & update_covariant_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)
+                    {
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
+                        for (unsigned int K = 0; K < dim; ++K)
+                          {
+                            tmp1[J][K] =
+                              data.covariant[q][i][0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] +=
+                                data.covariant[q][i][I] * input[q][I][J][K];
+                          }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] =
+                                data.covariant[q][k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] +=
+                                  data.covariant[q][k][K] * tmp2[K];
+                            }
+                        }
+                    }
+
+                return;
+              }
+
+            case mapping_piola_hessian:
+              {
+                Assert(
+                  data.update_each & update_covariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+                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"));
+
+                for (unsigned int q = 0; q < output.size(); ++q)
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    {
+                      double factor[dim];
+                      for (unsigned int I = 0; I < dim; ++I)
+                        factor[I] =
+                          data.contravariant[q][i][I] / data.volume_elements[q];
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
+                        for (unsigned int K = 0; K < dim; ++K)
+                          {
+                            tmp1[J][K] = factor[0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] += factor[I] * input[q][I][J][K];
+                          }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] =
+                                data.covariant[q][k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] +=
+                                  data.covariant[q][k][K] * tmp2[K];
+                            }
+                        }
+                    }
+
+                return;
+              }
+
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+      }
+
+
+
+      template <int dim, int spacedim, int rank>
+      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::MappingFE<dim, spacedim>::InternalData *>(
+             &mapping_data) != nullptr),
+          ExcInternalError());
+        const typename dealii::MappingFE<dim, spacedim>::InternalData &data =
+          static_cast<
+            const typename dealii::MappingFE<dim, spacedim>::InternalData &>(
+            mapping_data);
+
+        switch (mapping_kind)
+          {
+            case mapping_covariant:
+              {
+                Assert(
+                  data.update_each & update_contravariant_transformation,
+                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                    "update_covariant_transformation"));
+
+                for (unsigned int i = 0; i < output.size(); ++i)
+                  output[i] = apply_transformation(data.covariant[i], input[i]);
+
+                return;
+              }
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+      }
+    } // namespace
+  }   // namespace MappingFEImplementation
+} // namespace internal
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::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
+{
+  internal::MappingFEImplementation::transform_fields(input,
+                                                      mapping_kind,
+                                                      mapping_data,
+                                                      output);
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::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
+{
+  internal::MappingFEImplementation::transform_differential_forms(input,
+                                                                  mapping_kind,
+                                                                  mapping_data,
+                                                                  output);
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::transform(
+  const ArrayView<const Tensor<2, dim>> &                  input,
+  const MappingKind                                        mapping_kind,
+  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+  const ArrayView<Tensor<2, spacedim>> &                   output) const
+{
+  switch (mapping_kind)
+    {
+      case mapping_contravariant:
+        internal::MappingFEImplementation::transform_fields(input,
+                                                            mapping_kind,
+                                                            mapping_data,
+                                                            output);
+        return;
+
+      case mapping_piola_gradient:
+      case mapping_contravariant_gradient:
+      case mapping_covariant_gradient:
+        internal::MappingFEImplementation::transform_gradients(input,
+                                                               mapping_kind,
+                                                               mapping_data,
+                                                               output);
+        return;
+      default:
+        Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::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)
+                {
+                  double tmp[dim];
+                  for (unsigned int K = 0; K < dim; ++K)
+                    {
+                      tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
+                      for (unsigned int J = 1; J < dim; ++J)
+                        tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
+                    }
+                  for (unsigned int k = 0; k < spacedim; ++k)
+                    {
+                      output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
+                      for (unsigned int K = 1; K < dim; ++K)
+                        output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
+                    }
+                }
+          return;
+        }
+
+      default:
+        Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingFE<dim, spacedim>::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
+{
+  switch (mapping_kind)
+    {
+      case mapping_piola_hessian:
+      case mapping_contravariant_hessian:
+      case mapping_covariant_hessian:
+        internal::MappingFEImplementation::transform_hessians(input,
+                                                              mapping_kind,
+                                                              mapping_data,
+                                                              output);
+        return;
+      default:
+        Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<Point<spacedim>>
+MappingFE<dim, spacedim>::compute_mapping_support_points(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
+{
+  // get the vertices first
+  std::vector<Point<spacedim>> a(cell->n_vertices());
+
+  for (const unsigned int i : cell->vertex_indices())
+    a[i] = cell->vertex(i);
+
+  if (this->polynomial_degree > 1)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+  return a;
+}
+
+
+
+//--------------------------- Explicit instantiations -----------------------
+#include "mapping_fe.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/mapping_fe.inst.in b/source/fe/mapping_fe.inst.in
new file mode 100644 (file)
index 0000000..a29e0d5
--- /dev/null
@@ -0,0 +1,22 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+    template class MappingFE<deal_II_dimension, deal_II_space_dimension>;
+#endif
+  }

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.