]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Reorganize the internal mapping functions. 4592/head
authorDavid Wells <wellsd2@rpi.edu>
Sat, 8 Jul 2017 18:48:48 +0000 (14:48 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Sun, 9 Jul 2017 19:10:57 +0000 (15:10 -0400)
If we want to ultimately use a unity build then we will need to avoid duplicate
function names. This commit adds a new policy for avoiding these (a namespacing
scheme) and applies the fix to the mapping classes.

doc/doxygen/headers/coding_conventions.h
source/fe/mapping_fe_field.cc
source/fe/mapping_manifold.cc
source/fe/mapping_q_generic.cc

index 0ae51f4d933c464767bcd0f43436b3b1b95f3f18..3b2ce05ffe35aa9a43c181016737ba2ca6e4d6ed 100644 (file)
@@ -130,6 +130,26 @@ source files.
   handle run time errors (like I/O failures) which must be on
   in any case, not only in debug mode.</li>
 
+<li> Sometimes it makes sense to implement a class by using several
+  non-member functions that are not part of the public interface and are only
+  meant to be called in the current source file. Such free functions should be
+  put in an anonymous namespace structured in the following way:
+  <code>
+  <pre>
+  namespace internal
+  {
+    namespace ClassName
+    {
+      namespace
+      {
+        // free functions go here
+      }
+    }
+  }
+  </pre>
+  </code>
+  where <code>ClassName</code> is the name of the calling class.
+
 <li> Classes and types generally are named using uppercase letters to denote
   word beginnings (e.g. TriaIterator) &mdash; sometimes called
   <a href="http://en.wikipedia.org/wiki/Camel_case"><i>camel
index a91a96b8613a28a95475719568d0fde4a3c79393..e4cad13d82492bebd4e74f1091d1ad08cde73e80 100644 (file)
@@ -520,774 +520,777 @@ MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::get_subface_data
 
 namespace internal
 {
-  namespace
+  namespace MappingFEField
   {
-    /**
-     * Compute the locations of quadrature points on the object described by
-     * the first argument (and the cell for which the mapping support points
-     * have already been set), but only if the update_flags of the @p data
-     * argument indicate so.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_compute_q_points (const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-                            const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                            const FiniteElement<dim, spacedim>   &fe,
-                            const ComponentMask                  &fe_mask,
-                            const std::vector<unsigned int>      &fe_to_real,
-                            std::vector<Point<spacedim> >        &quadrature_points)
+    namespace
     {
-      const UpdateFlags update_flags = data.update_each;
+      /**
+       * Compute the locations of quadrature points on the object described by
+       * the first argument (and the cell for which the mapping support points
+       * have already been set), but only if the update_flags of the @p data
+       * argument indicate so.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_compute_q_points (const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+                              const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                              const FiniteElement<dim, spacedim>   &fe,
+                              const ComponentMask                  &fe_mask,
+                              const std::vector<unsigned int>      &fe_to_real,
+                              std::vector<Point<spacedim> >        &quadrature_points)
+      {
+        const UpdateFlags update_flags = data.update_each;
 
-      if (update_flags & update_quadrature_points)
-        {
-          for (unsigned int point=0; point<quadrature_points.size(); ++point)
-            {
-              Point<spacedim> result;
-              const double *shape = &data.shape(point+data_set,0);
+        if (update_flags & update_quadrature_points)
+          {
+            for (unsigned int point=0; point<quadrature_points.size(); ++point)
+              {
+                Point<spacedim> result;
+                const double *shape = &data.shape(point+data_set,0);
 
-              for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                {
-                  unsigned int comp_k = fe.system_to_component_index(k).first;
-                  if (fe_mask[comp_k])
-                    result[fe_to_real[comp_k]] += data.local_dof_values[k] * shape[k];
-                }
+                for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                  {
+                    unsigned int comp_k = fe.system_to_component_index(k).first;
+                    if (fe_mask[comp_k])
+                      result[fe_to_real[comp_k]] += data.local_dof_values[k] * shape[k];
+                  }
 
-              quadrature_points[point] = result;
-            }
-        }
-    }
+                quadrature_points[point] = result;
+              }
+          }
+      }
 
-    /**
-     * Update the co- and contravariant matrices as well as their determinant,
-     * for the cell described stored in the data object, but only if the
-     * update_flags of the @p data argument indicate so.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_update_Jacobians (const CellSimilarity::Similarity    cell_similarity,
-                            const typename dealii::QProjector<dim>::DataSetDescriptor  data_set,
-                            const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                            const FiniteElement<dim, spacedim> &fe,
-                            const ComponentMask                &fe_mask,
-                            const std::vector<unsigned int>    &fe_to_real)
-    {
-      const UpdateFlags update_flags = data.update_each;
+      /**
+       * Update the co- and contravariant matrices as well as their determinant,
+       * for the cell described stored in the data object, but only if the
+       * update_flags of the @p data argument indicate so.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_Jacobians (const CellSimilarity::Similarity    cell_similarity,
+                              const typename dealii::QProjector<dim>::DataSetDescriptor  data_set,
+                              const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                              const FiniteElement<dim, spacedim> &fe,
+                              const ComponentMask                &fe_mask,
+                              const std::vector<unsigned int>    &fe_to_real)
+      {
+        const UpdateFlags update_flags = data.update_each;
 
-      // then Jacobians
-      if (update_flags & update_contravariant_transformation)
-        {
+        // then Jacobians
+        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();
+            // 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();
 
-              Assert (data.n_shape_functions > 0, ExcInternalError());
+                Assert (data.n_shape_functions > 0, ExcInternalError());
 
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
-                  const Tensor<1,dim> *data_derv =
-                    &data.derivative(point+data_set, 0);
+                for (unsigned int point=0; point<n_q_points; ++point)
+                  {
+                    const Tensor<1,dim> *data_derv =
+                      &data.derivative(point+data_set, 0);
 
-                  Tensor<1, dim> result[spacedim];
+                    Tensor<1, dim> result[spacedim];
 
-                  for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                    {
-                      unsigned int comp_k = fe.system_to_component_index(k).first;
-                      if (fe_mask[comp_k])
-                        result[fe_to_real[comp_k]] += data.local_dof_values[k] * data_derv[k];
-                    }
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          result[fe_to_real[comp_k]] += data.local_dof_values[k] * data_derv[k];
+                      }
 
-                  // write result into contravariant data
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    {
-                      data.contravariant[point][i] = result[i];
-                    }
-                }
-            }
-        }
+                    // write result into contravariant data
+                    for (unsigned int i=0; i<spacedim; ++i)
+                      {
+                        data.contravariant[point][i] = result[i];
+                      }
+                  }
+              }
+          }
 
-      if (update_flags & update_covariant_transformation)
-        {
-          AssertDimension(data.covariant.size(), data.contravariant.size());
-          if (cell_similarity != CellSimilarity::translation)
-            for (unsigned int point=0; point<data.contravariant.size(); ++point)
-              data.covariant[point] = (data.contravariant[point]).covariant_form();
-        }
+        if (update_flags & update_covariant_transformation)
+          {
+            AssertDimension(data.covariant.size(), data.contravariant.size());
+            if (cell_similarity != CellSimilarity::translation)
+              for (unsigned int point=0; point<data.contravariant.size(); ++point)
+                data.covariant[point] = (data.contravariant[point]).covariant_form();
+          }
 
-      if (update_flags & update_volume_elements)
-        {
-          AssertDimension(data.covariant.size(), data.volume_elements.size());
-          if (cell_similarity != CellSimilarity::translation)
-            for (unsigned int point=0; point<data.contravariant.size(); ++point)
-              data.volume_elements[point] = data.contravariant[point].determinant();
-        }
-    }
+        if (update_flags & update_volume_elements)
+          {
+            AssertDimension(data.covariant.size(), data.volume_elements.size());
+            if (cell_similarity != CellSimilarity::translation)
+              for (unsigned int point=0; point<data.contravariant.size(); ++point)
+                data.volume_elements[point] = data.contravariant[point].determinant();
+          }
+      }
 
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_update_jacobian_grads (const CellSimilarity::Similarity              cell_similarity,
-                                 const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-                                 const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                                 const FiniteElement<dim, spacedim>           &fe,
-                                 const ComponentMask                          &fe_mask,
-                                 const std::vector<unsigned int>              &fe_to_real,
-                                 std::vector<DerivativeForm<2,dim,spacedim> > &jacobian_grads)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_grads)
-        {
-          const unsigned int n_q_points = jacobian_grads.size();
+      /**
+       * Update the Hessian of the transformation from unit to real cell, the
+       * Jacobian gradients.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_grads (const CellSimilarity::Similarity              cell_similarity,
+                                   const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+                                   const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                                   const FiniteElement<dim, spacedim>           &fe,
+                                   const ComponentMask                          &fe_mask,
+                                   const std::vector<unsigned int>              &fe_to_real,
+                                   std::vector<DerivativeForm<2,dim,spacedim> > &jacobian_grads)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_grads)
+          {
+            const unsigned int n_q_points = jacobian_grads.size();
 
-          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);
+            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);
 
-                  DerivativeForm<2,dim,spacedim> result;
+                    DerivativeForm<2,dim,spacedim> result;
 
-                  for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                    {
-                      unsigned int comp_k = fe.system_to_component_index(k).first;
-                      if (fe_mask[comp_k])
-                        for (unsigned int j=0; j<dim; ++j)
-                          for (unsigned int l=0; l<dim; ++l)
-                            result[fe_to_real[comp_k]][j][l] += (second[k][j][l]
-                                                                 * data.local_dof_values[k]);
-                    }
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              result[fe_to_real[comp_k]][j][l] += (second[k][j][l]
+                                                                   * data.local_dof_values[k]);
+                      }
 
-                  // never touch any data for j=dim in case dim<spacedim, so
-                  // it will always be zero as it was initialized
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        jacobian_grads[point][i][j][l] = result[i][j][l];
-                }
-            }
-        }
-    }
+                    // never touch any data for j=dim in case dim<spacedim, so
+                    // it will always be zero as it was initialized
+                    for (unsigned int i=0; i<spacedim; ++i)
+                      for (unsigned int j=0; j<dim; ++j)
+                        for (unsigned int l=0; l<dim; ++l)
+                          jacobian_grads[point][i][j][l] = result[i][j][l];
+                  }
+              }
+          }
+      }
 
-    /**
-     * Update the Hessian of the transformation from unit to real cell, the
-     * Jacobian gradients, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_update_jacobian_pushed_forward_grads (
-      const CellSimilarity::Similarity              cell_similarity,
-      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-      const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-      const FiniteElement<dim, spacedim>           &fe,
-      const ComponentMask                          &fe_mask,
-      const std::vector<unsigned int>              &fe_to_real,
-      std::vector<Tensor<3,spacedim> >             &jacobian_pushed_forward_grads )
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_grads)
-        {
-          const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
+      /**
+       * Update the Hessian of the transformation from unit to real cell, the
+       * Jacobian gradients, pushed forward to the real cell coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_pushed_forward_grads
+      (const CellSimilarity::Similarity              cell_similarity,
+       const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+       const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+       const FiniteElement<dim, spacedim>           &fe,
+       const ComponentMask                          &fe_mask,
+       const std::vector<unsigned int>              &fe_to_real,
+       std::vector<Tensor<3,spacedim> >             &jacobian_pushed_forward_grads )
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_grads)
+          {
+            const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
 
-          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);
+            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);
 
-                  DerivativeForm<2,dim,spacedim> result;
+                    DerivativeForm<2,dim,spacedim> result;
 
-                  for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                    {
-                      unsigned int comp_k = fe.system_to_component_index(k).first;
-                      if (fe_mask[comp_k])
-                        for (unsigned int j=0; j<dim; ++j)
-                          for (unsigned int l=0; l<dim; ++l)
-                            result[fe_to_real[comp_k]][j][l] += (second[k][j][l]
-                                                                 * data.local_dof_values[k]);
-                    }
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              result[fe_to_real[comp_k]][j][l] += (second[k][j][l]
+                                                                   * data.local_dof_values[k]);
+                      }
 
-                  // first push forward the j-components
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        {
-                          tmp[i][j][l] = result[i][0][l] *
-                                         data.covariant[point][j][0];
-                          for (unsigned int jr=1; jr<dim; ++jr)
-                            {
-                              tmp[i][j][l] += result[i][jr][l] *
-                                              data.covariant[point][j][jr];
-                            }
-                        }
+                    // 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];
-                            }
+                    // now, pushing forward the l-components
+                    for (unsigned int i=0; i<spacedim; ++i)
+                      for (unsigned int j=0; j<spacedim; ++j)
+                        for (unsigned int l=0; l<spacedim; ++l)
+                          {
+                            jacobian_pushed_forward_grads[point][i][j][l] = tmp[i][j][0] *
+                                                                            data.covariant[point][l][0];
+                            for (unsigned int lr=1; lr<dim; ++lr)
+                              {
+                                jacobian_pushed_forward_grads[point][i][j][l] += tmp[i][j][lr] *
+                                                                                 data.covariant[point][l][lr];
+                              }
 
-                        }
-                }
-            }
-        }
-    }
+                          }
+                  }
+              }
+          }
+      }
 
-    /**
-     * Update the third derivative of the transformation from unit to real
-     * cell, the Jacobian hessians.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity              cell_similarity,
-                                           const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-                                           const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                                           const FiniteElement<dim, spacedim>           &fe,
-                                           const ComponentMask                          &fe_mask,
-                                           const std::vector<unsigned int>              &fe_to_real,
-                                           std::vector<DerivativeForm<3,dim,spacedim> > &jacobian_2nd_derivatives)
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_2nd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_2nd_derivatives.size();
+      /**
+       * Update the third derivative of the transformation from unit to real
+       * cell, the Jacobian hessians.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity              cell_similarity,
+                                             const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+                                             const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                                             const FiniteElement<dim, spacedim>           &fe,
+                                             const ComponentMask                          &fe_mask,
+                                             const std::vector<unsigned int>              &fe_to_real,
+                                             std::vector<DerivativeForm<3,dim,spacedim> > &jacobian_2nd_derivatives)
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_2nd_derivatives)
+          {
+            const unsigned int n_q_points = jacobian_2nd_derivatives.size();
 
-          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);
+            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);
 
-                  DerivativeForm<3,dim,spacedim> result;
+                    DerivativeForm<3,dim,spacedim> result;
 
-                  for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                    {
-                      unsigned int comp_k = fe.system_to_component_index(k).first;
-                      if (fe_mask[comp_k])
-                        for (unsigned int j=0; j<dim; ++j)
-                          for (unsigned int l=0; l<dim; ++l)
-                            for (unsigned int m=0; m<dim; ++m)
-                              result[fe_to_real[comp_k]][j][l][m] += (third[k][j][l][m]
-                                                                      * data.local_dof_values[k]);
-                    }
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              for (unsigned int m=0; m<dim; ++m)
+                                result[fe_to_real[comp_k]][j][l][m] += (third[k][j][l][m]
+                                                                        * data.local_dof_values[k]);
+                      }
 
-                  // never touch any data for j=dim in case dim<spacedim, so
-                  // it will always be zero as it was initialized
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<dim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          jacobian_2nd_derivatives[point][i][j][l][m] = result[i][j][l][m];
-                }
-            }
-        }
-    }
+                    // never touch any data for j=dim in case dim<spacedim, so
+                    // it will always be zero as it was initialized
+                    for (unsigned int i=0; i<spacedim; ++i)
+                      for (unsigned int j=0; j<dim; ++j)
+                        for (unsigned int l=0; l<dim; ++l)
+                          for (unsigned int m=0; m<dim; ++m)
+                            jacobian_2nd_derivatives[point][i][j][l][m] = result[i][j][l][m];
+                  }
+              }
+          }
+      }
 
-    /**
-     * Update the third derivative of the transformation from unit to real cell,
-     * the Jacobian hessians, pushed forward to the real cell coordinates.
-     *
-     * Skip the computation if possible as indicated by the first argument.
-     */
-    template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-    void
-    maybe_update_jacobian_pushed_forward_2nd_derivatives (
-      const CellSimilarity::Similarity              cell_similarity,
-      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-      const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-      const FiniteElement<dim, spacedim>           &fe,
-      const ComponentMask                          &fe_mask,
-      const std::vector<unsigned int>              &fe_to_real,
-      std::vector<Tensor<4,spacedim> >             &jacobian_pushed_forward_2nd_derivatives )
-    {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
-        {
-          const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size();
+      /**
+       * Update the third derivative of the transformation from unit to real cell,
+       * the Jacobian hessians, pushed forward to the real cell coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_pushed_forward_2nd_derivatives
+      (const CellSimilarity::Similarity              cell_similarity,
+       const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+       const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+       const FiniteElement<dim, spacedim>           &fe,
+       const ComponentMask                          &fe_mask,
+       const std::vector<unsigned int>              &fe_to_real,
+       std::vector<Tensor<4,spacedim> >             &jacobian_pushed_forward_2nd_derivatives )
+      {
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
+          {
+            const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size();
 
-          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);
+            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);
 
-                  DerivativeForm<3,dim,spacedim> result;
+                    DerivativeForm<3,dim,spacedim> result;
 
-                  for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                    {
-                      unsigned int comp_k = fe.system_to_component_index(k).first;
-                      if (fe_mask[comp_k])
-                        for (unsigned int j=0; j<dim; ++j)
-                          for (unsigned int l=0; l<dim; ++l)
-                            for (unsigned int m=0; m<dim; ++m)
-                              result[fe_to_real[comp_k]][j][l][m] += (third[k][j][l][m]
-                                                                      * data.local_dof_values[k]);
-                    }
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              for (unsigned int m=0; m<dim; ++m)
+                                result[fe_to_real[comp_k]][j][l][m] += (third[k][j][l][m]
+                                                                        * data.local_dof_values[k]);
+                      }
 
-                  // push forward the j-coordinate
-                  for (unsigned int i=0; i<spacedim; ++i)
-                    for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<dim; ++l)
-                        for (unsigned int m=0; m<dim; ++m)
-                          {
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              = result[i][0][l][m]*
-                                data.covariant[point][j][0];
-                            for (unsigned int jr=1; jr<dim; ++jr)
+                    // 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][jr][l][m]*
-                                 data.covariant[point][j][jr];
-                          }
+                                = 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)
+                    // 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][lr][m]*
-                                 data.covariant[point][l][lr];
-                          }
+                                = 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)
+                    // 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][mr]*
-                                 data.covariant[point][m][mr];
-                          }
-                }
-            }
-        }
-    }
-  }
+                                = tmp[i][j][l][0]*
+                                  data.covariant[point][m][0];
+                              for (unsigned int mr=1; mr<dim; ++mr)
+                                jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
+                                += tmp[i][j][l][mr]*
+                                   data.covariant[point][m][mr];
+                            }
+                  }
+              }
+          }
+      }
 
-  /**
-   * Update the fourth derivative of the transformation from unit to real
-   * cell, the Jacobian hessian gradients.
-   *
-   * Skip the computation if possible as indicated by the first argument.
-   */
-  template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-  void
-  maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity              cell_similarity,
-                                         const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-                                         const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                                         const FiniteElement<dim, spacedim>           &fe,
-                                         const ComponentMask                          &fe_mask,
-                                         const std::vector<unsigned int>              &fe_to_real,
-                                         std::vector<DerivativeForm<4,dim,spacedim> > &jacobian_3rd_derivatives)
-  {
-    const UpdateFlags update_flags = data.update_each;
-    if (update_flags & update_jacobian_3rd_derivatives)
+      /**
+       * Update the fourth derivative of the transformation from unit to real
+       * cell, the Jacobian hessian gradients.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity              cell_similarity,
+                                             const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+                                             const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                                             const FiniteElement<dim, spacedim>           &fe,
+                                             const ComponentMask                          &fe_mask,
+                                             const std::vector<unsigned int>              &fe_to_real,
+                                             std::vector<DerivativeForm<4,dim,spacedim> > &jacobian_3rd_derivatives)
       {
-        const unsigned int n_q_points = jacobian_3rd_derivatives.size();
-
-        if (cell_similarity != CellSimilarity::translation)
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_3rd_derivatives)
           {
-            for (unsigned int point=0; point<n_q_points; ++point)
+            const unsigned int n_q_points = jacobian_3rd_derivatives.size();
+
+            if (cell_similarity != CellSimilarity::translation)
               {
-                const Tensor<4,dim> *fourth =
-                  &data.fourth_derivative(point+data_set, 0);
+                for (unsigned int point=0; point<n_q_points; ++point)
+                  {
+                    const Tensor<4,dim> *fourth =
+                      &data.fourth_derivative(point+data_set, 0);
 
-                DerivativeForm<4,dim,spacedim> result;
+                    DerivativeForm<4,dim,spacedim> result;
 
-                for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                  {
-                    unsigned int comp_k = fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              for (unsigned int m=0; m<dim; ++m)
+                                for (unsigned int n=0; n<dim; ++n)
+                                  result[fe_to_real[comp_k]][j][l][m][n] += (fourth[k][j][l][m][n]
+                                                                             * data.local_dof_values[k]);
+                      }
+
+                    // never touch any data for j,l,m,n=dim in case dim<spacedim, so
+                    // it will always be zero as it was initialized
+                    for (unsigned int i=0; i<spacedim; ++i)
                       for (unsigned int j=0; j<dim; ++j)
                         for (unsigned int l=0; l<dim; ++l)
                           for (unsigned int m=0; m<dim; ++m)
                             for (unsigned int n=0; n<dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n] += (fourth[k][j][l][m][n]
-                                                                         * data.local_dof_values[k]);
+                              jacobian_3rd_derivatives[point][i][j][l][m][n] = result[i][j][l][m][n];
                   }
-
-                // never touch any data for j,l,m,n=dim in case dim<spacedim, so
-                // it will always be zero as it was initialized
-                for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<dim; ++j)
-                    for (unsigned int l=0; l<dim; ++l)
-                      for (unsigned int m=0; m<dim; ++m)
-                        for (unsigned int n=0; n<dim; ++n)
-                          jacobian_3rd_derivatives[point][i][j][l][m][n] = result[i][j][l][m][n];
               }
           }
       }
-  }
 
-  /**
-   * Update the fourth derivative of the transformation from unit to real cell,
-   * the Jacobian hessian gradients, pushed forward to the real cell
-   * coordinates.
-   *
-   * Skip the computation if possible as indicated by the first argument.
-   */
-  template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-  void
-  maybe_update_jacobian_pushed_forward_3rd_derivatives (
-    const CellSimilarity::Similarity              cell_similarity,
-    const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-    const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-    const FiniteElement<dim, spacedim>           &fe,
-    const ComponentMask                          &fe_mask,
-    const std::vector<unsigned int>              &fe_to_real,
-    std::vector<Tensor<5,spacedim> >             &jacobian_pushed_forward_3rd_derivatives )
-  {
-    const UpdateFlags update_flags = data.update_each;
-    if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
+      /**
+       * Update the fourth derivative of the transformation from unit to real cell,
+       * the Jacobian hessian gradients, pushed forward to the real cell
+       * coordinates.
+       *
+       * Skip the computation if possible as indicated by the first argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_update_jacobian_pushed_forward_3rd_derivatives
+      (const CellSimilarity::Similarity              cell_similarity,
+       const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+       const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+       const FiniteElement<dim, spacedim>           &fe,
+       const ComponentMask                          &fe_mask,
+       const std::vector<unsigned int>              &fe_to_real,
+       std::vector<Tensor<5,spacedim> >             &jacobian_pushed_forward_3rd_derivatives )
       {
-        const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size();
-
-        if (cell_similarity != CellSimilarity::translation)
+        const UpdateFlags update_flags = data.update_each;
+        if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
           {
-            double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
-            for (unsigned int point=0; point<n_q_points; ++point)
+            const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size();
+
+            if (cell_similarity != CellSimilarity::translation)
               {
-                const Tensor<4,dim> *fourth =
-                  &data.fourth_derivative(point+data_set, 0);
+                double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
+                for (unsigned int point=0; point<n_q_points; ++point)
+                  {
+                    const Tensor<4,dim> *fourth =
+                      &data.fourth_derivative(point+data_set, 0);
 
-                DerivativeForm<4,dim,spacedim> result;
+                    DerivativeForm<4,dim,spacedim> result;
 
-                for (unsigned int k=0; k<data.n_shape_functions; ++k)
-                  {
-                    unsigned int comp_k = fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j=0; j<dim; ++j)
+                    for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                      {
+                        unsigned int comp_k = fe.system_to_component_index(k).first;
+                        if (fe_mask[comp_k])
+                          for (unsigned int j=0; j<dim; ++j)
+                            for (unsigned int l=0; l<dim; ++l)
+                              for (unsigned int m=0; m<dim; ++m)
+                                for (unsigned int n=0; n<dim; ++n)
+                                  result[fe_to_real[comp_k]][j][l][m][n]
+                                  += (fourth[k][j][l][m][n]
+                                      * data.local_dof_values[k]);
+                      }
+
+                    // push-forward the j-coordinate
+                    for (unsigned int i=0; i<spacedim; ++i)
+                      for (unsigned int j=0; j<spacedim; ++j)
                         for (unsigned int l=0; l<dim; ++l)
                           for (unsigned int m=0; m<dim; ++m)
                             for (unsigned int n=0; n<dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n]
-                              += (fourth[k][j][l][m][n]
-                                  * data.local_dof_values[k]);
+                              {
+                                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];
+                              }
                   }
-
-                // push-forward the j-coordinate
-                for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<spacedim; ++j)
-                    for (unsigned int l=0; l<dim; ++l)
-                      for (unsigned int m=0; m<dim; ++m)
-                        for (unsigned int n=0; n<dim; ++n)
-                          {
-                            tmp[i][j][l][m][n] = result[i][0][l][m][n] *
-                                                 data.covariant[point][j][0];
-                            for (unsigned int jr=1; jr<dim; ++jr)
-                              tmp[i][j][l][m][n] += result[i][jr][l][m][n] *
-                                                    data.covariant[point][j][jr];
-                          }
-
-                // push-forward the l-coordinate
-                for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<spacedim; ++j)
-                    for (unsigned int l=0; l<spacedim; ++l)
-                      for (unsigned int m=0; m<dim; ++m)
-                        for (unsigned int n=0; n<dim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                              = tmp[i][j][0][m][n] *
-                                data.covariant[point][l][0];
-                            for (unsigned int lr=1; lr<dim; ++lr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                              += tmp[i][j][lr][m][n] *
-                                 data.covariant[point][l][lr];
-                          }
-
-                // push-forward the m-coordinate
-                for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<spacedim; ++j)
-                    for (unsigned int l=0; l<spacedim; ++l)
-                      for (unsigned int m=0; m<spacedim; ++m)
-                        for (unsigned int n=0; n<dim; ++n)
-                          {
-                            tmp[i][j][l][m][n]
-                              = jacobian_pushed_forward_3rd_derivatives[point][i][j][l][0][n] *
-                                data.covariant[point][m][0];
-                            for (unsigned int mr=1; mr<dim; ++mr)
-                              tmp[i][j][l][m][n]
-                              += jacobian_pushed_forward_3rd_derivatives[point][i][j][l][mr][n] *
-                                 data.covariant[point][m][mr];
-                          }
-
-                // push-forward the n-coordinate
-                for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<spacedim; ++j)
-                    for (unsigned int l=0; l<spacedim; ++l)
-                      for (unsigned int m=0; m<spacedim; ++m)
-                        for (unsigned int n=0; n<spacedim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                              = tmp[i][j][l][m][0] *
-                                data.covariant[point][n][0];
-                            for (unsigned int nr=1; nr<dim; ++nr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                              += tmp[i][j][l][m][nr] *
-                                 data.covariant[point][n][nr];
-                          }
               }
           }
       }
-  }
-
 
-  /**
-   * Depending on what information is called for in the update flags of the
-   * @p data object, compute the various pieces of information that is
-   * required by the fill_fe_face_values() and fill_fe_subface_values()
-   * functions.  This function simply unifies the work that would be done by
-   * those two functions.
-   *
-   * The resulting data is put into the @p output_data argument.
-   */
-  template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-  void
-  maybe_compute_face_data (const dealii::Mapping<dim,spacedim> &mapping,
-                           const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
-                           const unsigned int               face_no,
-                           const unsigned int               subface_no,
-                           const std::vector<double>        &weights,
-                           const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
-                           internal::FEValues::MappingRelatedData<dim,spacedim>         &output_data)
-  {
-    const UpdateFlags update_flags = data.update_each;
 
-    if (update_flags & update_boundary_forms)
+      /**
+       * Depending on what information is called for in the update flags of the
+       * @p data object, compute the various pieces of information that is
+       * required by the fill_fe_face_values() and fill_fe_subface_values()
+       * functions.  This function simply unifies the work that would be done by
+       * those two functions.
+       *
+       * The resulting data is put into the @p output_data argument.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      maybe_compute_face_data (const dealii::Mapping<dim,spacedim> &mapping,
+                               const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
+                               const unsigned int               face_no,
+                               const unsigned int               subface_no,
+                               const std::vector<double>        &weights,
+                               const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &data,
+                               internal::FEValues::MappingRelatedData<dim,spacedim>         &output_data)
       {
-        const unsigned int n_q_points = output_data.boundary_forms.size();
-        if (update_flags & update_normal_vectors)
-          AssertDimension (output_data.normal_vectors.size(), n_q_points);
-        if (update_flags & update_JxW_values)
-          AssertDimension (output_data.JxW_values.size(), n_q_points);
-
-        // map the unit tangentials to the real cell. checking for d!=dim-1
-        // eliminates compiler warnings regarding unsigned int expressions <
-        // 0.
-        for (unsigned int d=0; d!=dim-1; ++d)
-          {
-            Assert (face_no+GeometryInfo<dim>::faces_per_cell*d <
-                    data.unit_tangentials.size(),
-                    ExcInternalError());
-            Assert (data.aux[d].size() <=
-                    data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
-                    ExcInternalError());
-
-            mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
-                               mapping_contravariant,
-                               data,
-                               make_array_view(data.aux[d]));
-          }
+        const UpdateFlags update_flags = data.update_each;
 
-        // if dim==spacedim, we can use the unit tangentials to compute the
-        // boundary form by simply taking the cross product
-        if (dim == spacedim)
+        if (update_flags & update_boundary_forms)
           {
-            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)
+            const unsigned int n_q_points = output_data.boundary_forms.size();
+            if (update_flags & update_normal_vectors)
+              AssertDimension (output_data.normal_vectors.size(), n_q_points);
+            if (update_flags & update_JxW_values)
+              AssertDimension (output_data.JxW_values.size(), n_q_points);
+
+            // map the unit tangentials to the real cell. checking for d!=dim-1
+            // eliminates compiler warnings regarding unsigned int expressions <
+            // 0.
+            for (unsigned int d=0; d!=dim-1; ++d)
               {
-                if (dim==1)
+                Assert (face_no+GeometryInfo<dim>::faces_per_cell*d <
+                        data.unit_tangentials.size(),
+                        ExcInternalError());
+                Assert (data.aux[d].size() <=
+                        data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
+                        ExcInternalError());
+
+                mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
+                                   mapping_contravariant,
+                                   data,
+                                   make_array_view(data.aux[d]));
+              }
+
+            // if dim==spacedim, we can use the unit tangentials to compute the
+            // boundary form by simply taking the cross product
+            if (dim == spacedim)
+              {
+                for (unsigned int i=0; i<n_q_points; ++i)
+                  switch (dim)
+                    {
+                    case 1:
+                      // in 1d, we don't have access to any of the data.aux
+                      // fields (because it has only dim-1 components), but we
+                      // can still compute the boundary form by simply looking
+                      // at the number of the face
+                      output_data.boundary_forms[i][0] = (face_no == 0 ?
+                                                          -1 : +1);
+                      break;
+                    case 2:
+                      output_data.boundary_forms[i] = cross_product_2d(data.aux[0][i]);
+                      break;
+                    case 3:
+                      output_data.boundary_forms[i] =
+                        cross_product_3d(data.aux[0][i], data.aux[1][i]);
+                      break;
+                    default:
+                      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)
                   {
-                    // 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==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();
+                    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();
+                        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);
-                  }
+                        // then compute the face normal from the face tangent
+                        // and the cell normal:
+                        output_data.boundary_forms[point] =
+                          cross_product_3d(data.aux[0][point], cell_normal);
+                      }
 
+                  }
               }
-          }
 
-        if (update_flags & (update_normal_vectors | update_JxW_values))
-          for (unsigned int i=0; i<output_data.boundary_forms.size(); ++i)
-            {
-              if (update_flags & update_JxW_values)
+            if (update_flags & (update_normal_vectors | 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() * weights[i];
-
-                  if (subface_no != numbers::invalid_unsigned_int)
+                  if (update_flags & update_JxW_values)
                     {
-                      const double area_ratio=GeometryInfo<dim>::subface_ratio(
-                                                cell->subface_case(face_no), subface_no);
-                      output_data.JxW_values[i] *= area_ratio;
+                      output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * weights[i];
+
+                      if (subface_no != numbers::invalid_unsigned_int)
+                        {
+                          const double area_ratio=GeometryInfo<dim>::subface_ratio(
+                                                    cell->subface_case(face_no), subface_no);
+                          output_data.JxW_values[i] *= area_ratio;
+                        }
                     }
-                }
 
-              if (update_flags & update_normal_vectors)
-                output_data.normal_vectors[i] = Point<spacedim>(output_data.boundary_forms[i] / output_data.boundary_forms[i].norm());
-            }
+                  if (update_flags & update_normal_vectors)
+                    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_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();
+            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 MappingFEField::fill_fe_face_values() and
-   * MappingFEField::fill_fe_subface_values() in a generic way, using the
-   * 'data_set' to differentiate whether we will work on a face (and if so,
-   * which one) or subface.
-   */
-  template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
-  void
-  do_fill_fe_face_values (const dealii::Mapping<dim,spacedim>                               &mapping,
-                          const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
-                          const unsigned int                                                 face_no,
-                          const unsigned int                                                 subface_no,
-                          const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
-                          const Quadrature<dim-1>                                           &quadrature,
-                          const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData    &data,
-                          const FiniteElement<dim, spacedim>                                &fe,
-                          const ComponentMask                                               &fe_mask,
-                          const std::vector<unsigned int>                                   &fe_to_real,
-                          internal::FEValues::MappingRelatedData<dim,spacedim>              &output_data)
-  {
-    maybe_compute_q_points<dim,spacedim,VectorType,DoFHandlerType>
-    (data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.quadrature_points);
-
-    maybe_update_Jacobians<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real);
-
-    maybe_update_jacobian_grads<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_grads);
-
-    maybe_update_jacobian_pushed_forward_grads<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_pushed_forward_grads);
-
-    maybe_update_jacobian_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_2nd_derivatives);
-
-    maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_pushed_forward_2nd_derivatives);
-
-    maybe_update_jacobian_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_3rd_derivatives);
-
-    maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
-    (CellSimilarity::none,
-     data_set,
-     data,
-     fe, fe_mask, fe_to_real,
-     output_data.jacobian_pushed_forward_3rd_derivatives);
-
-    maybe_compute_face_data<dim,spacedim,VectorType,DoFHandlerType>
-    (mapping,
-     cell, face_no, subface_no,
-     quadrature.get_weights(), data,
-     output_data);
+      /**
+       * Do the work of MappingFEField::fill_fe_face_values() and
+       * MappingFEField::fill_fe_subface_values() in a generic way, using the
+       * 'data_set' to differentiate whether we will work on a face (and if so,
+       * which one) or subface.
+       */
+      template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+      void
+      do_fill_fe_face_values (const dealii::Mapping<dim,spacedim>                               &mapping,
+                              const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
+                              const unsigned int                                                 face_no,
+                              const unsigned int                                                 subface_no,
+                              const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
+                              const Quadrature<dim-1>                                           &quadrature,
+                              const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData    &data,
+                              const FiniteElement<dim, spacedim>                                &fe,
+                              const ComponentMask                                               &fe_mask,
+                              const std::vector<unsigned int>                                   &fe_to_real,
+                              internal::FEValues::MappingRelatedData<dim,spacedim>              &output_data)
+      {
+        maybe_compute_q_points<dim,spacedim,VectorType,DoFHandlerType>
+        (data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.quadrature_points);
+
+        maybe_update_Jacobians<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real);
+
+        maybe_update_jacobian_grads<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_grads);
+
+        maybe_update_jacobian_pushed_forward_grads<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_pushed_forward_grads);
+
+        maybe_update_jacobian_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_2nd_derivatives);
+
+        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_pushed_forward_2nd_derivatives);
+
+        maybe_update_jacobian_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_3rd_derivatives);
+
+        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+        (CellSimilarity::none,
+         data_set,
+         data,
+         fe, fe_mask, fe_to_real,
+         output_data.jacobian_pushed_forward_3rd_derivatives);
+
+        maybe_compute_face_data<dim,spacedim,VectorType,DoFHandlerType>
+        (mapping,
+         cell, face_no, subface_no,
+         quadrature.get_weights(), data,
+         output_data);
+      }
+    }
   }
 }
 
@@ -1318,12 +1321,12 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
   update_internal_dofs(cell, data);
 
-  internal::maybe_compute_q_points<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_compute_q_points<dim,spacedim,VectorType,DoFHandlerType>
   (QProjector<dim>::DataSetDescriptor::cell (),
    data, *fe, fe_mask, fe_to_real,
    output_data.quadrature_points);
 
-  internal::maybe_update_Jacobians<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_Jacobians<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell (),
    data, *fe, fe_mask, fe_to_real);
@@ -1427,35 +1430,35 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
     }
 
   // calculate derivatives of the Jacobians
-  internal::maybe_update_jacobian_grads<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_grads<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
    output_data.jacobian_grads);
 
   // calculate derivatives of the Jacobians pushed forward to real cell coordinates
-  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_pushed_forward_grads<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
    output_data.jacobian_pushed_forward_grads);
 
   // calculate hessians of the Jacobians
-  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
    output_data.jacobian_2nd_derivatives);
 
   // calculate hessians of the Jacobians pushed forward to real cell coordinates
-  internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
    output_data.jacobian_pushed_forward_2nd_derivatives);
 
   // calculate gradients of the hessians of the Jacobians
-  internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
@@ -1463,7 +1466,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
   // calculate gradients of the hessians of the Jacobians pushed forward to real
   // cell coordinates
-  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim,VectorType,DoFHandlerType>
   (cell_similarity,
    QProjector<dim>::DataSetDescriptor::cell(),
    data, *fe, fe_mask, fe_to_real,
@@ -1491,7 +1494,7 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
 
   update_internal_dofs(cell, data);
 
-  internal::do_fill_fe_face_values<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::do_fill_fe_face_values<dim,spacedim,VectorType,DoFHandlerType>
   (*this,
    cell, face_no, numbers::invalid_unsigned_int,
    QProjector<dim>::DataSetDescriptor::
@@ -1525,7 +1528,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 
   update_internal_dofs(cell, data);
 
-  internal::do_fill_fe_face_values<dim,spacedim,VectorType,DoFHandlerType>
+  internal::MappingFEField::do_fill_fe_face_values<dim,spacedim,VectorType,DoFHandlerType>
   (*this,
    cell, face_no, numbers::invalid_unsigned_int,
    QProjector<dim>::DataSetDescriptor::
@@ -1542,101 +1545,107 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 }
 
 
-namespace
+namespace internal
 {
-  template <int dim, int spacedim, int rank, typename VectorType, typename DoFHandlerType>
-  void
-  transform_fields(const ArrayView<const Tensor<rank,dim> >                &input,
-                   const MappingType                                        mapping_type,
-                   const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-                   const ArrayView<Tensor<rank,spacedim> >                 &output)
+  namespace MappingFEField
   {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData
-    &data = static_cast<const typename MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &>(mapping_data);
-
-    switch (mapping_type)
-      {
-      case mapping_contravariant:
+    namespace
+    {
+      template <int dim, int spacedim, int rank, typename VectorType, typename DoFHandlerType>
+      void
+      transform_fields(const ArrayView<const Tensor<rank,dim> >                &input,
+                       const MappingType                                        mapping_type,
+                       const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+                       const ArrayView<Tensor<rank,spacedim> >                 &output)
       {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData
+        &data = static_cast<const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.contravariant[i], input[i]);
+        switch (mapping_type)
+          {
+          case mapping_contravariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
 
-        return;
-      }
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.contravariant[i], input[i]);
 
-      case mapping_piola:
-      {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
-        Assert (data.update_each & update_volume_elements,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
-        Assert (rank==1, ExcMessage("Only for rank 1"));
-        for (unsigned int i=0; i<output.size(); ++i)
+            return;
+          }
+
+          case mapping_piola:
           {
-            output[i] = apply_transformation(data.contravariant[i], input[i]);
-            output[i] /= data.volume_elements[i];
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+            Assert (data.update_each & update_volume_elements,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+            Assert (rank==1, ExcMessage("Only for rank 1"));
+            for (unsigned int i=0; i<output.size(); ++i)
+              {
+                output[i] = apply_transformation(data.contravariant[i], input[i]);
+                output[i] /= data.volume_elements[i];
+              }
+            return;
           }
-        return;
-      }
 
 
-      //We still allow this operation as in the
-      //reference cell Derivatives are Tensor
-      //rather than DerivativeForm
-      case mapping_covariant:
-      {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+          //We still allow this operation as in the
+          //reference cell Derivatives are Tensor
+          //rather than DerivativeForm
+          case mapping_covariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.covariant[i], input[i]);
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.covariant[i], input[i]);
 
-        return;
-      }
+            return;
+          }
 
-      default:
-        Assert(false, ExcNotImplemented());
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
-  }
 
 
-  template <int dim, int spacedim, int rank, typename VectorType, typename DoFHandlerType>
-  void
-  transform_differential_forms
-  (const ArrayView<const DerivativeForm<rank, dim,spacedim> >  &input,
-   const MappingType                                            mapping_type,
-   const typename Mapping<dim,spacedim>::InternalDataBase      &mapping_data,
-   const ArrayView<Tensor<rank+1, spacedim> >                  &output)
-  {
+      template <int dim, int spacedim, int rank, typename VectorType, typename DoFHandlerType>
+      void
+      transform_differential_forms
+      (const ArrayView<const DerivativeForm<rank, dim,spacedim> >  &input,
+       const MappingType                                            mapping_type,
+       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 MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData
-    &data = static_cast<const typename MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &>(mapping_data);
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData
+        &data = static_cast<const typename dealii::MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData &>(mapping_data);
 
-    switch (mapping_type)
-      {
-      case mapping_covariant:
-      {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+        switch (mapping_type)
+          {
+          case mapping_covariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.covariant[i], input[i]);
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.covariant[i], input[i]);
 
-        return;
-      }
-      default:
-        Assert(false, ExcNotImplemented());
-      }
+            return;
+          }
+          default:
+            Assert(false, ExcNotImplemented());
+          }
 
+      }
+    }
   }
 }
 
@@ -1652,7 +1661,7 @@ transform (const ArrayView<const Tensor<1,dim> >                  &input,
 {
   AssertDimension (input.size(), output.size());
 
-  transform_fields<dim,spacedim,1,VectorType,DoFHandlerType>(input, mapping_type, mapping_data, output);
+  internal::MappingFEField::transform_fields<dim,spacedim,1,VectorType,DoFHandlerType>(input, mapping_type, mapping_data, output);
 }
 
 
@@ -1667,7 +1676,7 @@ transform (const ArrayView<const DerivativeForm<1, dim, spacedim> > &input,
 {
   AssertDimension (input.size(), output.size());
 
-  transform_differential_forms<dim,spacedim,1,VectorType,DoFHandlerType>(input, mapping_type, mapping_data, output);
+  internal::MappingFEField::transform_differential_forms<dim,spacedim,1,VectorType,DoFHandlerType>(input, mapping_type, mapping_data, output);
 }
 
 
index 2ac37d2cb988651377f8be7f622c61064c8ecc30..8949115094496fbc64451e17bd81284c64331b1b 100644 (file)
@@ -349,121 +349,126 @@ MappingManifold<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
 
 namespace internal
 {
-  namespace
+  namespace MappingManifold
   {
-    /**
-     * 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::MappingManifold<dim,spacedim>::InternalData      &data,
-                            std::vector<Point<spacedim> >                                     &quadrature_points)
+    namespace
     {
-      const UpdateFlags update_flags = data.update_each;
+      /**
+       * 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::MappingManifold<dim,spacedim>::InternalData &data,
+       std::vector<Point<spacedim> >                                      &quadrature_points)
+      {
+        const UpdateFlags update_flags = data.update_each;
 
-      AssertDimension(data.vertices.size(), GeometryInfo<dim>::vertices_per_cell);
+        AssertDimension(data.vertices.size(), GeometryInfo<dim>::vertices_per_cell);
 
-      if (update_flags & update_quadrature_points)
-        {
-          for (unsigned int point=0; point<quadrature_points.size(); ++point)
-            {
-              quadrature_points[point] = data.manifold->
-                                         get_new_point(data.vertices,
-                                                       data.cell_manifold_quadrature_weights[point+data_set]);
-            }
-        }
-    }
+        if (update_flags & update_quadrature_points)
+          {
+            for (unsigned int point=0; point<quadrature_points.size(); ++point)
+              {
+                quadrature_points[point] = data.manifold->
+                                           get_new_point(data.vertices,
+                                                         data.cell_manifold_quadrature_weights[point+data_set]);
+              }
+          }
+      }
 
 
 
-    /**
-     * 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.
-     */
-    template <int dim, int spacedim>
-    void
-    maybe_update_Jacobians (const typename dealii::QProjector<dim>::DataSetDescriptor          data_set,
-                            const typename dealii::MappingManifold<dim,spacedim>::InternalData      &data)
-    {
-      const UpdateFlags update_flags = data.update_each;
+      /**
+       * 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.
+       */
+      template <int dim, int spacedim>
+      void
+      maybe_update_Jacobians
+      (const typename dealii::QProjector<dim>::DataSetDescriptor           data_set,
+       const typename dealii::MappingManifold<dim,spacedim>::InternalData &data)
+      {
+        const UpdateFlags update_flags = data.update_each;
 
-      if (update_flags & update_contravariant_transformation)
-        {
-          const unsigned int n_q_points = data.contravariant.size();
+        if (update_flags & update_contravariant_transformation)
+          {
+            const unsigned int n_q_points = data.contravariant.size();
 
-          std::fill(data.contravariant.begin(), data.contravariant.end(),
-                    DerivativeForm<1,dim,spacedim>());
+            std::fill(data.contravariant.begin(), data.contravariant.end(),
+                      DerivativeForm<1,dim,spacedim>());
 
-          AssertDimension(GeometryInfo<dim>::vertices_per_cell,
-                          data.vertices.size());
-          for (unsigned int point=0; point<n_q_points; ++point)
-            {
-              // Start by figuring out how to compute the direction in
-              // the reference space:
-              const Point<dim> &p = data.quad.point(point+data_set);
-
-              // And get its image on the manifold:
-              const Point<spacedim> P = data.manifold->
-                                        get_new_point(data.vertices,
-                                                      data.cell_manifold_quadrature_weights[point+data_set]);
-
-              // To compute the Jacobian, we choose dim points aligned
-              // with the dim reference axes, which are still in the
-              // given cell, and ask for the tangent vector in these
-              // directions. Choosing the points is somewhat arbitrary,
-              // so we try to be smart and we pick points which are
-              // on the opposite quadrant w.r.t. the evaluation
-              // point.
-              for (unsigned int i=0; i<dim; ++i)
-                {
-                  const Point<dim> ei = Point<dim>::unit_vector(i);
-                  const double pi = p[i];
-                  Assert(pi >=0 && pi <= 1.0,
-                         ExcInternalError("Was expecting a quadrature point "
-                                          "inside the unit reference element."));
+            AssertDimension(GeometryInfo<dim>::vertices_per_cell,
+                            data.vertices.size());
+            for (unsigned int point=0; point<n_q_points; ++point)
+              {
+                // Start by figuring out how to compute the direction in
+                // the reference space:
+                const Point<dim> &p = data.quad.point(point+data_set);
+
+                // And get its image on the manifold:
+                const Point<spacedim> P = data.manifold->
+                                          get_new_point(data.vertices,
+                                                        data.cell_manifold_quadrature_weights[point+data_set]);
+
+                // To compute the Jacobian, we choose dim points aligned
+                // with the dim reference axes, which are still in the
+                // given cell, and ask for the tangent vector in these
+                // directions. Choosing the points is somewhat arbitrary,
+                // so we try to be smart and we pick points which are
+                // on the opposite quadrant w.r.t. the evaluation
+                // point.
+                for (unsigned int i=0; i<dim; ++i)
+                  {
+                    const Point<dim> ei = Point<dim>::unit_vector(i);
+                    const double pi = p[i];
+                    Assert(pi >=0 && pi <= 1.0,
+                           ExcInternalError("Was expecting a quadrature point "
+                                            "inside the unit reference element."));
 
-                  // In the length L, we store also the direction sign,
-                  // which is positive, if the coordinate is < .5,
-                  const double L = pi > .5 ? -pi: 1-pi;
+                    // In the length L, we store also the direction sign,
+                    // which is positive, if the coordinate is < .5,
+                    const double L = pi > .5 ? -pi: 1-pi;
 
-                  const Point<dim> np(p + L*ei);
+                    const Point<dim> np(p + L*ei);
 
-                  // Get the weights to compute the np point in real space
-                  for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
-                    data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);
+                    // Get the weights to compute the np point in real space
+                    for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+                      data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);
 
-                  const Point<spacedim> NP=
-                    data.manifold->get_new_point(data.vertices,
-                                                 data.vertex_weights);
+                    const Point<spacedim> NP=
+                      data.manifold->get_new_point(data.vertices,
+                                                   data.vertex_weights);
 
-                  const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);
+                    const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);
 
-                  for (unsigned int d=0; d<spacedim; ++d)
-                    data.contravariant[point][d][i] = T[d]/L;
-                }
-            }
+                    for (unsigned int d=0; d<spacedim; ++d)
+                      data.contravariant[point][d][i] = T[d]/L;
+                  }
+              }
 
-          if (update_flags & update_covariant_transformation)
-            {
-              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_covariant_transformation)
+              {
+                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)
-            {
-              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();
-            }
-        }
+            if (update_flags & update_volume_elements)
+              {
+                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();
+              }
+          }
+      }
     }
   }
 }
@@ -489,12 +494,14 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   data.store_vertices(cell);
   data.manifold = &(cell->get_manifold());
 
-  internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
-                                                  data,
-                                                  output_data.quadrature_points);
+  internal::MappingManifold::maybe_compute_q_points<dim,spacedim>
+  (QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.quadrature_points);
 
-  internal::maybe_update_Jacobians<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
-                                                  data);
+  internal::MappingManifold::maybe_update_Jacobians<dim,spacedim>
+  (QProjector<dim>::DataSetDescriptor::cell (),
+   data);
 
   const UpdateFlags update_flags = data.update_each;
   const std::vector<double> &weights=quadrature.get_weights();
@@ -616,580 +623,582 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
 namespace internal
 {
-  namespace
+  namespace MappingManifold
   {
-    /**
-     * 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::MappingManifold<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 std::vector<double>        &weights,
-                             const typename dealii::MappingManifold<dim,spacedim>::InternalData &data,
-                             internal::FEValues::MappingRelatedData<dim,spacedim>         &output_data)
+    namespace
     {
-      const UpdateFlags update_flags = data.update_each;
-
-      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);
-
-          // 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+GeometryInfo<dim>::faces_per_cell*d <
-                      data.unit_tangentials.size(),
-                      ExcInternalError());
-              Assert (data.aux[d].size() <=
-                      data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
-                      ExcInternalError());
-
-              mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
-                                 mapping_contravariant,
-                                 data,
-                                 make_array_view(data.aux[d]));
-            }
+      /**
+       * 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::MappingManifold<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 std::vector<double>                                          &weights,
+       const typename dealii::MappingManifold<dim,spacedim>::InternalData &data,
+       internal::FEValues::MappingRelatedData<dim,spacedim>               &output_data)
+      {
+        const UpdateFlags update_flags = data.update_each;
 
-          // 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);
+        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);
+
+            // 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+GeometryInfo<dim>::faces_per_cell*d <
+                        data.unit_tangentials.size(),
+                        ExcInternalError());
+                Assert (data.aux[d].size() <=
+                        data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
+                        ExcInternalError());
+
+                mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
+                                   mapping_contravariant,
+                                   data,
+                                   make_array_view(data.aux[d]));
+              }
 
-              for (unsigned int point=0; point<n_q_points; ++point)
-                {
+            // 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:
-                    {
-                      // 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();
-
+                      // 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:
-                    {
-                      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);
-
+                      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());
+                      Assert(false, ExcNotImplemented());
                     }
-                }
-            }
-
-          if (update_flags & (update_normal_vectors
-                              | update_JxW_values))
-            for (unsigned int i=0; i<output_data.boundary_forms.size(); ++i)
+              }
+            else //(dim < spacedim)
               {
-                if (update_flags & update_JxW_values)
+                // 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)
                   {
-                    output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * weights[i];
-
-                    if (subface_no!=numbers::invalid_unsigned_int)
+                    switch (dim)
                       {
-                        const double area_ratio=GeometryInfo<dim>::subface_ratio(cell->subface_case(face_no),
-                                                                                 subface_no);
-                        output_data.JxW_values[i] *= area_ratio;
-                      }
-                  }
-
-                if (update_flags & update_normal_vectors)
-                  output_data.normal_vectors[i] = Point<spacedim>(output_data.boundary_forms[i] /
-                                                                  output_data.boundary_forms[i].norm());
-              }
-
-          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();
-        }
-    }
-
+                      case 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();
 
-    /**
-     * Do the work of MappingManifold::fill_fe_face_values() and
-     * MappingManifold::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::MappingManifold<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::MappingManifold<dim,spacedim>::InternalData      &data,
-                            internal::FEValues::MappingRelatedData<dim,spacedim>              &output_data)
-    {
-      data.store_vertices(cell);
+                        break;
+                      }
 
-      data.manifold = &cell->face(face_no)->get_manifold();
+                      case 2:
+                      {
+                        const DerivativeForm<1,spacedim,dim> DX_t =
+                          data.contravariant[point].transpose();
 
-      maybe_compute_q_points<dim,spacedim> (data_set,
-                                            data,
-                                            output_data.quadrature_points);
-      maybe_update_Jacobians<dim,spacedim> (data_set,
-                                            data);
+                        Tensor<1, spacedim> cell_normal =
+                          cross_product_3d(DX_t[0], DX_t[1]);
+                        cell_normal /= cell_normal.norm();
 
-      maybe_compute_face_data (mapping,
-                               cell, face_no, subface_no, quadrature.size(),
-                               quadrature.get_weights(), data,
-                               output_data);
-    }
-  }
-}
+                        // 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);
 
+                        break;
+                      }
 
+                      default:
+                        Assert (false, ExcNotImplemented());
+                      }
+                  }
+              }
 
-template <int dim, int spacedim>
-void
-MappingManifold<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::FEValues::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 (update_flags & (update_normal_vectors
+                                | update_JxW_values))
+              for (unsigned int i=0; i<output_data.boundary_forms.size(); ++i)
+                {
+                  if (update_flags & update_JxW_values)
+                    {
+                      output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * weights[i];
+
+                      if (subface_no!=numbers::invalid_unsigned_int)
+                        {
+                          const double area_ratio=GeometryInfo<dim>::subface_ratio(cell->subface_case(face_no),
+                                                                                   subface_no);
+                          output_data.JxW_values[i] *= area_ratio;
+                        }
+                    }
 
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, numbers::invalid_unsigned_int,
-                                    QProjector<dim>::DataSetDescriptor::face (face_no,
-                                        cell->face_orientation(face_no),
-                                        cell->face_flip(face_no),
-                                        cell->face_rotation(face_no),
-                                        quadrature.size()),
-                                    quadrature,
-                                    data,
-                                    output_data);
-}
+                  if (update_flags & update_normal_vectors)
+                    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();
+          }
+      }
 
-template <int dim, int spacedim>
-void
-MappingManifold<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::FEValues::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);
 
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, subface_no,
-                                    QProjector<dim>::DataSetDescriptor::subface (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);
-}
+      /**
+       * Do the work of MappingManifold::fill_fe_face_values() and
+       * MappingManifold::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::MappingManifold<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::MappingManifold<dim,spacedim>::InternalData &data,
+       internal::FEValues::MappingRelatedData<dim,spacedim>               &output_data)
+      {
+        data.store_vertices(cell);
 
+        data.manifold = &cell->face(face_no)->get_manifold();
 
+        maybe_compute_q_points<dim,spacedim> (data_set,
+                                              data,
+                                              output_data.quadrature_points);
+        maybe_update_Jacobians<dim,spacedim> (data_set,
+                                              data);
 
-namespace
-{
-  template <int dim, int spacedim, int rank>
-  void
-  transform_fields(const ArrayView<const Tensor<rank,dim> >               &input,
-                   const MappingType                                       mapping_type,
-                   const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
-                   const ArrayView<Tensor<rank,spacedim> >                &output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingManifold<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
+        maybe_compute_face_data (mapping,
+                                 cell, face_no, subface_no, quadrature.size(),
+                                 quadrature.get_weights(), data,
+                                 output_data);
+      }
 
-    switch (mapping_type)
+      template <int dim, int spacedim, int rank>
+      void
+      transform_fields(const ArrayView<const Tensor<rank,dim> >               &input,
+                       const MappingType                                       mapping_type,
+                       const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                       const ArrayView<Tensor<rank,spacedim> >                &output)
       {
-      case mapping_contravariant:
-      {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingManifold<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.contravariant[i], input[i]);
+        switch (mapping_type)
+          {
+          case mapping_contravariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
 
-        return;
-      }
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.contravariant[i], input[i]);
 
-      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)
+            return;
+          }
+
+          case mapping_piola:
           {
-            output[i] = apply_transformation(data.contravariant[i], input[i]);
-            output[i] /= data.volume_elements[i];
+            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;
           }
-        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"));
+          //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]);
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.covariant[i], input[i]);
 
-        return;
-      }
+            return;
+          }
 
-      default:
-        Assert(false, ExcNotImplemented());
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
-  }
-
 
-  template <int dim, int spacedim, int rank>
-  void
-  transform_gradients(const ArrayView<const Tensor<rank,dim> >                &input,
-                      const MappingType                                        mapping_type,
-                      const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-                      const ArrayView<Tensor<rank,spacedim> >                 &output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingManifold<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-    switch (mapping_type)
+      template <int dim, int spacedim, int rank>
+      void
+      transform_gradients(const ArrayView<const Tensor<rank,dim> >                &input,
+                          const MappingType                                        mapping_type,
+                          const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+                          const ArrayView<Tensor<rank,spacedim> >                 &output)
       {
-      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"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingManifold<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
+        switch (mapping_type)
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.contravariant[i], transpose(input[i]) );
-            output[i] = apply_transformation(data.covariant[i], A.transpose() );
-          }
+          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"));
 
-        return;
-      }
+            for (unsigned int i=0; i<output.size(); ++i)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.contravariant[i], transpose(input[i]) );
+                output[i] = apply_transformation(data.covariant[i], A.transpose() );
+              }
 
-      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"));
+            return;
+          }
 
-        for (unsigned int i=0; i<output.size(); ++i)
+          case mapping_covariant_gradient:
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.covariant[i], transpose(input[i]) );
-            output[i] = apply_transformation(data.covariant[i], A.transpose() );
-          }
+            Assert (data.update_each & update_covariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+            Assert (rank==2, ExcMessage("Only for rank 2"));
 
-        return;
-      }
+            for (unsigned int i=0; i<output.size(); ++i)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.covariant[i], transpose(input[i]) );
+                output[i] = apply_transformation(data.covariant[i], A.transpose() );
+              }
 
-      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)
+            return;
+          }
+
+          case mapping_piola_gradient:
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.covariant[i], input[i] );
-            Tensor<2,spacedim> T =
-              apply_transformation(data.contravariant[i], A.transpose() );
+            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)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.covariant[i], input[i] );
+                Tensor<2,spacedim> T =
+                  apply_transformation(data.contravariant[i], A.transpose() );
 
-            output[i] = transpose(T);
-            output[i] /= data.volume_elements[i];
-          }
+                output[i] = transpose(T);
+                output[i] /= data.volume_elements[i];
+              }
 
-        return;
-      }
+            return;
+          }
 
-      default:
-        Assert(false, ExcNotImplemented());
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
-  }
 
 
 
 
-  template <int dim, int spacedim>
-  void
-  transform_hessians(const ArrayView<const Tensor<3,dim> >                  &input,
-                     const MappingType                                       mapping_type,
-                     const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
-                     const ArrayView<Tensor<3,spacedim> >                   &output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingManifold<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
-
-    switch (mapping_type)
-      {
-      case mapping_contravariant_hessian:
+      template <int dim, int spacedim>
+      void
+      transform_hessians(const ArrayView<const Tensor<3,dim> >                  &input,
+                         const MappingType                                       mapping_type,
+                         const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                         const ArrayView<Tensor<3,spacedim> >                   &output)
       {
-        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"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingManifold<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-        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)
+        switch (mapping_type)
+          {
+          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 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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            return;
+          }
 
-      case mapping_covariant_hessian:
-      {
-        Assert (data.update_each & update_covariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+          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)
+            for (unsigned int q=0; q<output.size(); ++q)
+              for (unsigned int i=0; i<spacedim; ++i)
                 {
-                  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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            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)
+          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 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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            return;
+          }
 
-      default:
-        Assert(false, ExcNotImplemented());
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
-  }
-
 
 
 
-  template <int dim, int spacedim, int rank>
-  void
-  transform_differential_forms(const ArrayView<const DerivativeForm<rank, dim,spacedim> >   &input,
-                               const MappingType                                             mapping_type,
-                               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 MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingManifold<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-    switch (mapping_type)
+      template <int dim, int spacedim, int rank>
+      void
+      transform_differential_forms(const ArrayView<const DerivativeForm<rank, dim,spacedim> >   &input,
+                                   const MappingType                                             mapping_type,
+                                   const typename Mapping<dim,spacedim>::InternalDataBase       &mapping_data,
+                                   const ArrayView<Tensor<rank+1, spacedim> >                   &output)
       {
-      case mapping_covariant:
-      {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingManifold<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingManifold<dim,spacedim>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.covariant[i], input[i]);
+        switch (mapping_type)
+          {
+          case mapping_covariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
 
-        return;
-      }
-      default:
-        Assert(false, ExcNotImplemented());
+            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>
+void
+MappingManifold<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::FEValues::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);
+
+  internal::MappingManifold::do_fill_fe_face_values
+  (*this,
+   cell, face_no, numbers::invalid_unsigned_int,
+   QProjector<dim>::DataSetDescriptor::face (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
+MappingManifold<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::FEValues::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);
+
+  internal::MappingManifold::do_fill_fe_face_values
+  (*this,
+   cell, face_no, subface_no,
+   QProjector<dim>::DataSetDescriptor::subface (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);
+}
+
+
+
 template <int dim, int spacedim>
 void
 MappingManifold<dim,spacedim>::
@@ -1198,7 +1207,7 @@ transform (const ArrayView<const Tensor<1, dim> >                  &input,
            const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
            const ArrayView<Tensor<1, spacedim> >                   &output) const
 {
-  transform_fields(input, mapping_type, mapping_data, output);
+  internal::MappingManifold::transform_fields(input, mapping_type, mapping_data, output);
 }
 
 
@@ -1211,7 +1220,7 @@ transform (const ArrayView<const DerivativeForm<1, dim,spacedim> >  &input,
            const typename Mapping<dim,spacedim>::InternalDataBase   &mapping_data,
            const ArrayView<Tensor<2, spacedim> >                    &output) const
 {
-  transform_differential_forms(input, mapping_type, mapping_data, output);
+  internal::MappingManifold::transform_differential_forms(input, mapping_type, mapping_data, output);
 }
 
 
@@ -1227,13 +1236,13 @@ transform (const ArrayView<const Tensor<2, dim> >                  &input,
   switch (mapping_type)
     {
     case mapping_contravariant:
-      transform_fields(input, mapping_type, mapping_data, output);
+      internal::MappingManifold::transform_fields(input, mapping_type, mapping_data, output);
       return;
 
     case mapping_piola_gradient:
     case mapping_contravariant_gradient:
     case mapping_covariant_gradient:
-      transform_gradients(input, mapping_type, mapping_data, output);
+      internal::MappingManifold::transform_gradients(input, mapping_type, mapping_data, output);
       return;
     default:
       Assert(false, ExcNotImplemented());
@@ -1304,7 +1313,7 @@ transform (const ArrayView<const  Tensor<3,dim> >                  &input,
     case mapping_piola_hessian:
     case mapping_contravariant_hessian:
     case mapping_covariant_hessian:
-      transform_hessians(input, mapping_type, mapping_data, output);
+      internal::MappingManifold::transform_hessians(input, mapping_type, mapping_data, output);
       return;
     default:
       Assert(false, ExcNotImplemented());
@@ -1316,4 +1325,3 @@ transform (const ArrayView<const  Tensor<3,dim> >                  &input,
 
 
 DEAL_II_NAMESPACE_CLOSE
-
index 97f3d31512eb7720ff111b04e1c69225348cc833..292af5d8b90b4f326215d1e52918c775f58df755 100644 (file)
@@ -18,6 +18,7 @@
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/qprojector.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/base/memory_consumption.h>
 #include <deal.II/lac/full_matrix.h>
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace
+namespace internal
 {
-  template <int dim>
-  std::vector<unsigned int>
-  get_dpo_vector (const unsigned int degree)
+  namespace MappingQGeneric
   {
-    std::vector<unsigned int> dpo(dim+1, 1U);
-    for (unsigned int i=1; i<dpo.size(); ++i)
-      dpo[i]=dpo[i-1]*(degree-1);
-    return dpo;
+    namespace
+    {
+      template <int dim>
+      std::vector<unsigned int>
+      get_dpo_vector (const unsigned int degree)
+      {
+        std::vector<unsigned int> dpo(dim+1, 1U);
+        for (unsigned int i=1; i<dpo.size(); ++i)
+          dpo[i]=dpo[i-1]*(degree-1);
+        return dpo;
+      }
+    }
   }
-}
 
-namespace internal
-{
   namespace MappingQ1
   {
     namespace
     {
-
       // These are left as templates on the spatial dimension (even though dim
       // == spacedim must be true for them to make sense) because templates are
       // expanded before the compiler eliminates code due to the 'if (dim ==
@@ -205,9 +208,9 @@ namespace internal
         // then also construct the mapping from lexicographic to the Qp shape function numbering
         const std::vector<unsigned int>
         renumber (FETools::
-                  lexicographic_to_hierarchic_numbering (
-                    FiniteElementData<dim> (get_dpo_vector<dim>(data.polynomial_degree), 1,
-                                            data.polynomial_degree)));
+                  lexicographic_to_hierarchic_numbering
+                  (FiniteElementData<dim> (internal::MappingQGeneric::get_dpo_vector<dim>
+                                           (data.polynomial_degree), 1, data.polynomial_degree)));
 
         std::vector<double> values;
         std::vector<Tensor<1,dim> > grads;
@@ -830,633 +833,501 @@ compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
 }
 
 
-namespace
+namespace internal
 {
-  /**
-   * Compute the <tt>support_point_weights_on_quad(hex)</tt> arrays.
-   *
-   * Called by the <tt>compute_support_point_weights_on_quad(hex)</tt> functions if the
-   * data is not yet hardcoded.
-   *
-   * For the definition of the <tt>support_point_weights_on_quad(hex)</tt> please
-   * refer to equation (8) of the `mapping' report.
-   */
-  template <int dim>
-  Table<2,double>
-  compute_laplace_vector(const unsigned int polynomial_degree)
+  namespace MappingQGeneric
   {
-    Table<2,double> lvs;
-
-    Assert(lvs.n_rows()==0, ExcInternalError());
-    Assert(dim==2 || dim==3, ExcNotImplemented());
-
-    // for degree==1, we shouldn't have to compute any support points, since all
-    // of them are on the vertices
-    Assert(polynomial_degree>1, ExcInternalError());
-
-    const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
-    const unsigned int n_outer = (dim==1) ? 2 :
-                                 ((dim==2) ?
-                                  4+4*(polynomial_degree-1) :
-                                  8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1));
-
-
-    // compute the shape gradients at the quadrature points on the unit cell
-    const QGauss<dim> quadrature(polynomial_degree+1);
-    const unsigned int n_q_points=quadrature.size();
-
-    typename MappingQGeneric<dim>::InternalData quadrature_data(polynomial_degree);
-    quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions *
-                                             n_q_points);
-    quadrature_data.compute_shape_function_values(quadrature.get_points());
+    namespace
+    {
+      /**
+       * Compute the <tt>support_point_weights_on_quad(hex)</tt> arrays.
+       *
+       * Called by the <tt>compute_support_point_weights_on_quad(hex)</tt> functions if the
+       * data is not yet hardcoded.
+       *
+       * For the definition of the <tt>support_point_weights_on_quad(hex)</tt> please
+       * refer to equation (8) of the `mapping' report.
+       */
+      template <int dim>
+      dealii::Table<2,double>
+      compute_laplace_vector(const unsigned int polynomial_degree)
+      {
+        dealii::Table<2,double> lvs;
 
-    // Compute the stiffness matrix of the inner dofs
-    FullMatrix<long double> S(n_inner);
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<n_inner; ++i)
-        for (unsigned int j=0; j<n_inner; ++j)
-          {
-            long double res = 0.;
-            for (unsigned int l=0; l<dim; ++l)
-              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
-                     (long double)quadrature_data.derivative(point, n_outer+j)[l];
+        Assert(lvs.n_rows()==0, ExcInternalError());
+        Assert(dim==2 || dim==3, ExcNotImplemented());
 
-            S(i,j) += res * (long double)quadrature.weight(point);
-          }
+        // for degree==1, we shouldn't have to compute any support points, since all
+        // of them are on the vertices
+        Assert(polynomial_degree>1, ExcInternalError());
 
-    // Compute the components of T to be the product of gradients of inner and
-    // outer shape functions.
-    FullMatrix<long double> T(n_inner, n_outer);
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<n_inner; ++i)
-        for (unsigned int k=0; k<n_outer; ++k)
-          {
-            long double res = 0.;
-            for (unsigned int l=0; l<dim; ++l)
-              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
-                     (long double)quadrature_data.derivative(point, k)[l];
+        const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
+        const unsigned int n_outer = (dim==1) ? 2 :
+                                     ((dim==2) ?
+                                      4+4*(polynomial_degree-1) :
+                                      8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1));
 
-            T(i,k) += res *(long double)quadrature.weight(point);
-          }
 
-    FullMatrix<long double> S_1(n_inner);
-    S_1.invert(S);
+        // compute the shape gradients at the quadrature points on the unit cell
+        const QGauss<dim> quadrature(polynomial_degree+1);
+        const unsigned int n_q_points=quadrature.size();
 
-    FullMatrix<long double> S_1_T(n_inner, n_outer);
+        typename dealii::MappingQGeneric<dim>::InternalData quadrature_data(polynomial_degree);
+        quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions *
+                                                 n_q_points);
+        quadrature_data.compute_shape_function_values(quadrature.get_points());
 
-    // S:=S_1*T
-    S_1.mmult(S_1_T,T);
+        // Compute the stiffness matrix of the inner dofs
+        FullMatrix<long double> S(n_inner);
+        for (unsigned int point=0; point<n_q_points; ++point)
+          for (unsigned int i=0; i<n_inner; ++i)
+            for (unsigned int j=0; j<n_inner; ++j)
+              {
+                long double res = 0.;
+                for (unsigned int l=0; l<dim; ++l)
+                  res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                         (long double)quadrature_data.derivative(point, n_outer+j)[l];
 
-    // Resize and initialize the lvs
-    lvs.reinit (n_inner, n_outer);
-    for (unsigned int i=0; i<n_inner; ++i)
-      for (unsigned int k=0; k<n_outer; ++k)
-        lvs(i,k) = -S_1_T(i,k);
+                S(i,j) += res * (long double)quadrature.weight(point);
+              }
 
-    return lvs;
-  }
+        // Compute the components of T to be the product of gradients of inner and
+        // outer shape functions.
+        FullMatrix<long double> T(n_inner, n_outer);
+        for (unsigned int point=0; point<n_q_points; ++point)
+          for (unsigned int i=0; i<n_inner; ++i)
+            for (unsigned int k=0; k<n_outer; ++k)
+              {
+                long double res = 0.;
+                for (unsigned int l=0; l<dim; ++l)
+                  res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                         (long double)quadrature_data.derivative(point, k)[l];
 
+                T(i,k) += res *(long double)quadrature.weight(point);
+              }
 
-  /**
-   * This function is needed by the constructor of
-   * <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 2 and 3.
-   *
-   * For <tt>degree<4</tt> this function sets the @p support_point_weights_on_quad to
-   * the hardcoded data. For <tt>degree>=4</tt> and MappingQ<2> this vector is
-   * computed.
-   *
-   * For the definition of the @p support_point_weights_on_quad please refer to
-   * equation (8) of the `mapping' report.
-   */
-  Table<2,double>
-  compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
-  {
-    Table<2,double> loqvs;
+        FullMatrix<long double> S_1(n_inner);
+        S_1.invert(S);
 
-    // we are asked to compute weights for interior support points, but
-    // there are no interior points if degree==1
-    if (polynomial_degree == 1)
-      return loqvs;
+        FullMatrix<long double> S_1_T(n_inner, n_outer);
 
-    const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1);
-    const unsigned int n_outer_2d=4+4*(polynomial_degree-1);
+        // S:=S_1*T
+        S_1.mmult(S_1_T,T);
 
-    // first check whether we have precomputed the values for some polynomial
-    // degree; the sizes of arrays is n_inner_2d*n_outer_2d
-    if (polynomial_degree == 2)
-      {
-        // (checked these values against the output of compute_laplace_vector
-        // again, and found they're indeed right -- just in case someone wonders
-        // where they come from -- WB)
-        static const double loqv2[1*8]
-          = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
-        Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
-                n_inner_2d * n_outer_2d,
-                ExcInternalError());
+        // Resize and initialize the lvs
+        lvs.reinit (n_inner, n_outer);
+        for (unsigned int i=0; i<n_inner; ++i)
+          for (unsigned int k=0; k<n_outer; ++k)
+            lvs(i,k) = -S_1_T(i,k);
 
-        // copy and return
-        loqvs.reinit(n_inner_2d, n_outer_2d);
-        for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
-          for (unsigned int k=0; k<n_outer_2d; ++k)
-            loqvs[unit_point][k] = loqv2[unit_point*n_outer_2d+k];
-      }
-    else
-      {
-        // not precomputed, then do so now
-        loqvs = compute_laplace_vector<2>(polynomial_degree);
+        return lvs;
       }
 
-    // the sum of weights of the points at the outer rim should be one. check
-    // this
-    for (unsigned int unit_point=0; unit_point<loqvs.n_rows(); ++unit_point)
-      Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
-                                       loqvs[unit_point].end(),0.)-1)<1e-13*polynomial_degree,
-             ExcInternalError());
-
-    return loqvs;
-  }
-
-
-
-  /**
-   * This function is needed by the constructor of <tt>MappingQ<3></tt>.
-   *
-   * For <tt>degree==2</tt> this function sets the @p support_point_weights_on_hex to
-   * the hardcoded data. For <tt>degree>2</tt> this vector is computed.
-   *
-   * For the definition of the @p support_point_weights_on_hex please refer to
-   * equation (8) of the `mapping' report.
-   */
-  Table<2,double>
-  compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
-  {
-    Table<2,double> lohvs;
-
-    // we are asked to compute weights for interior support points, but
-    // there are no interior points if degree==1
-    if (polynomial_degree == 1)
-      return lohvs;
 
-    const unsigned int n_inner = Utilities::fixed_power<3>(polynomial_degree-1);
-    const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1);
 
-    // first check whether we have precomputed the values for some polynomial
-    // degree; the sizes of arrays is n_inner_2d*n_outer_2d
-    if (polynomial_degree == 2)
+      /**
+       * This function is needed by the constructor of
+       * <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 2 and 3.
+       *
+       * For <tt>degree<4</tt> this function sets the @p support_point_weights_on_quad to
+       * the hardcoded data. For <tt>degree>=4</tt> and MappingQ<2> this vector is
+       * computed.
+       *
+       * For the definition of the @p support_point_weights_on_quad please refer to
+       * equation (8) of the `mapping' report.
+       */
+      dealii::Table<2,double>
+      compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
       {
-        static const double lohv2[26]
-          = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
-             7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192.,
-             7/192., 7/192., 7/192., 7/192.,
-             1/12., 1/12., 1/12., 1/12., 1/12., 1/12.
-            };
-
-        // copy and return
-        lohvs.reinit(n_inner, n_outer);
-        for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
-          for (unsigned int k=0; k<n_outer; ++k)
-            lohvs[unit_point][k] = lohv2[unit_point*n_outer+k];
-      }
-    else
-      {
-        // not precomputed, then do so now
-        lohvs = compute_laplace_vector<3>(polynomial_degree);
-      }
-
-    // the sum of weights of the points at the outer rim should be one. check
-    // this
-    for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
-      Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
-                                       lohvs[unit_point].end(),0.) - 1)<1e-13*polynomial_degree,
-             ExcInternalError());
-
-    return lohvs;
-  }
+        dealii::Table<2,double> loqvs;
 
-  /**
-   * This function collects the output of
-   * compute_support_point_weights_on_{quad,hex} in a single data structure.
-   */
-  std::vector<Table<2,double> >
-  compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree,
-                                                      const unsigned int dim)
-  {
-    Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
-    std::vector<Table<2,double> > output(dim);
-    if (polynomial_degree <= 1)
-      return output;
-
-    // fill the 1D interior weights
-    QGaussLobatto<1> quadrature(polynomial_degree+1);
-    output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell);
-    for (unsigned int q=0; q<polynomial_degree-1; ++q)
-      for (unsigned int i=0; i<GeometryInfo<1>::vertices_per_cell; ++i)
-        output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1),
-                                                                  i);
-
-    if (dim > 1)
-      output[1] = compute_support_point_weights_on_quad(polynomial_degree);
-
-    if (dim > 2)
-      output[2] = compute_support_point_weights_on_hex(polynomial_degree);
-
-    return output;
-  }
-
-  /**
-   * Collects all interior points for the various dimensions.
-   */
-  template <int dim>
-  Table<2,double>
-  compute_support_point_weights_cell(const unsigned int polynomial_degree)
-  {
-    Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
-    if (polynomial_degree <= 1)
-      return Table<2,double>();
+        // we are asked to compute weights for interior support points, but
+        // there are no interior points if degree==1
+        if (polynomial_degree == 1)
+          return loqvs;
 
-    QGaussLobatto<dim> quadrature(polynomial_degree+1);
-    std::vector<unsigned int> h2l(quadrature.size());
-    FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree, h2l);
+        const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1);
+        const unsigned int n_outer_2d=4+4*(polynomial_degree-1);
 
-    Table<2,double> output(quadrature.size() - GeometryInfo<dim>::vertices_per_cell,
-                           GeometryInfo<dim>::vertices_per_cell);
-    for (unsigned int q=0; q<output.size(0); ++q)
-      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-        output(q,i) = GeometryInfo<dim>::d_linear_shape_function(quadrature.point(h2l[q+GeometryInfo<dim>::vertices_per_cell]),
-                                                                 i);
+        // first check whether we have precomputed the values for some polynomial
+        // degree; the sizes of arrays is n_inner_2d*n_outer_2d
+        if (polynomial_degree == 2)
+          {
+            // (checked these values against the output of compute_laplace_vector
+            // again, and found they're indeed right -- just in case someone wonders
+            // where they come from -- WB)
+            static const double loqv2[1*8]
+              = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
+            Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
+                    n_inner_2d * n_outer_2d,
+                    ExcInternalError());
+
+            // copy and return
+            loqvs.reinit(n_inner_2d, n_outer_2d);
+            for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
+              for (unsigned int k=0; k<n_outer_2d; ++k)
+                loqvs[unit_point][k] = loqv2[unit_point*n_outer_2d+k];
+          }
+        else
+          {
+            // not precomputed, then do so now
+            loqvs = compute_laplace_vector<2>(polynomial_degree);
+          }
 
-    return output;
-  }
-}
+        // the sum of weights of the points at the outer rim should be one. check
+        // this
+        for (unsigned int unit_point=0; unit_point<loqvs.n_rows(); ++unit_point)
+          Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
+                                           loqvs[unit_point].end(),0.)-1)<1e-13*polynomial_degree,
+                 ExcInternalError());
 
+        return loqvs;
+      }
 
 
 
-template <int dim, int spacedim>
-MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
-  :
-  polynomial_degree(p),
-  line_support_points(this->polynomial_degree+1),
-  fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : nullptr),
-  support_point_weights_perimeter_to_interior (compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)),
-  support_point_weights_cell (compute_support_point_weights_cell<dim>(this->polynomial_degree))
-{
-  Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings "
-                              "with a polynomial degree greater or equal to one."));
-}
+      /**
+       * This function is needed by the constructor of <tt>MappingQ<3></tt>.
+       *
+       * For <tt>degree==2</tt> this function sets the @p support_point_weights_on_hex to
+       * the hardcoded data. For <tt>degree>2</tt> this vector is computed.
+       *
+       * For the definition of the @p support_point_weights_on_hex please refer to
+       * equation (8) of the `mapping' report.
+       */
+      dealii::Table<2,double>
+      compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
+      {
+        dealii::Table<2,double> lohvs;
 
+        // we are asked to compute weights for interior support points, but
+        // there are no interior points if degree==1
+        if (polynomial_degree == 1)
+          return lohvs;
 
+        const unsigned int n_inner = Utilities::fixed_power<3>(polynomial_degree-1);
+        const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1);
 
-template <int dim, int spacedim>
-MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
-  :
-  polynomial_degree(mapping.polynomial_degree),
-  line_support_points(mapping.line_support_points),
-  fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : nullptr),
-  support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior),
-  support_point_weights_cell (mapping.support_point_weights_cell)
-{}
+        // first check whether we have precomputed the values for some polynomial
+        // degree; the sizes of arrays is n_inner_2d*n_outer_2d
+        if (polynomial_degree == 2)
+          {
+            static const double lohv2[26]
+              = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
+                 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192.,
+                 7/192., 7/192., 7/192., 7/192.,
+                 1/12., 1/12., 1/12., 1/12., 1/12., 1/12.
+                };
+
+            // copy and return
+            lohvs.reinit(n_inner, n_outer);
+            for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
+              for (unsigned int k=0; k<n_outer; ++k)
+                lohvs[unit_point][k] = lohv2[unit_point*n_outer+k];
+          }
+        else
+          {
+            // not precomputed, then do so now
+            lohvs = compute_laplace_vector<3>(polynomial_degree);
+          }
 
+        // the sum of weights of the points at the outer rim should be one. check
+        // this
+        for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
+          Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
+                                           lohvs[unit_point].end(),0.) - 1)<1e-13*polynomial_degree,
+                 ExcInternalError());
 
+        return lohvs;
+      }
 
 
-template <int dim, int spacedim>
-Mapping<dim,spacedim> *
-MappingQGeneric<dim,spacedim>::clone () const
-{
-  return new MappingQGeneric<dim,spacedim>(*this);
-}
 
+      /**
+       * This function collects the output of
+       * compute_support_point_weights_on_{quad,hex} in a single data structure.
+       */
+      std::vector<dealii::Table<2,double> >
+      compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree,
+                                                          const unsigned int dim)
+      {
+        Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+        std::vector<dealii::Table<2,double> > output(dim);
+        if (polynomial_degree <= 1)
+          return output;
+
+        // fill the 1D interior weights
+        QGaussLobatto<1> quadrature(polynomial_degree+1);
+        output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell);
+        for (unsigned int q=0; q<polynomial_degree-1; ++q)
+          for (unsigned int i=0; i<GeometryInfo<1>::vertices_per_cell; ++i)
+            output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1),
+                                                                      i);
+
+        if (dim > 1)
+          output[1] = compute_support_point_weights_on_quad(polynomial_degree);
+
+        if (dim > 2)
+          output[2] = compute_support_point_weights_on_hex(polynomial_degree);
+
+        return output;
+      }
 
+      /**
+       * Collects all interior points for the various dimensions.
+       */
+      template <int dim>
+      dealii::Table<2,double>
+      compute_support_point_weights_cell(const unsigned int polynomial_degree)
+      {
+        Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+        if (polynomial_degree <= 1)
+          return dealii::Table<2,double>();
+
+        QGaussLobatto<dim> quadrature(polynomial_degree+1);
+        std::vector<unsigned int> h2l(quadrature.size());
+        FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree, h2l);
+
+        dealii::Table<2,double> output(quadrature.size() - GeometryInfo<dim>::vertices_per_cell,
+                                       GeometryInfo<dim>::vertices_per_cell);
+        for (unsigned int q=0; q<output.size(0); ++q)
+          for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+            output(q,i) = GeometryInfo<dim>::d_linear_shape_function(quadrature.point(h2l[q+GeometryInfo<dim>::vertices_per_cell]),
+                                                                     i);
+
+        return output;
+      }
 
 
-template <int dim, int spacedim>
-unsigned int
-MappingQGeneric<dim,spacedim>::get_degree() const
-{
-  return polynomial_degree;
-}
 
+      /**
+       * Using the relative weights of the shape functions evaluated at
+       * one point on the reference cell (and stored in data.shape_values
+       * and accessed via data.shape(0,i)) and the locations of mapping
+       * support points (stored in data.mapping_support_points), compute
+       * the mapped location of that point in real space.
+       */
+      template <int dim, int spacedim>
+      Point<spacedim>
+      compute_mapped_location_of_point
+      (const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &data)
+      {
+        AssertDimension (data.shape_values.size(),
+                         data.mapping_support_points.size());
 
+        // use now the InternalData to compute the point in real space.
+        Point<spacedim> p_real;
+        for (unsigned int i=0; i<data.mapping_support_points.size(); ++i)
+          p_real += data.mapping_support_points[i] * data.shape(0,i);
 
-template <int dim, int spacedim>
-Point<spacedim>
-MappingQGeneric<dim,spacedim>::
-transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                             const Point<dim> &p) const
-{
-  // set up the polynomial space
-  const TensorProductPolynomials<dim>
-  tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
-  Assert (tensor_pols.n() == Utilities::fixed_power<dim>(polynomial_degree+1),
-          ExcInternalError());
+        return p_real;
+      }
 
-  // then also construct the mapping from lexicographic to the Qp shape function numbering
-  const std::vector<unsigned int>
-  renumber (FETools::
-            lexicographic_to_hierarchic_numbering (
-              FiniteElementData<dim> (get_dpo_vector<dim>(polynomial_degree), 1,
-                                      polynomial_degree)));
 
-  const std::vector<Point<spacedim> > support_points
-    = this->compute_mapping_support_points(cell);
 
-  Point<spacedim> mapped_point;
-  for (unsigned int i=0; i<tensor_pols.n(); ++i)
-    mapped_point += support_points[renumber[i]] * tensor_pols.compute_value (i, p);
+      /**
+       * Implementation of transform_real_to_unit_cell for dim==spacedim
+       */
+      template <int dim>
+      Point<dim>
+      do_transform_real_to_unit_cell_internal
+      (const typename dealii::Triangulation<dim,dim>::cell_iterator &cell,
+       const Point<dim>                                             &p,
+       const Point<dim>                                             &initial_p_unit,
+       typename dealii::MappingQGeneric<dim,dim>::InternalData      &mdata)
+      {
+        const unsigned int spacedim = dim;
 
-  return mapped_point;
-}
+        const unsigned int n_shapes=mdata.shape_values.size();
+        (void)n_shapes;
+        Assert(n_shapes!=0, ExcInternalError());
+        AssertDimension (mdata.shape_derivatives.size(), n_shapes);
 
+        std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
+        AssertDimension (points.size(), n_shapes);
 
-// In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
-// seeing which of the overloaded versions of
-// do_transform_real_to_unit_cell_internal() to call. This leads to bad
-// error messages and, generally, nothing very good. Avoid this by ensuring
-// that this class exists, but does not have an inner InternalData
-// type, thereby ruling out the codim-1 version of the function
-// below when doing overload resolution.
-template <>
-class MappingQGeneric<3,4>
-{};
 
-namespace
-{
-  /**
-   * Using the relative weights of the shape functions evaluated at
-   * one point on the reference cell (and stored in data.shape_values
-   * and accessed via data.shape(0,i)) and the locations of mapping
-   * support points (stored in data.mapping_support_points), compute
-   * the mapped location of that point in real space.
-   */
-  template <int dim, int spacedim>
-  Point<spacedim>
-  compute_mapped_location_of_point (const typename MappingQGeneric<dim,spacedim>::InternalData &data)
-  {
-    AssertDimension (data.shape_values.size(),
-                     data.mapping_support_points.size());
+        // Newton iteration to solve
+        //    f(x)=p(x)-p=0
+        // where we are looking for 'x' and p(x) is the forward transformation
+        // from unit to real cell. We solve this using a Newton iteration
+        //    x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+        // The start value is set to be the linear approximation to the cell
 
-    // use now the InternalData to compute the point in real space.
-    Point<spacedim> p_real;
-    for (unsigned int i=0; i<data.mapping_support_points.size(); ++i)
-      p_real += data.mapping_support_points[i] * data.shape(0,i);
+        // The shape values and derivatives of the mapping at this point are
+        // previously computed.
 
-    return p_real;
-  }
+        Point<dim> p_unit = initial_p_unit;
 
+        mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
-  /**
-   * Implementation of transform_real_to_unit_cell for dim==spacedim
-   */
-  template <int dim>
-  Point<dim>
-  do_transform_real_to_unit_cell_internal
-  (const typename Triangulation<dim,dim>::cell_iterator &cell,
-   const Point<dim>                                     &p,
-   const Point<dim>                                     &initial_p_unit,
-   typename MappingQGeneric<dim,dim>::InternalData      &mdata)
-  {
-    const unsigned int spacedim = dim;
-
-    const unsigned int n_shapes=mdata.shape_values.size();
-    (void)n_shapes;
-    Assert(n_shapes!=0, ExcInternalError());
-    AssertDimension (mdata.shape_derivatives.size(), n_shapes);
-
-    std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
-    AssertDimension (points.size(), n_shapes);
-
-
-    // Newton iteration to solve
-    //    f(x)=p(x)-p=0
-    // where we are looking for 'x' and p(x) is the forward transformation
-    // from unit to real cell. We solve this using a Newton iteration
-    //    x_{n+1}=x_n-[f'(x)]^{-1}f(x)
-    // The start value is set to be the linear approximation to the cell
-
-    // The shape values and derivatives of the mapping at this point are
-    // previously computed.
-
-    Point<dim> p_unit = initial_p_unit;
-
-    mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
-
-    Point<spacedim> p_real = compute_mapped_location_of_point<dim,spacedim>(mdata);
-    Tensor<1,spacedim> f = p_real-p;
-
-    // early out if we already have our point
-    if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter())
-      return p_unit;
-
-    // we need to compare the position of the computed p(x) against the given
-    // point 'p'. We will terminate the iteration and return 'x' if they are
-    // less than eps apart. The question is how to choose eps -- or, put maybe
-    // more generally: in which norm we want these 'p' and 'p(x)' to be eps
-    // apart.
-    //
-    // the question is difficult since we may have to deal with very elongated
-    // cells where we may achieve 1e-12*h for the distance of these two points
-    // in the 'long' direction, but achieving this tolerance in the 'short'
-    // direction of the cell may not be possible
-    //
-    // what we do instead is then to terminate iterations if
-    //    \| p(x) - p \|_A < eps
-    // where the A-norm is somehow induced by the transformation of the cell.
-    // in particular, we want to measure distances relative to the sizes of
-    // the cell in its principal directions.
-    //
-    // to define what exactly A should be, note that to first order we have
-    // the following (assuming that x* is the solution of the problem, i.e.,
-    // p(x*)=p):
-    //    p(x) - p = p(x) - p(x*)
-    //             = -grad p(x) * (x*-x) + higher order terms
-    // This suggest to measure with a norm that corresponds to
-    //    A = {[grad p(x]^T [grad p(x)]}^{-1}
-    // because then
-    //    \| p(x) - p \|_A  \approx  \| x - x* \|
-    // Consequently, we will try to enforce that
-    //    \| p(x) - p \|_A  =  \| f \|  <=  eps
-    //
-    // Note that using this norm is a bit dangerous since the norm changes
-    // in every iteration (A isn't fixed by depends on xk). However, if the
-    // cell is not too deformed (it may be stretched, but not twisted) then
-    // the mapping is almost linear and A is indeed constant or nearly so.
-    const double eps = 1.e-11;
-    const unsigned int newton_iteration_limit = 20;
-
-    unsigned int newton_iteration = 0;
-    double last_f_weighted_norm;
-    do
-      {
+        Point<spacedim> p_real = compute_mapped_location_of_point<dim,spacedim>(mdata);
+        Tensor<1,spacedim> f = p_real-p;
+
+        // early out if we already have our point
+        if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter())
+          return p_unit;
+
+        // we need to compare the position of the computed p(x) against the given
+        // point 'p'. We will terminate the iteration and return 'x' if they are
+        // less than eps apart. The question is how to choose eps -- or, put maybe
+        // more generally: in which norm we want these 'p' and 'p(x)' to be eps
+        // apart.
+        //
+        // the question is difficult since we may have to deal with very elongated
+        // cells where we may achieve 1e-12*h for the distance of these two points
+        // in the 'long' direction, but achieving this tolerance in the 'short'
+        // direction of the cell may not be possible
+        //
+        // what we do instead is then to terminate iterations if
+        //    \| p(x) - p \|_A < eps
+        // where the A-norm is somehow induced by the transformation of the cell.
+        // in particular, we want to measure distances relative to the sizes of
+        // the cell in its principal directions.
+        //
+        // to define what exactly A should be, note that to first order we have
+        // the following (assuming that x* is the solution of the problem, i.e.,
+        // p(x*)=p):
+        //    p(x) - p = p(x) - p(x*)
+        //             = -grad p(x) * (x*-x) + higher order terms
+        // This suggest to measure with a norm that corresponds to
+        //    A = {[grad p(x]^T [grad p(x)]}^{-1}
+        // because then
+        //    \| p(x) - p \|_A  \approx  \| x - x* \|
+        // Consequently, we will try to enforce that
+        //    \| p(x) - p \|_A  =  \| f \|  <=  eps
+        //
+        // Note that using this norm is a bit dangerous since the norm changes
+        // in every iteration (A isn't fixed by depends on xk). However, if the
+        // cell is not too deformed (it may be stretched, but not twisted) then
+        // the mapping is almost linear and A is indeed constant or nearly so.
+        const double eps = 1.e-11;
+        const unsigned int newton_iteration_limit = 20;
+
+        unsigned int newton_iteration = 0;
+        double last_f_weighted_norm;
+        do
+          {
 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
-        std::cout << "Newton iteration " << newton_iteration << std::endl;
+            std::cout << "Newton iteration " << newton_iteration << std::endl;
 #endif
 
-        // f'(x)
-        Tensor<2,spacedim> df;
-        for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
-          {
-            const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
-            const Point<spacedim> &point=points[k];
+            // f'(x)
+            Tensor<2,spacedim> df;
+            for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
+              {
+                const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
+                const Point<spacedim> &point=points[k];
 
-            for (unsigned int i=0; i<spacedim; ++i)
-              for (unsigned int j=0; j<dim; ++j)
-                df[i][j]+=point[i]*grad_transform[j];
-          }
+                for (unsigned int i=0; i<spacedim; ++i)
+                  for (unsigned int j=0; j<dim; ++j)
+                    df[i][j]+=point[i]*grad_transform[j];
+              }
 
-        // Solve  [f'(x)]d=f(x)
-        AssertThrow(determinant(df) > 0,
-                    (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
-        Tensor<2,spacedim> df_inverse = invert(df);
-        const Tensor<1,spacedim> delta = df_inverse * static_cast<const Tensor<1,spacedim>&>(f);
+            // Solve  [f'(x)]d=f(x)
+            AssertThrow(determinant(df) > 0,
+                        (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+            Tensor<2,spacedim> df_inverse = invert(df);
+            const Tensor<1,spacedim> delta = df_inverse * static_cast<const Tensor<1,spacedim>&>(f);
 
 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
-        std::cout << "   delta=" << delta  << std::endl;
+            std::cout << "   delta=" << delta  << std::endl;
 #endif
 
-        // do a line search
-        double step_length = 1;
-        do
-          {
-            // update of p_unit. The spacedim-th component of transformed point
-            // is simply ignored in codimension one case. When this component is
-            // not zero, then we are projecting the point to the surface or
-            // curve identified by the cell.
-            Point<dim> p_unit_trial = p_unit;
-            for (unsigned int i=0; i<dim; ++i)
-              p_unit_trial[i] -= step_length * delta[i];
-
-            // shape values and derivatives
-            // at new p_unit point
-            mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit_trial));
-
-            // f(x)
-            Point<spacedim> p_real_trial = compute_mapped_location_of_point<dim,spacedim>(mdata);
-            const Tensor<1,spacedim> f_trial = p_real_trial-p;
+            // do a line search
+            double step_length = 1;
+            do
+              {
+                // update of p_unit. The spacedim-th component of transformed point
+                // is simply ignored in codimension one case. When this component is
+                // not zero, then we are projecting the point to the surface or
+                // curve identified by the cell.
+                Point<dim> p_unit_trial = p_unit;
+                for (unsigned int i=0; i<dim; ++i)
+                  p_unit_trial[i] -= step_length * delta[i];
+
+                // shape values and derivatives
+                // at new p_unit point
+                mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit_trial));
+
+                // f(x)
+                Point<spacedim> p_real_trial = internal::MappingQGeneric::compute_mapped_location_of_point<dim,spacedim>(mdata);
+                const Tensor<1,spacedim> f_trial = p_real_trial-p;
 
 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
-            std::cout << "     step_length=" << step_length << std::endl
-                      << "       ||f ||   =" << f.norm() << std::endl
-                      << "       ||f*||   =" << f_trial.norm() << std::endl
-                      << "       ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
+                std::cout << "     step_length=" << step_length << std::endl
+                          << "       ||f ||   =" << f.norm() << std::endl
+                          << "       ||f*||   =" << f_trial.norm() << std::endl
+                          << "       ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
 #endif
 
-            // see if we are making progress with the current step length
-            // and if not, reduce it by a factor of two and try again
-            //
-            // strictly speaking, we should probably use the same norm as we use
-            // for the outer algorithm. in practice, line search is just a
-            // crutch to find a "reasonable" step length, and so using the l2
-            // norm is probably just fine
-            if (f_trial.norm() < f.norm())
-              {
-                p_real = p_real_trial;
-                p_unit = p_unit_trial;
-                f = f_trial;
-                break;
+                // see if we are making progress with the current step length
+                // and if not, reduce it by a factor of two and try again
+                //
+                // strictly speaking, we should probably use the same norm as we use
+                // for the outer algorithm. in practice, line search is just a
+                // crutch to find a "reasonable" step length, and so using the l2
+                // norm is probably just fine
+                if (f_trial.norm() < f.norm())
+                  {
+                    p_real = p_real_trial;
+                    p_unit = p_unit_trial;
+                    f = f_trial;
+                    break;
+                  }
+                else if (step_length > 0.05)
+                  step_length /= 2;
+                else
+                  AssertThrow (false,
+                               (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
               }
-            else if (step_length > 0.05)
-              step_length /= 2;
-            else
+            while (true);
+
+            ++newton_iteration;
+            if (newton_iteration > newton_iteration_limit)
               AssertThrow (false,
                            (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+            last_f_weighted_norm = (df_inverse * f).norm();
           }
-        while (true);
-
-        ++newton_iteration;
-        if (newton_iteration > newton_iteration_limit)
-          AssertThrow (false,
-                       (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
-        last_f_weighted_norm = (df_inverse * f).norm();
-      }
-    while (last_f_weighted_norm > eps);
-
-    return p_unit;
-  }
-
-
-
-  /**
-   * Implementation of transform_real_to_unit_cell for dim==spacedim-1
-   */
-  template <int dim>
-  Point<dim>
-  do_transform_real_to_unit_cell_internal_codim1
-  (const typename Triangulation<dim,dim+1>::cell_iterator &cell,
-   const Point<dim+1>                                       &p,
-   const Point<dim>                                         &initial_p_unit,
-   typename MappingQGeneric<dim,dim+1>::InternalData       &mdata)
-  {
-    const unsigned int spacedim = dim+1;
-
-    const unsigned int n_shapes=mdata.shape_values.size();
-    (void)n_shapes;
-    Assert(n_shapes!=0, ExcInternalError());
-    Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
-    Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError());
+        while (last_f_weighted_norm > eps);
 
-    std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
-    Assert(points.size()==n_shapes, ExcInternalError());
-
-    Point<spacedim> p_minus_F;
-
-    Tensor<1,spacedim>  DF[dim];
-    Tensor<1,spacedim>  D2F[dim][dim];
-
-    Point<dim> p_unit = initial_p_unit;
-    Point<dim> f;
-    Tensor<2,dim>  df;
-
-    // Evaluate first and second derivatives
-    mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
-
-    for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
-      {
-        const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
-        const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
-        const Point<spacedim> &point_k = points[k];
-
-        for (unsigned int j=0; j<dim; ++j)
-          {
-            DF[j] += grad_phi_k[j] * point_k;
-            for (unsigned int l=0; l<dim; ++l)
-              D2F[j][l] += hessian_k[j][l] * point_k;
-          }
+        return p_unit;
       }
 
-    p_minus_F = p;
-    p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
-
 
-    for (unsigned int j=0; j<dim; ++j)
-      f[j] = DF[j] * p_minus_F;
 
-    for (unsigned int j=0; j<dim; ++j)
+      /**
+       * Implementation of transform_real_to_unit_cell for dim==spacedim-1
+       */
+      template <int dim>
+      Point<dim>
+      do_transform_real_to_unit_cell_internal_codim1
+      (const typename dealii::Triangulation<dim,dim+1>::cell_iterator &cell,
+       const Point<dim+1>                                             &p,
+       const Point<dim>                                               &initial_p_unit,
+       typename dealii::MappingQGeneric<dim,dim+1>::InternalData      &mdata)
       {
-        f[j] = DF[j] * p_minus_F;
-        for (unsigned int l=0; l<dim; ++l)
-          df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
-      }
+        const unsigned int spacedim = dim+1;
 
+        const unsigned int n_shapes=mdata.shape_values.size();
+        (void)n_shapes;
+        Assert(n_shapes!=0, ExcInternalError());
+        Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
+        Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError());
 
-    const double eps = 1.e-12*cell->diameter();
-    const unsigned int loop_limit = 10;
+        std::vector<Point<spacedim> > &points=mdata.mapping_support_points;
+        Assert(points.size()==n_shapes, ExcInternalError());
 
-    unsigned int loop=0;
+        Point<spacedim> p_minus_F;
 
-    while (f.norm()>eps && loop++<loop_limit)
-      {
-        // Solve  [df(x)]d=f(x)
-        const Tensor<1,dim> d = invert(df) * static_cast<const Tensor<1,dim>&>(f);
-        p_unit -= d;
+        Tensor<1,spacedim>  DF[dim];
+        Tensor<1,spacedim>  D2F[dim][dim];
 
-        for (unsigned int j=0; j<dim; ++j)
-          {
-            DF[j].clear();
-            for (unsigned int l=0; l<dim; ++l)
-              D2F[j][l].clear();
-          }
+        Point<dim> p_unit = initial_p_unit;
+        Point<dim> f;
+        Tensor<2,dim>  df;
 
+        // Evaluate first and second derivatives
         mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
         for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
@@ -1473,11 +1344,13 @@ namespace
               }
           }
 
-        //TODO: implement a line search here in much the same way as for
-        // the corresponding function above that does so for dim==spacedim
         p_minus_F = p;
         p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
 
+
+        for (unsigned int j=0; j<dim; ++j)
+          f[j] = DF[j] * p_minus_F;
+
         for (unsigned int j=0; j<dim; ++j)
           {
             f[j] = DF[j] * p_minus_F;
@@ -1485,66 +1358,747 @@ namespace
               df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
           }
 
-      }
 
+        const double eps = 1.e-12*cell->diameter();
+        const unsigned int loop_limit = 10;
 
-    // Here we check that in the last execution of while the first
-    // condition was already wrong, meaning the residual was below
-    // eps. Only if the first condition failed, loop will have been
-    // increased and tested, and thus have reached the limit.
-    AssertThrow (loop<loop_limit, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+        unsigned int loop=0;
 
-    return p_unit;
-  }
+        while (f.norm()>eps && loop++<loop_limit)
+          {
+            // Solve  [df(x)]d=f(x)
+            const Tensor<1,dim> d = invert(df) * static_cast<const Tensor<1,dim>&>(f);
+            p_unit -= d;
 
+            for (unsigned int j=0; j<dim; ++j)
+              {
+                DF[j].clear();
+                for (unsigned int l=0; l<dim; ++l)
+                  D2F[j][l].clear();
+              }
 
-}
+            mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
+            for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
+              {
+                const Tensor<1,dim>   &grad_phi_k = mdata.derivative(0,k);
+                const Tensor<2,dim>   &hessian_k  = mdata.second_derivative(0,k);
+                const Point<spacedim> &point_k = points[k];
 
+                for (unsigned int j=0; j<dim; ++j)
+                  {
+                    DF[j] += grad_phi_k[j] * point_k;
+                    for (unsigned int l=0; l<dim; ++l)
+                      D2F[j][l] += hessian_k[j][l] * point_k;
+                  }
+              }
 
-// visual studio freaks out when trying to determine if
-// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
-// candidate. So instead of letting the compiler pick the correct overload, we
-// use template specialization to make sure we pick up the right function to
-// call:
+            //TODO: implement a line search here in much the same way as for
+            // the corresponding function above that does so for dim==spacedim
+            p_minus_F = p;
+            p_minus_F -= compute_mapped_location_of_point<dim,spacedim>(mdata);
 
-template <int dim, int spacedim>
-Point<dim>
-MappingQGeneric<dim,spacedim>::
-transform_real_to_unit_cell_internal
-(const typename Triangulation<dim,spacedim>::cell_iterator &,
- const Point<spacedim> &,
- const Point<dim> &) const
-{
-  // default implementation (should never be called)
-  Assert(false, ExcInternalError());
-  return Point<dim>();
-}
+            for (unsigned int j=0; j<dim; ++j)
+              {
+                f[j] = DF[j] * p_minus_F;
+                for (unsigned int l=0; l<dim; ++l)
+                  df[j][l] = -DF[j]*DF[l] + D2F[j][l] * p_minus_F;
+              }
 
-template <>
-Point<1>
-MappingQGeneric<1,1>::
-transform_real_to_unit_cell_internal
-(const Triangulation<1,1>::cell_iterator &cell,
- const Point<1>                            &p,
- const Point<1>                                 &initial_p_unit) const
-{
-  const int dim = 1;
-  const int spacedim = 1;
+          }
 
-  const Quadrature<dim> point_quadrature(initial_p_unit);
 
-  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
-  if (spacedim>dim)
-    update_flags |= update_jacobian_grads;
-  std::unique_ptr<InternalData> mdata (get_data(update_flags,
-                                                point_quadrature));
+        // Here we check that in the last execution of while the first
+        // condition was already wrong, meaning the residual was below
+        // eps. Only if the first condition failed, loop will have been
+        // increased and tested, and thus have reached the limit.
+        AssertThrow (loop<loop_limit, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
 
-  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
+        return p_unit;
+      }
+
+      /**
+       * 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::MappingQGeneric<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::MappingQGeneric<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[0];
+
+              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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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];
+                              }
+                  }
+              }
+          }
+      }
+    }
+  }
+}
+
+
+
+template <int dim, int spacedim>
+MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
+  :
+  polynomial_degree(p),
+  line_support_points(this->polynomial_degree+1),
+  fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : nullptr),
+  support_point_weights_perimeter_to_interior (internal::MappingQGeneric::compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)),
+  support_point_weights_cell (internal::MappingQGeneric::compute_support_point_weights_cell<dim>(this->polynomial_degree))
+{
+  Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings "
+                              "with a polynomial degree greater or equal to one."));
+}
+
+
+
+template <int dim, int spacedim>
+MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
+  :
+  polynomial_degree(mapping.polynomial_degree),
+  line_support_points(mapping.line_support_points),
+  fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : nullptr),
+  support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior),
+  support_point_weights_cell (mapping.support_point_weights_cell)
+{}
+
+
+
+
+template <int dim, int spacedim>
+Mapping<dim,spacedim> *
+MappingQGeneric<dim,spacedim>::clone () const
+{
+  return new MappingQGeneric<dim,spacedim>(*this);
+}
+
+
+
+
+template <int dim, int spacedim>
+unsigned int
+MappingQGeneric<dim,spacedim>::get_degree() const
+{
+  return polynomial_degree;
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+MappingQGeneric<dim,spacedim>::
+transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                             const Point<dim> &p) const
+{
+  // set up the polynomial space
+  const TensorProductPolynomials<dim>
+  tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
+  Assert (tensor_pols.n() == Utilities::fixed_power<dim>(polynomial_degree+1),
+          ExcInternalError());
+
+  // then also construct the mapping from lexicographic to the Qp shape function numbering
+  const std::vector<unsigned int>
+  renumber (FETools::
+            lexicographic_to_hierarchic_numbering
+            (FiniteElementData<dim> (internal::MappingQGeneric::get_dpo_vector<dim>
+                                     (polynomial_degree), 1, polynomial_degree)));
+
+  const std::vector<Point<spacedim> > support_points
+    = this->compute_mapping_support_points(cell);
+
+  Point<spacedim> mapped_point;
+  for (unsigned int i=0; i<tensor_pols.n(); ++i)
+    mapped_point += support_points[renumber[i]] * tensor_pols.compute_value (i, p);
+
+  return mapped_point;
+}
+
+
+// In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
+// seeing which of the overloaded versions of
+// do_transform_real_to_unit_cell_internal() to call. This leads to bad
+// error messages and, generally, nothing very good. Avoid this by ensuring
+// that this class exists, but does not have an inner InternalData
+// type, thereby ruling out the codim-1 version of the function
+// below when doing overload resolution.
+template <>
+class MappingQGeneric<3,4>
+{};
+
+
+
+// visual studio freaks out when trying to determine if
+// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
+// candidate. So instead of letting the compiler pick the correct overload, we
+// use template specialization to make sure we pick up the right function to
+// call:
+
+template <int dim, int spacedim>
+Point<dim>
+MappingQGeneric<dim,spacedim>::
+transform_real_to_unit_cell_internal
+(const typename Triangulation<dim,spacedim>::cell_iterator &,
+ const Point<spacedim> &,
+ const Point<dim> &) const
+{
+  // default implementation (should never be called)
+  Assert(false, ExcInternalError());
+  return Point<dim>();
+}
+
+template <>
+Point<1>
+MappingQGeneric<1,1>::
+transform_real_to_unit_cell_internal
+(const Triangulation<1,1>::cell_iterator &cell,
+ const Point<1>                            &p,
+ const Point<1>                                 &initial_p_unit) const
+{
+  const int dim = 1;
+  const int spacedim = 1;
+
+  const Quadrature<dim> point_quadrature(initial_p_unit);
+
+  UpdateFlags update_flags = update_quadrature_points | update_jacobians;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std::unique_ptr<InternalData> mdata (get_data(update_flags,
+                                                point_quadrature));
+
+  mdata->mapping_support_points = this->compute_mapping_support_points (cell);
 
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
+  return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
 }
 
 template <>
@@ -1570,7 +2124,7 @@ transform_real_to_unit_cell_internal
 
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
+  return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
 }
 
 template <>
@@ -1596,7 +2150,7 @@ transform_real_to_unit_cell_internal
 
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
+  return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
 }
 
 template <>
@@ -1622,7 +2176,7 @@ transform_real_to_unit_cell_internal
 
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata);
+  return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata);
 }
 
 template <>
@@ -1648,7 +2202,7 @@ transform_real_to_unit_cell_internal
 
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
-  return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata);
+  return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata);
 }
 
 template <>
@@ -1866,584 +2420,39 @@ MappingQGeneric<dim,spacedim>::get_data (const UpdateFlags update_flags,
   data->initialize (this->requires_update_flags(update_flags), q, q.size());
 
   return data;
-}
-
-
-
-template <int dim, int spacedim>
-typename MappingQGeneric<dim,spacedim>::InternalData *
-MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
-                                              const Quadrature<dim-1> &quadrature) const
-{
-  InternalData *data = new InternalData(polynomial_degree);
-  data->initialize_face (this->requires_update_flags(update_flags),
-                         QProjector<dim>::project_to_all_faces(quadrature),
-                         quadrature.size());
-
-  return data;
-}
-
-
-
-template <int dim, int spacedim>
-typename MappingQGeneric<dim,spacedim>::InternalData *
-MappingQGeneric<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
-                                                 const Quadrature<dim-1>& quadrature) const
-{
-  InternalData *data = new InternalData(polynomial_degree);
-  data->initialize_face (this->requires_update_flags(update_flags),
-                         QProjector<dim>::project_to_all_subfaces(quadrature),
-                         quadrature.size());
-
-  return data;
-}
-
-
-
-namespace internal
-{
-  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::MappingQGeneric<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::MappingQGeneric<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[0];
-
-            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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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];
-                            }
+template <int dim, int spacedim>
+typename MappingQGeneric<dim,spacedim>::InternalData *
+MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
+                                              const Quadrature<dim-1> &quadrature) const
+{
+  InternalData *data = new InternalData(polynomial_degree);
+  data->initialize_face (this->requires_update_flags(update_flags),
+                         QProjector<dim>::project_to_all_faces(quadrature),
+                         quadrature.size());
 
-                  // 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];
-                            }
-                }
-            }
-        }
-    }
-  }
+  return data;
 }
 
 
 
+template <int dim, int spacedim>
+typename MappingQGeneric<dim,spacedim>::InternalData *
+MappingQGeneric<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
+                                                 const Quadrature<dim-1>& quadrature) const
+{
+  InternalData *data = new InternalData(polynomial_degree);
+  data->initialize_face (this->requires_update_flags(update_flags),
+                         QProjector<dim>::project_to_all_subfaces(quadrature),
+                         quadrature.size());
+
+  return data;
+}
+
+
 
 template <int dim, int spacedim>
 CellSimilarity::Similarity
@@ -2481,12 +2490,14 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const CellSimilarity::Similarity computed_cell_similarity =
     (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
 
-  internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
-                                                  data,
-                                                  output_data.quadrature_points);
-  internal::maybe_update_Jacobians<dim,spacedim> (computed_cell_similarity,
-                                                  QProjector<dim>::DataSetDescriptor::cell (),
-                                                  data);
+  internal::MappingQGeneric::maybe_compute_q_points<dim,spacedim>
+  (QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.quadrature_points);
+  internal::MappingQGeneric::maybe_update_Jacobians<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data);
 
   const UpdateFlags update_flags = data.update_each;
   const std::vector<double> &weights=quadrature.get_weights();
@@ -2599,35 +2610,41 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
           output_data.inverse_jacobians[point] = data.covariant[point].transpose();
     }
 
-  internal::maybe_update_jacobian_grads<dim,spacedim> (computed_cell_similarity,
-                                                       QProjector<dim>::DataSetDescriptor::cell (),
-                                                       data,
-                                                       output_data.jacobian_grads);
-
-  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_pushed_forward_grads);
-
-  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_2nd_derivatives);
-
-  internal::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::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_3rd_derivatives);
-
-  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell (),
-      data,
-      output_data.jacobian_pushed_forward_3rd_derivatives);
+  internal::MappingQGeneric::maybe_update_jacobian_grads<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.jacobian_grads);
+
+  internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_grads<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.jacobian_pushed_forward_grads);
+
+  internal::MappingQGeneric::maybe_update_jacobian_2nd_derivatives<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.jacobian_2nd_derivatives);
+
+  internal::MappingQGeneric::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::MappingQGeneric::maybe_update_jacobian_3rd_derivatives<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.jacobian_3rd_derivatives);
+
+  internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim>
+  (computed_cell_similarity,
+   QProjector<dim>::DataSetDescriptor::cell (),
+   data,
+   output_data.jacobian_pushed_forward_3rd_derivatives);
 
   return computed_cell_similarity;
 }
@@ -2639,214 +2656,217 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
 namespace internal
 {
-  namespace
+  namespace MappingQGeneric
   {
-    /**
-     * 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::MappingQGeneric<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 std::vector<double>        &weights,
-                             const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &data,
-                             internal::FEValues::MappingRelatedData<dim,spacedim>         &output_data)
+    namespace
     {
-      const UpdateFlags update_flags = data.update_each;
+      /**
+       * 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::MappingQGeneric<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 std::vector<double>        &weights,
+                               const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &data,
+                               internal::FEValues::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+GeometryInfo<dim>::faces_per_cell*d <
-                      data.unit_tangentials.size(),
-                      ExcInternalError());
-              Assert (data.aux[d].size() <=
-                      data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
-                      ExcInternalError());
-
-              mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d]),
-                                 mapping_contravariant,
-                                 data,
-                                 make_array_view(data.aux[d]));
-            }
+        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+GeometryInfo<dim>::faces_per_cell*d <
+                        data.unit_tangentials.size(),
+                        ExcInternalError());
+                Assert (data.aux[d].size() <=
+                        data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
+                        ExcInternalError());
+
+                mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*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)
+            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)
                         {
-                          // 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();
+                        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 (dim==2)
+                          {
+                            const DerivativeForm<1,spacedim,dim> DX_t =
+                              data.contravariant[point].transpose();
 
-          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() * weights[i];
+                            Tensor<1, spacedim> cell_normal =
+                              cross_product_3d(DX_t[0], DX_t[1]);
+                            cell_normal /= cell_normal.norm();
 
-                if (subface_no != numbers::invalid_unsigned_int)
-                  {
-                    const double area_ratio = GeometryInfo<dim>::subface_ratio(cell->subface_case(face_no),
-                                                                               subface_no);
-                    output_data.JxW_values[i] *= area_ratio;
+                            // then compute the face normal from the face tangent
+                            // and the cell normal:
+                            output_data.boundary_forms[point] =
+                              cross_product_3d(data.aux[0][point], cell_normal);
+                          }
+                      }
                   }
               }
 
-          if (update_flags & update_normal_vectors)
-            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_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() * weights[i];
+
+                  if (subface_no != numbers::invalid_unsigned_int)
+                    {
+                      const double area_ratio = GeometryInfo<dim>::subface_ratio(cell->subface_case(face_no),
+                                                                                 subface_no);
+                      output_data.JxW_values[i] *= area_ratio;
+                    }
+                }
+
+            if (update_flags & update_normal_vectors)
+              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_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();
-        }
-    }
+            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 MappingQGeneric::fill_fe_face_values() and
-     * MappingQGeneric::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::MappingQGeneric<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::MappingQGeneric<dim,spacedim>::InternalData      &data,
-                            internal::FEValues::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(),
-                               quadrature.get_weights(), data,
-                               output_data);
+      /**
+       * Do the work of MappingQGeneric::fill_fe_face_values() and
+       * MappingQGeneric::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::MappingQGeneric<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::MappingQGeneric<dim,spacedim>::InternalData      &data,
+                              internal::FEValues::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(),
+                                 quadrature.get_weights(), data,
+                                 output_data);
+      }
     }
   }
 }
@@ -2883,16 +2903,17 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
       data.cell_of_current_support_points = cell;
     }
 
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, numbers::invalid_unsigned_int,
-                                    QProjector<dim>::DataSetDescriptor::face (face_no,
-                                        cell->face_orientation(face_no),
-                                        cell->face_flip(face_no),
-                                        cell->face_rotation(face_no),
-                                        quadrature.size()),
-                                    quadrature,
-                                    data,
-                                    output_data);
+  internal::MappingQGeneric::do_fill_fe_face_values
+  (*this,
+   cell, face_no, numbers::invalid_unsigned_int,
+   QProjector<dim>::DataSetDescriptor::face (face_no,
+                                             cell->face_orientation(face_no),
+                                             cell->face_flip(face_no),
+                                             cell->face_rotation(face_no),
+                                             quadrature.size()),
+   quadrature,
+   data,
+   output_data);
 }
 
 
@@ -2928,337 +2949,363 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
       data.cell_of_current_support_points = cell;
     }
 
-  internal::do_fill_fe_face_values (*this,
-                                    cell, face_no, subface_no,
-                                    QProjector<dim>::DataSetDescriptor::subface (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);
+  internal::MappingQGeneric::do_fill_fe_face_values
+  (*this,
+   cell, face_no, subface_no,
+   QProjector<dim>::DataSetDescriptor::subface (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
+namespace internal
 {
-  template <int dim, int spacedim, int rank>
-  void
-  transform_fields(const ArrayView<const Tensor<rank,dim> >               &input,
-                   const MappingType                                       mapping_type,
-                   const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
-                   const ArrayView<Tensor<rank,spacedim> >                &output)
+  namespace MappingQGeneric
   {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingQGeneric<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
-
-    switch (mapping_type)
-      {
-      case mapping_contravariant:
+    namespace
+    {
+      // We cannot query a manifold from the faces of a 1D elements (i.e.,
+      // vertices), which is why we add a specialization for the 3D case here
+      template <typename Iterator>
+      bool check_identical_manifolds_of_quads(const Iterator &)
       {
-        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;
+        Assert(false, ExcNotImplemented());
+        return true;
       }
 
-      case mapping_piola:
+      bool check_identical_manifolds_of_quads(const dealii::Triangulation<3,3>::cell_iterator &cell)
       {
-        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;
+        for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+          if (&cell->face(f)->get_manifold() != &cell->get_manifold())
+            return false;
+        return true;
       }
-      //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_fields(const ArrayView<const Tensor<rank,dim> >               &input,
+                       const MappingType                                       mapping_type,
+                       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::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingQGeneric<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
+        switch (mapping_type)
+          {
+          case mapping_contravariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
 
-  template <int dim, int spacedim, int rank>
-  void
-  transform_gradients(const ArrayView<const Tensor<rank,dim> >                &input,
-                      const MappingType                                        mapping_type,
-                      const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
-                      const ArrayView<Tensor<rank,spacedim> >                 &output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingQGeneric<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.contravariant[i], input[i]);
 
-    switch (mapping_type)
-      {
-      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"));
+            return;
+          }
 
-        for (unsigned int i=0; i<output.size(); ++i)
+          case mapping_piola:
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.contravariant[i], transpose(input[i]) );
-            output[i] = apply_transformation(data.covariant[i], A.transpose() );
+            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;
           }
 
-        return;
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
 
-      case mapping_covariant_gradient:
+
+      template <int dim, int spacedim, int rank>
+      void
+      transform_gradients(const ArrayView<const Tensor<rank,dim> >                &input,
+                          const MappingType                                        mapping_type,
+                          const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
+                          const ArrayView<Tensor<rank,spacedim> >                 &output)
       {
-        Assert (data.update_each & update_covariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
-        Assert (rank==2, ExcMessage("Only for rank 2"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingQGeneric<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
+        switch (mapping_type)
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.covariant[i], transpose(input[i]) );
-            output[i] = apply_transformation(data.covariant[i], A.transpose() );
-          }
+          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"));
 
-        return;
-      }
+            for (unsigned int i=0; i<output.size(); ++i)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.contravariant[i], transpose(input[i]) );
+                output[i] = apply_transformation(data.covariant[i], A.transpose() );
+              }
 
-      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)
+            return;
+          }
+
+          case mapping_covariant_gradient:
           {
-            DerivativeForm<1,spacedim,dim> A =
-              apply_transformation(data.covariant[i], input[i] );
-            Tensor<2,spacedim> T =
-              apply_transformation(data.contravariant[i], A.transpose() );
+            Assert (data.update_each & update_covariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+            Assert (rank==2, ExcMessage("Only for rank 2"));
 
-            output[i] = transpose(T);
-            output[i] /= data.volume_elements[i];
+            for (unsigned int i=0; i<output.size(); ++i)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.covariant[i], transpose(input[i]) );
+                output[i] = apply_transformation(data.covariant[i], A.transpose() );
+              }
+
+            return;
           }
 
-        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)
+              {
+                DerivativeForm<1,spacedim,dim> A =
+                  apply_transformation(data.covariant[i], input[i] );
+                Tensor<2,spacedim> T =
+                  apply_transformation(data.contravariant[i], A.transpose() );
 
-      default:
-        Assert(false, ExcNotImplemented());
-      }
-  }
+                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 MappingType                                       mapping_type,
-                     const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
-                     const ArrayView<Tensor<3,spacedim> >                   &output)
-  {
-    AssertDimension (input.size(), output.size());
-    Assert ((dynamic_cast<const typename MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingQGeneric<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
-    switch (mapping_type)
-      {
-      case mapping_contravariant_hessian:
+      template <int dim, int spacedim>
+      void
+      transform_hessians(const ArrayView<const Tensor<3,dim> >                  &input,
+                         const MappingType                                       mapping_type,
+                         const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                         const ArrayView<Tensor<3,spacedim> >                   &output)
       {
-        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"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingQGeneric<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
-        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)
+        switch (mapping_type)
+          {
+          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 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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            return;
+          }
 
-      case mapping_covariant_hessian:
-      {
-        Assert (data.update_each & update_covariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+          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)
+            for (unsigned int q=0; q<output.size(); ++q)
+              for (unsigned int i=0; i<spacedim; ++i)
                 {
-                  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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            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)
+          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 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)
+                  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)
                     {
-                      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];
+                      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;
-      }
+            return;
+          }
 
-      default:
-        Assert(false, ExcNotImplemented());
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
-  }
-
 
 
 
-  template <int dim, int spacedim, int rank>
-  void
-  transform_differential_forms(const ArrayView<const DerivativeForm<rank, dim,spacedim> >   &input,
-                               const MappingType                                             mapping_type,
-                               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 MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
-            ExcInternalError());
-    const typename MappingQGeneric<dim,spacedim>::InternalData
-    &data = static_cast<const typename MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
-    switch (mapping_type)
-      {
-      case mapping_covariant:
+      template <int dim, int spacedim, int rank>
+      void
+      transform_differential_forms(const ArrayView<const DerivativeForm<rank, dim,spacedim> >   &input,
+                                   const MappingType                                             mapping_type,
+                                   const typename Mapping<dim,spacedim>::InternalDataBase       &mapping_data,
+                                   const ArrayView<Tensor<rank+1, spacedim> >                   &output)
       {
-        Assert (data.update_each & update_contravariant_transformation,
-                typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+        AssertDimension (input.size(), output.size());
+        Assert ((dynamic_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData *>(&mapping_data) != nullptr),
+                ExcInternalError());
+        const typename dealii::MappingQGeneric<dim,spacedim>::InternalData
+        &data = static_cast<const typename dealii::MappingQGeneric<dim,spacedim>::InternalData &>(mapping_data);
 
-        for (unsigned int i=0; i<output.size(); ++i)
-          output[i] = apply_transformation(data.covariant[i], input[i]);
+        switch (mapping_type)
+          {
+          case mapping_covariant:
+          {
+            Assert (data.update_each & update_contravariant_transformation,
+                    typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
 
-        return;
-      }
-      default:
-        Assert(false, ExcNotImplemented());
+            for (unsigned int i=0; i<output.size(); ++i)
+              output[i] = apply_transformation(data.covariant[i], input[i]);
+
+            return;
+          }
+          default:
+            Assert(false, ExcNotImplemented());
+          }
       }
+    }
   }
 }
 
@@ -3272,7 +3319,7 @@ transform (const ArrayView<const Tensor<1, dim> >                  &input,
            const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
            const ArrayView<Tensor<1, spacedim> >                   &output) const
 {
-  transform_fields(input, mapping_type, mapping_data, output);
+  internal::MappingQGeneric::transform_fields(input, mapping_type, mapping_data, output);
 }
 
 
@@ -3285,7 +3332,7 @@ transform (const ArrayView<const DerivativeForm<1, dim,spacedim> >  &input,
            const typename Mapping<dim,spacedim>::InternalDataBase   &mapping_data,
            const ArrayView<Tensor<2, spacedim> >                    &output) const
 {
-  transform_differential_forms(input, mapping_type, mapping_data, output);
+  internal::MappingQGeneric::transform_differential_forms(input, mapping_type, mapping_data, output);
 }
 
 
@@ -3301,13 +3348,13 @@ transform (const ArrayView<const Tensor<2, dim> >                  &input,
   switch (mapping_type)
     {
     case mapping_contravariant:
-      transform_fields(input, mapping_type, mapping_data, output);
+      internal::MappingQGeneric::transform_fields(input, mapping_type, mapping_data, output);
       return;
 
     case mapping_piola_gradient:
     case mapping_contravariant_gradient:
     case mapping_covariant_gradient:
-      transform_gradients(input, mapping_type, mapping_data, output);
+      internal::MappingQGeneric::transform_gradients(input, mapping_type, mapping_data, output);
       return;
     default:
       Assert(false, ExcNotImplemented());
@@ -3378,7 +3425,7 @@ transform (const ArrayView<const  Tensor<3,dim> >                  &input,
     case mapping_piola_hessian:
     case mapping_contravariant_hessian:
     case mapping_covariant_hessian:
-      transform_hessians(input, mapping_type, mapping_data, output);
+      internal::MappingQGeneric::transform_hessians(input, mapping_type, mapping_data, output);
       return;
     default:
       Assert(false, ExcNotImplemented());
@@ -3387,28 +3434,6 @@ transform (const ArrayView<const  Tensor<3,dim> >                  &input,
 
 
 
-namespace
-{
-  // We cannot query a manifold from the faces of a 1D elements (i.e.,
-  // vertices), which is why we add a specialization for the 3D case here
-  template <typename Iterator>
-  bool check_identical_manifolds_of_quads(const Iterator &)
-  {
-    Assert(false, ExcNotImplemented());
-    return true;
-  }
-
-  bool check_identical_manifolds_of_quads(const Triangulation<3,3>::cell_iterator &cell)
-  {
-    for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
-      if (&cell->face(f)->get_manifold() != &cell->get_manifold())
-        return false;
-    return true;
-  }
-}
-
-
-
 template <int dim, int spacedim>
 void
 MappingQGeneric<dim,spacedim>::
@@ -3642,7 +3667,7 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
               if (&cell->line(l)->get_manifold() != &cell->get_manifold())
                 all_manifold_ids_are_equal = false;
           if (dim == 3)
-            if (check_identical_manifolds_of_quads(cell) == false)
+            if (internal::MappingQGeneric::check_identical_manifolds_of_quads(cell) == false)
               all_manifold_ids_are_equal = false;
         }
 

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.