]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid the use of anonymous namespaces in a header file. 18505/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 26 May 2025 23:08:11 +0000 (17:08 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 26 May 2025 23:08:43 +0000 (17:08 -0600)
include/deal.II/fe/mapping_fe_field.templates.h

index ea7b935647d8b10732f4d9b87f46fe205c9a1bde..c132cced577bf490eef6edbf6c29cad20cac0931 100644 (file)
@@ -687,819 +687,797 @@ namespace internal
 {
   namespace MappingFEFieldImplementation
   {
-    namespace
+    /**
+     * Compute the locations of quadrature points on the object described by
+     * the first argument (and the cell for which the mapping support points
+     * have already been set), but only if the update_flags of the @p data
+     * argument indicate so.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_compute_q_points(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real,
+      std::vector<Point<spacedim>>       &quadrature_points)
     {
-      /**
-       * Compute the locations of quadrature points on the object described by
-       * the first argument (and the cell for which the mapping support points
-       * have already been set), but only if the update_flags of the @p data
-       * argument indicate so.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_compute_q_points(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Point<spacedim>>       &quadrature_points)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        if (update_flags & update_quadrature_points)
-          {
-            for (unsigned int point = 0; point < quadrature_points.size();
-                 ++point)
-              {
-                Point<spacedim> result;
-                const double   *shape = &data.shape(point + data_set, 0);
+      const UpdateFlags update_flags = data.update_each;
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      result[fe_to_real[comp_k]] +=
-                        data.local_dof_values[k] * shape[k];
-                  }
+      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);
 
-                quadrature_points[point] = result;
-              }
-          }
-      }
-
-      /**
-       * Update the co- and contravariant matrices as well as their determinant,
-       * for the cell described stored in the data object, but only if the
-       * update_flags of the @p data argument indicate so.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_Jacobians(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        // then Jacobians
-        if (update_flags & update_contravariant_transformation)
-          {
-            const unsigned int n_q_points = data.contravariant.size();
-
-            Assert(data.n_shape_functions > 0, ExcInternalError());
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<1, dim> *data_derv =
-                  &data.derivative(point + data_set, 0);
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    result[fe_to_real[comp_k]] +=
+                      data.local_dof_values[k] * shape[k];
+                }
 
-                Tensor<1, dim> result[spacedim];
+              quadrature_points[point] = result;
+            }
+        }
+    }
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      result[fe_to_real[comp_k]] +=
-                        data.local_dof_values[k] * data_derv[k];
-                  }
+    /**
+     * Update the co- and contravariant matrices as well as their determinant,
+     * for the cell described stored in the data object, but only if the
+     * update_flags of the @p data argument indicate so.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_Jacobians(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real)
+    {
+      const UpdateFlags update_flags = data.update_each;
 
-                // write result into contravariant data
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  {
-                    data.contravariant[point][i] = result[i];
-                  }
-              }
-          }
-
-        if (update_flags & update_covariant_transformation)
-          {
-            AssertDimension(data.covariant.size(), data.contravariant.size());
-            for (unsigned int point = 0; point < data.contravariant.size();
-                 ++point)
-              data.covariant[point] =
-                (data.contravariant[point]).covariant_form();
-          }
-
-        if (update_flags & update_volume_elements)
-          {
-            AssertDimension(data.contravariant.size(),
-                            data.volume_elements.size());
-            for (unsigned int point = 0; point < data.contravariant.size();
-                 ++point)
-              data.volume_elements[point] =
-                data.contravariant[point].determinant();
-          }
-      }
-
-      /**
-       * Update the Hessian of the transformation from unit to real cell, the
-       * Jacobian gradients.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_grads(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_grads)
-          {
-            const unsigned int n_q_points = jacobian_grads.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<2, dim> *second =
-                  &data.second_derivative(point + data_set, 0);
+      // then Jacobians
+      if (update_flags & update_contravariant_transformation)
+        {
+          const unsigned int n_q_points = data.contravariant.size();
 
-                DerivativeForm<2, dim, spacedim> result;
+          Assert(data.n_shape_functions > 0, ExcInternalError());
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          result[fe_to_real[comp_k]][j][l] +=
-                            (second[k][j][l] * data.local_dof_values[k]);
-                  }
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<1, dim> *data_derv =
+                &data.derivative(point + data_set, 0);
 
-                // never touch any data for j=dim in case dim<spacedim, so
-                // it will always be zero as it was initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      jacobian_grads[point][i][j][l] = result[i][j][l];
-              }
-          }
-      }
-
-      /**
-       * Update the Hessian of the transformation from unit to real cell, the
-       * Jacobian gradients, pushed forward to the real cell coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_grads(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<3, spacedim>>   &jacobian_pushed_forward_grads)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_grads)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_grads.size();
-
-            double tmp[spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<2, dim> *second =
-                  &data.second_derivative(point + data_set, 0);
+              Tensor<1, dim> result[spacedim];
 
-                DerivativeForm<2, dim, spacedim> result;
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    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)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          result[fe_to_real[comp_k]][j][l] +=
-                            (second[k][j][l] * data.local_dof_values[k]);
-                  }
+              // write result into contravariant data
+              for (unsigned int i = 0; i < spacedim; ++i)
+                {
+                  data.contravariant[point][i] = result[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)
+      if (update_flags & update_covariant_transformation)
+        {
+          AssertDimension(data.covariant.size(), data.contravariant.size());
+          for (unsigned int point = 0; point < data.contravariant.size();
+               ++point)
+            data.covariant[point] =
+              (data.contravariant[point]).covariant_form();
+        }
+
+      if (update_flags & update_volume_elements)
+        {
+          AssertDimension(data.contravariant.size(),
+                          data.volume_elements.size());
+          for (unsigned int point = 0; point < data.contravariant.size();
+               ++point)
+            data.volume_elements[point] =
+              data.contravariant[point].determinant();
+        }
+    }
+
+    /**
+     * Update the Hessian of the transformation from unit to real cell, the
+     * Jacobian gradients.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_grads(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                                &data,
+      const FiniteElement<dim, spacedim>            &fe,
+      const ComponentMask                           &fe_mask,
+      const std::vector<unsigned int>               &fe_to_real,
+      std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_grads)
+        {
+          const unsigned int n_q_points = jacobian_grads.size();
+
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<2, dim> *second =
+                &data.second_derivative(point + data_set, 0);
+
+              DerivativeForm<2, dim, spacedim> result;
+
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        result[fe_to_real[comp_k]][j][l] +=
+                          (second[k][j][l] * data.local_dof_values[k]);
+                }
+
+              // never touch any data for j=dim in case dim<spacedim, so
+              // it will always be zero as it was initialized
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < dim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    jacobian_grads[point][i][j][l] = result[i][j][l];
+            }
+        }
+    }
+
+    /**
+     * Update the Hessian of the transformation from unit to real cell, the
+     * Jacobian gradients, pushed forward to the real cell coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_pushed_forward_grads(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real,
+      std::vector<Tensor<3, spacedim>>   &jacobian_pushed_forward_grads)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_pushed_forward_grads)
+        {
+          const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
+
+          double tmp[spacedim][spacedim][spacedim];
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<2, dim> *second =
+                &data.second_derivative(point + data_set, 0);
+
+              DerivativeForm<2, dim, spacedim> result;
+
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        result[fe_to_real[comp_k]][j][l] +=
+                          (second[k][j][l] * data.local_dof_values[k]);
+                }
+
+              // first push forward the j-components
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < spacedim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    {
+                      tmp[i][j][l] =
+                        result[i][0][l] * data.covariant[point][j][0];
+                      for (unsigned int jr = 1; jr < dim; ++jr)
+                        {
+                          tmp[i][j][l] +=
+                            result[i][jr][l] * data.covariant[point][j][jr];
+                        }
+                    }
+
+              // now, pushing forward the l-components
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < spacedim; ++j)
+                  for (unsigned int l = 0; l < spacedim; ++l)
+                    {
+                      jacobian_pushed_forward_grads[point][i][j][l] =
+                        tmp[i][j][0] * data.covariant[point][l][0];
+                      for (unsigned int lr = 1; lr < dim; ++lr)
+                        {
+                          jacobian_pushed_forward_grads[point][i][j][l] +=
+                            tmp[i][j][lr] * data.covariant[point][l][lr];
+                        }
+                    }
+            }
+        }
+    }
+
+    /**
+     * Update the third derivative of the transformation from unit to real
+     * cell, the Jacobian hessians.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_2nd_derivatives(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                                &data,
+      const FiniteElement<dim, spacedim>            &fe,
+      const ComponentMask                           &fe_mask,
+      const std::vector<unsigned int>               &fe_to_real,
+      std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_2nd_derivatives)
+        {
+          const unsigned int n_q_points = jacobian_2nd_derivatives.size();
+
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<3, dim> *third =
+                &data.third_derivative(point + data_set, 0);
+
+              DerivativeForm<3, dim, spacedim> result;
+
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        for (unsigned int m = 0; m < dim; ++m)
+                          result[fe_to_real[comp_k]][j][l][m] +=
+                            (third[k][j][l][m] * data.local_dof_values[k]);
+                }
+
+              // never touch any data for j=dim in case dim<spacedim, so
+              // it will always be zero as it was initialized
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < dim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    for (unsigned int m = 0; m < dim; ++m)
+                      jacobian_2nd_derivatives[point][i][j][l][m] =
+                        result[i][j][l][m];
+            }
+        }
+    }
+
+    /**
+     * Update the third derivative of the transformation from unit to real
+     * cell, the Jacobian hessians, pushed forward to the real cell
+     * coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_pushed_forward_2nd_derivatives(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real,
+      std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
+        {
+          const unsigned int n_q_points =
+            jacobian_pushed_forward_2nd_derivatives.size();
+
+          double tmp[spacedim][spacedim][spacedim][spacedim];
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<3, dim> *third =
+                &data.third_derivative(point + data_set, 0);
+
+              DerivativeForm<3, dim, spacedim> result;
+
+              for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        for (unsigned int m = 0; m < dim; ++m)
+                          result[fe_to_real[comp_k]][j][l][m] +=
+                            (third[k][j][l][m] * data.local_dof_values[k]);
+                }
+
+              // push forward the j-coordinate
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < spacedim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    for (unsigned int m = 0; m < dim; ++m)
                       {
-                        tmp[i][j][l] =
-                          result[i][0][l] * data.covariant[point][j][0];
+                        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)
-                          {
-                            tmp[i][j][l] +=
-                              result[i][jr][l] * data.covariant[point][j][jr];
-                          }
+                          jacobian_pushed_forward_2nd_derivatives[point][i][j]
+                                                                 [l][m] +=
+                            result[i][jr][l][m] * 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)
+              // 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)
                       {
-                        jacobian_pushed_forward_grads[point][i][j][l] =
-                          tmp[i][j][0] * data.covariant[point][l][0];
+                        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)
-                          {
-                            jacobian_pushed_forward_grads[point][i][j][l] +=
-                              tmp[i][j][lr] * data.covariant[point][l][lr];
-                          }
+                          tmp[i][j][l][m] +=
+                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
+                                                                   [lr][m] *
+                            data.covariant[point][l][lr];
                       }
-              }
-          }
-      }
-
-      /**
-       * Update the third derivative of the transformation from unit to real
-       * cell, the Jacobian hessians.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_2nd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_2nd_derivatives)
-          {
-            const unsigned int n_q_points = jacobian_2nd_derivatives.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<3, dim> *third =
-                  &data.third_derivative(point + data_set, 0);
 
-                DerivativeForm<3, dim, spacedim> result;
+              // 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];
+                      }
+            }
+        }
+    }
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            result[fe_to_real[comp_k]][j][l][m] +=
-                              (third[k][j][l][m] * data.local_dof_values[k]);
-                  }
+    /**
+     * Update the fourth derivative of the transformation from unit to real
+     * cell, the Jacobian hessian gradients.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_3rd_derivatives(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                                &data,
+      const FiniteElement<dim, spacedim>            &fe,
+      const ComponentMask                           &fe_mask,
+      const std::vector<unsigned int>               &fe_to_real,
+      std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_3rd_derivatives)
+        {
+          const unsigned int n_q_points = jacobian_3rd_derivatives.size();
 
-                // never touch any data for j=dim in case dim<spacedim, so
-                // it will always be zero as it was initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        jacobian_2nd_derivatives[point][i][j][l][m] =
-                          result[i][j][l][m];
-              }
-          }
-      }
-
-      /**
-       * Update the third derivative of the transformation from unit to real
-       * cell, the Jacobian hessians, pushed forward to the real cell
-       * coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_2nd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<4, spacedim>>
-          &jacobian_pushed_forward_2nd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_2nd_derivatives.size();
-
-            double tmp[spacedim][spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<3, dim> *third =
-                  &data.third_derivative(point + data_set, 0);
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<4, dim> *fourth =
+                &data.fourth_derivative(point + data_set, 0);
 
-                DerivativeForm<3, dim, spacedim> result;
+              DerivativeForm<4, dim, spacedim> result;
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            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)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        for (unsigned int m = 0; m < dim; ++m)
+                          for (unsigned int n = 0; n < dim; ++n)
+                            result[fe_to_real[comp_k]][j][l][m][n] +=
+                              (fourth[k][j][l][m][n] *
+                               data.local_dof_values[k]);
+                }
+
+              // never touch any data for j,l,m,n=dim in case
+              // dim<spacedim, so it will always be zero as it was
+              // initialized
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < dim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    for (unsigned int m = 0; m < dim; ++m)
+                      for (unsigned int n = 0; n < dim; ++n)
+                        jacobian_3rd_derivatives[point][i][j][l][m][n] =
+                          result[i][j][l][m][n];
+            }
+        }
+    }
+
+    /**
+     * Update the fourth derivative of the transformation from unit to real
+     * cell, the Jacobian hessian gradients, pushed forward to the real cell
+     * coordinates.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_update_jacobian_pushed_forward_3rd_derivatives(
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real,
+      std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
+    {
+      const UpdateFlags update_flags = data.update_each;
+      if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
+        {
+          const unsigned int n_q_points =
+            jacobian_pushed_forward_3rd_derivatives.size();
+
+          double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
+          for (unsigned int point = 0; point < n_q_points; ++point)
+            {
+              const Tensor<4, dim> *fourth =
+                &data.fourth_derivative(point + data_set, 0);
+
+              DerivativeForm<4, dim, spacedim> result;
 
-                // 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 k = 0; k < data.n_shape_functions; ++k)
+                {
+                  const unsigned int comp_k =
+                    fe.system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j = 0; j < dim; ++j)
+                      for (unsigned int l = 0; l < dim; ++l)
+                        for (unsigned int m = 0; m < dim; ++m)
+                          for (unsigned int n = 0; n < dim; ++n)
+                            result[fe_to_real[comp_k]][j][l][m][n] +=
+                              (fourth[k][j][l][m][n] *
+                               data.local_dof_values[k]);
+                }
+
+              // push-forward the j-coordinate
+              for (unsigned int i = 0; i < spacedim; ++i)
+                for (unsigned int j = 0; j < spacedim; ++j)
+                  for (unsigned int l = 0; l < dim; ++l)
+                    for (unsigned int m = 0; m < dim; ++m)
+                      for (unsigned int n = 0; n < dim; ++n)
                         {
-                          jacobian_pushed_forward_2nd_derivatives
-                            [point][i][j][l][m] =
-                              result[i][0][l][m] * data.covariant[point][j][0];
+                          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)
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [l][m] +=
-                              result[i][jr][l][m] *
-                              data.covariant[point][j][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)
+              // 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)
                         {
-                          tmp[i][j][l][m] =
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [0][m] *
-                            data.covariant[point][l][0];
+                          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)
-                            tmp[i][j][l][m] +=
-                              jacobian_pushed_forward_2nd_derivatives[point][i]
-                                                                     [j][lr]
-                                                                     [m] *
+                            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)
+              // 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)
                         {
-                          jacobian_pushed_forward_2nd_derivatives
-                            [point][i][j][l][m] =
-                              tmp[i][j][l][0] * data.covariant[point][m][0];
+                          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)
-                            jacobian_pushed_forward_2nd_derivatives[point][i][j]
-                                                                   [l][m] +=
-                              tmp[i][j][l][mr] * data.covariant[point][m][mr];
+                            tmp[i][j][l][m][n] +=
+                              jacobian_pushed_forward_3rd_derivatives[point][i]
+                                                                     [j][l][mr]
+                                                                     [n] *
+                              data.covariant[point][m][mr];
                         }
-              }
-          }
-      }
-
-      /**
-       * Update the fourth derivative of the transformation from unit to real
-       * cell, the Jacobian hessian gradients.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_3rd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                                &data,
-        const FiniteElement<dim, spacedim>            &fe,
-        const ComponentMask                           &fe_mask,
-        const std::vector<unsigned int>               &fe_to_real,
-        std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_3rd_derivatives)
-          {
-            const unsigned int n_q_points = jacobian_3rd_derivatives.size();
-
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<4, dim> *fourth =
-                  &data.fourth_derivative(point + data_set, 0);
 
-                DerivativeForm<4, dim, spacedim> result;
+              // 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];
+                        }
+            }
+        }
+    }
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
-                  {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.local_dof_values[k]);
-                  }
 
-                // never touch any data for j,l,m,n=dim in case
-                // dim<spacedim, so it will always be zero as it was
-                // initialized
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          jacobian_3rd_derivatives[point][i][j][l][m][n] =
-                            result[i][j][l][m][n];
-              }
-          }
-      }
-
-      /**
-       * Update the fourth derivative of the transformation from unit to real
-       * cell, the Jacobian hessian gradients, pushed forward to the real cell
-       * coordinates.
-       *
-       * Skip the computation if possible as indicated by the first argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_update_jacobian_pushed_forward_3rd_derivatives(
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        std::vector<Tensor<5, spacedim>>
-          &jacobian_pushed_forward_3rd_derivatives)
-      {
-        const UpdateFlags update_flags = data.update_each;
-        if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
-          {
-            const unsigned int n_q_points =
-              jacobian_pushed_forward_3rd_derivatives.size();
-
-            double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              {
-                const Tensor<4, dim> *fourth =
-                  &data.fourth_derivative(point + data_set, 0);
+    /**
+     * Depending on what information is called for in the update flags of the
+     * @p data object, compute the various pieces of information that is
+     * required by the fill_fe_face_values() and fill_fe_subface_values()
+     * functions.  This function simply unifies the work that would be done by
+     * those two functions.
+     *
+     * The resulting data is put into the @p output_data argument.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    maybe_compute_face_data(
+      const dealii::Mapping<dim, spacedim> &mapping,
+      const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
+      const unsigned int                                face_no,
+      const unsigned int                                subface_no,
+      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData &data,
+      internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+        &output_data)
+    {
+      const UpdateFlags update_flags = data.update_each;
 
-                DerivativeForm<4, dim, spacedim> result;
+      if (update_flags & update_boundary_forms)
+        {
+          const unsigned int n_q_points = output_data.boundary_forms.size();
+          if (update_flags & update_normal_vectors)
+            AssertDimension(output_data.normal_vectors.size(), n_q_points);
+          if (update_flags & update_JxW_values)
+            AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+          // map the unit tangentials to the real cell. checking for d!=dim-1
+          // eliminates compiler warnings regarding unsigned int expressions <
+          // 0.
+          for (unsigned int d = 0; d != dim - 1; ++d)
+            {
+              Assert(face_no + cell->n_faces() * d <
+                       data.unit_tangentials.size(),
+                     ExcInternalError());
+              Assert(
+                data.aux[d].size() <=
+                  data.unit_tangentials[face_no + cell->n_faces() * d].size(),
+                ExcInternalError());
+
+              mapping.transform(
+                make_array_view(
+                  data.unit_tangentials[face_no + cell->n_faces() * d]),
+                mapping_contravariant,
+                data,
+                make_array_view(data.aux[d]));
+            }
 
-                for (unsigned int k = 0; k < data.n_shape_functions; ++k)
+          // 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)
                   {
-                    const unsigned int comp_k =
-                      fe.system_to_component_index(k).first;
-                    if (fe_mask[comp_k])
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[fe_to_real[comp_k]][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.local_dof_values[k]);
+                    case 1:
+                      // in 1d, we don't have access to any of the data.aux
+                      // fields (because it has only dim-1 components), but we
+                      // can still compute the boundary form by simply looking
+                      // at the number of the face
+                      output_data.boundary_forms[i][0] =
+                        (face_no == 0 ? -1 : +1);
+                      break;
+                    case 2:
+                      output_data.boundary_forms[i] =
+                        cross_product_2d(data.aux[0][i]);
+                      break;
+                    case 3:
+                      output_data.boundary_forms[i] =
+                        cross_product_3d(data.aux[0][i], data.aux[1][i]);
+                      break;
+                    default:
+                      DEAL_II_NOT_IMPLEMENTED();
                   }
-
-                // push-forward the j-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            tmp[i][j][l][m][n] = result[i][0][l][m][n] *
-                                                 data.covariant[point][j][0];
-                            for (unsigned int jr = 1; jr < dim; ++jr)
-                              tmp[i][j][l][m][n] +=
-                                result[i][jr][l][m][n] *
-                                data.covariant[point][j][jr];
-                          }
-
-                // push-forward the l-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < dim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives
-                              [point][i][j][l][m][n] =
-                                tmp[i][j][0][m][n] *
-                                data.covariant[point][l][0];
-                            for (unsigned int lr = 1; lr < dim; ++lr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][m]
-                                                                     [n] +=
-                                tmp[i][j][lr][m][n] *
-                                data.covariant[point][l][lr];
-                          }
-
-                // push-forward the m-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < spacedim; ++m)
-                        for (unsigned int n = 0; n < dim; ++n)
-                          {
-                            tmp[i][j][l][m][n] =
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][0]
-                                                                     [n] *
-                              data.covariant[point][m][0];
-                            for (unsigned int mr = 1; mr < dim; ++mr)
-                              tmp[i][j][l][m][n] +=
-                                jacobian_pushed_forward_3rd_derivatives[point]
-                                                                       [i][j][l]
-                                                                       [mr][n] *
-                                data.covariant[point][m][mr];
-                          }
-
-                // push-forward the n-coordinate
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < spacedim; ++j)
-                    for (unsigned int l = 0; l < spacedim; ++l)
-                      for (unsigned int m = 0; m < spacedim; ++m)
-                        for (unsigned int n = 0; n < spacedim; ++n)
-                          {
-                            jacobian_pushed_forward_3rd_derivatives
-                              [point][i][j][l][m][n] =
-                                tmp[i][j][l][m][0] *
-                                data.covariant[point][n][0];
-                            for (unsigned int nr = 1; nr < dim; ++nr)
-                              jacobian_pushed_forward_3rd_derivatives[point][i]
-                                                                     [j][l][m]
-                                                                     [n] +=
-                                tmp[i][j][l][m][nr] *
-                                data.covariant[point][n][nr];
-                          }
-              }
-          }
-      }
-
-
-      /**
-       * Depending on what information is called for in the update flags of the
-       * @p data object, compute the various pieces of information that is
-       * required by the fill_fe_face_values() and fill_fe_subface_values()
-       * functions.  This function simply unifies the work that would be done by
-       * those two functions.
-       *
-       * The resulting data is put into the @p output_data argument.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      maybe_compute_face_data(
-        const dealii::Mapping<dim, spacedim> &mapping,
-        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
-                                                         &cell,
-        const unsigned int                                face_no,
-        const unsigned int                                subface_no,
-        const typename QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data,
-        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-          &output_data)
-      {
-        const UpdateFlags update_flags = data.update_each;
-
-        if (update_flags & update_boundary_forms)
-          {
-            const unsigned int n_q_points = output_data.boundary_forms.size();
-            if (update_flags & update_normal_vectors)
-              AssertDimension(output_data.normal_vectors.size(), n_q_points);
-            if (update_flags & update_JxW_values)
-              AssertDimension(output_data.JxW_values.size(), n_q_points);
-
-            // map the unit tangentials to the real cell. checking for d!=dim-1
-            // eliminates compiler warnings regarding unsigned int expressions <
-            // 0.
-            for (unsigned int d = 0; d != dim - 1; ++d)
-              {
-                Assert(face_no + cell->n_faces() * d <
-                         data.unit_tangentials.size(),
-                       ExcInternalError());
-                Assert(
-                  data.aux[d].size() <=
-                    data.unit_tangentials[face_no + cell->n_faces() * d].size(),
-                  ExcInternalError());
-
-                mapping.transform(
-                  make_array_view(
-                    data.unit_tangentials[face_no + cell->n_faces() * d]),
-                  mapping_contravariant,
-                  data,
-                  make_array_view(data.aux[d]));
-              }
-
-            // if dim==spacedim, we can use the unit tangentials to compute the
-            // boundary form by simply taking the cross product
-            if (dim == spacedim)
-              {
-                for (unsigned int i = 0; i < n_q_points; ++i)
-                  switch (dim)
+            }
+          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)
                     {
-                      case 1:
-                        // in 1d, we don't have access to any of the data.aux
-                        // fields (because it has only dim-1 components), but we
-                        // can still compute the boundary form by simply looking
-                        // at the number of the face
-                        output_data.boundary_forms[i][0] =
-                          (face_no == 0 ? -1 : +1);
-                        break;
-                      case 2:
-                        output_data.boundary_forms[i] =
-                          cross_product_2d(data.aux[0][i]);
-                        break;
-                      case 3:
-                        output_data.boundary_forms[i] =
-                          cross_product_3d(data.aux[0][i], data.aux[1][i]);
-                        break;
-                      default:
-                        DEAL_II_NOT_IMPLEMENTED();
+                      // 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();
                     }
-              }
-            else //(dim < spacedim)
-              {
-                // in the codim-one case, the boundary form results from the
-                // cross product of all the face tangential vectors and the cell
-                // normal vector
-                //
-                // to compute the cell normal, use the same method used in
-                // fill_fe_values for cells above
-                AssertDimension(data.contravariant.size(), n_q_points);
-
-                for (unsigned int point = 0; point < n_q_points; ++point)
-                  {
-                    if (dim == 1)
-                      {
-                        // J is a tangent vector
-                        output_data.boundary_forms[point] =
-                          data.contravariant[point].transpose()[0];
-                        output_data.boundary_forms[point] /=
-                          (face_no == 0 ? -1. : +1.) *
-                          output_data.boundary_forms[point].norm();
-                      }
-
-                    if (dim == 2)
-                      {
-                        const DerivativeForm<1, spacedim, dim> DX_t =
-                          data.contravariant[point].transpose();
-
-                        Tensor<1, spacedim> cell_normal =
-                          cross_product_3d(DX_t[0], DX_t[1]);
-                        cell_normal /= cell_normal.norm();
-
-                        // then compute the face normal from the face tangent
-                        // and the cell normal:
-                        output_data.boundary_forms[point] =
-                          cross_product_3d(data.aux[0][point], cell_normal);
-                      }
-                  }
-              }
 
-            if (update_flags & (update_normal_vectors | update_JxW_values))
-              for (unsigned int i = 0; i < output_data.boundary_forms.size();
-                   ++i)
-                {
-                  if (update_flags & update_JxW_values)
+                  if (dim == 2)
                     {
-                      output_data.JxW_values[i] =
-                        output_data.boundary_forms[i].norm() *
-                        data.quadrature_weights[i + data_set];
+                      const DerivativeForm<1, spacedim, dim> DX_t =
+                        data.contravariant[point].transpose();
 
-                      if (subface_no != numbers::invalid_unsigned_int)
-                        {
-                          // TODO
-                          const double area_ratio =
-                            GeometryInfo<dim>::subface_ratio(
-                              cell->subface_case(face_no), subface_no);
-                          output_data.JxW_values[i] *= area_ratio;
-                        }
-                    }
+                      Tensor<1, spacedim> cell_normal =
+                        cross_product_3d(DX_t[0], DX_t[1]);
+                      cell_normal /= cell_normal.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());
+                      // 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);
+                    }
                 }
-          }
-      }
-
-      /**
-       * Do the work of MappingFEField::fill_fe_face_values() and
-       * MappingFEField::fill_fe_subface_values() in a generic way, using the
-       * 'data_set' to differentiate whether we will work on a face (and if so,
-       * which one) or subface.
-       */
-      template <int dim, int spacedim, typename VectorType>
-      void
-      do_fill_fe_face_values(
-        const dealii::Mapping<dim, spacedim> &mapping,
-        const typename dealii::Triangulation<dim, spacedim>::cell_iterator
-                                                                 &cell,
-        const unsigned int                                        face_no,
-        const unsigned int                                        subface_no,
-        const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData                     &data,
-        const FiniteElement<dim, spacedim> &fe,
-        const ComponentMask                &fe_mask,
-        const std::vector<unsigned int>    &fe_to_real,
-        internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-          &output_data)
-      {
-        maybe_compute_q_points<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.quadrature_points);
-
-        maybe_update_Jacobians<dim, spacedim, VectorType>(
-          data_set, data, fe, fe_mask, fe_to_real);
-
-        const UpdateFlags  update_flags = data.update_each;
-        const unsigned int n_q_points   = data.contravariant.size();
-
-        if (update_flags & update_jacobians)
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.jacobians[point] = data.contravariant[point];
+            }
 
-        if (update_flags & update_inverse_jacobians)
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            output_data.inverse_jacobians[point] =
-              data.covariant[point].transpose();
+          if (update_flags & (update_normal_vectors | update_JxW_values))
+            for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
+              {
+                if (update_flags & update_JxW_values)
+                  {
+                    output_data.JxW_values[i] =
+                      output_data.boundary_forms[i].norm() *
+                      data.quadrature_weights[i + data_set];
 
-        maybe_update_jacobian_grads<dim, spacedim, VectorType>(
-          data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
+                    if (subface_no != numbers::invalid_unsigned_int)
+                      {
+                        // TODO
+                        const double area_ratio =
+                          GeometryInfo<dim>::subface_ratio(
+                            cell->subface_case(face_no), subface_no);
+                        output_data.JxW_values[i] *= area_ratio;
+                      }
+                  }
 
-        maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_grads);
+                if (update_flags & update_normal_vectors)
+                  output_data.normal_vectors[i] =
+                    Point<spacedim>(output_data.boundary_forms[i] /
+                                    output_data.boundary_forms[i].norm());
+              }
+        }
+    }
 
-        maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_2nd_derivatives);
+    /**
+     * Do the work of MappingFEField::fill_fe_face_values() and
+     * MappingFEField::fill_fe_subface_values() in a generic way, using the
+     * 'data_set' to differentiate whether we will work on a face (and if so,
+     * which one) or subface.
+     */
+    template <int dim, int spacedim, typename VectorType>
+    void
+    do_fill_fe_face_values(
+      const dealii::Mapping<dim, spacedim> &mapping,
+      const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
+      const unsigned int                                        face_no,
+      const unsigned int                                        subface_no,
+      const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData                     &data,
+      const FiniteElement<dim, spacedim> &fe,
+      const ComponentMask                &fe_mask,
+      const std::vector<unsigned int>    &fe_to_real,
+      internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+        &output_data)
+    {
+      maybe_compute_q_points<dim, spacedim, VectorType>(
+        data_set, data, fe, fe_mask, fe_to_real, output_data.quadrature_points);
 
-        maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_2nd_derivatives);
+      maybe_update_Jacobians<dim, spacedim, VectorType>(
+        data_set, data, fe, fe_mask, fe_to_real);
 
-        maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_3rd_derivatives);
+      const UpdateFlags  update_flags = data.update_each;
+      const unsigned int n_q_points   = data.contravariant.size();
 
-        maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
-                                                             spacedim,
-                                                             VectorType>(
-          data_set,
-          data,
-          fe,
-          fe_mask,
-          fe_to_real,
-          output_data.jacobian_pushed_forward_3rd_derivatives);
+      if (update_flags & update_jacobians)
+        for (unsigned int point = 0; point < n_q_points; ++point)
+          output_data.jacobians[point] = data.contravariant[point];
 
-        maybe_compute_face_data<dim, spacedim, VectorType>(
-          mapping, cell, face_no, subface_no, data_set, data, output_data);
-      }
-    } // namespace
-  }   // namespace MappingFEFieldImplementation
+      if (update_flags & update_inverse_jacobians)
+        for (unsigned int point = 0; point < n_q_points; ++point)
+          output_data.inverse_jacobians[point] =
+            data.covariant[point].transpose();
+
+      maybe_update_jacobian_grads<dim, spacedim, VectorType>(
+        data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
+
+      maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
+        data_set,
+        data,
+        fe,
+        fe_mask,
+        fe_to_real,
+        output_data.jacobian_pushed_forward_grads);
+
+      maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
+        data_set,
+        data,
+        fe,
+        fe_mask,
+        fe_to_real,
+        output_data.jacobian_2nd_derivatives);
+
+      maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
+                                                           spacedim,
+                                                           VectorType>(
+        data_set,
+        data,
+        fe,
+        fe_mask,
+        fe_to_real,
+        output_data.jacobian_pushed_forward_2nd_derivatives);
+
+      maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
+        data_set,
+        data,
+        fe,
+        fe_mask,
+        fe_to_real,
+        output_data.jacobian_3rd_derivatives);
+
+      maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
+                                                           spacedim,
+                                                           VectorType>(
+        data_set,
+        data,
+        fe,
+        fe_mask,
+        fe_to_real,
+        output_data.jacobian_pushed_forward_3rd_derivatives);
+
+      maybe_compute_face_data<dim, spacedim, VectorType>(
+        mapping, cell, face_no, subface_no, data_set, data, output_data);
+    }
+  } // namespace MappingFEFieldImplementation
 } // namespace internal
 
 
@@ -1957,125 +1935,117 @@ namespace internal
 {
   namespace MappingFEFieldImplementation
   {
-    namespace
+    template <int dim, int spacedim, int rank, typename VectorType>
+    void
+    transform_fields(
+      const ArrayView<const Tensor<rank, dim>>                &input,
+      const MappingKind                                        mapping_kind,
+      const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
+      const ArrayView<Tensor<rank, spacedim>>                 &output)
     {
-      template <int dim, int spacedim, int rank, typename VectorType>
-      void
-      transform_fields(
-        const ArrayView<const Tensor<rank, dim>>                &input,
-        const MappingKind                                        mapping_kind,
-        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-        const ArrayView<Tensor<rank, spacedim>>                 &output)
-      {
-        AssertDimension(input.size(), output.size());
-        Assert((dynamic_cast<
-                  const typename dealii::
-                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
-                  &mapping_data) != nullptr),
-               ExcInternalError());
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data = static_cast<
-            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-              InternalData &>(mapping_data);
-
-        switch (mapping_kind)
-          {
-            case mapping_contravariant:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
+      AssertDimension(input.size(), output.size());
+      Assert((dynamic_cast<
+                const typename dealii::
+                  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
+                &mapping_data) != nullptr),
+             ExcInternalError());
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData &data = static_cast<
+          const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+            InternalData &>(mapping_data);
+
+      switch (mapping_kind)
+        {
+          case mapping_contravariant:
+            {
+              Assert(data.update_each & update_contravariant_transformation,
+                     typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                       "update_contravariant_transformation"));
+
+              for (unsigned int i = 0; i < output.size(); ++i)
+                output[i] =
+                  apply_transformation(data.contravariant[i], input[i]);
 
-                for (unsigned int i = 0; i < output.size(); ++i)
+              return;
+            }
+
+          case mapping_piola:
+            {
+              Assert(data.update_each & update_contravariant_transformation,
+                     typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                       "update_contravariant_transformation"));
+              Assert(data.update_each & update_volume_elements,
+                     typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                       "update_volume_elements"));
+              Assert(rank == 1, ExcMessage("Only for rank 1"));
+              for (unsigned int i = 0; i < output.size(); ++i)
+                {
                   output[i] =
                     apply_transformation(data.contravariant[i], input[i]);
+                  output[i] /= data.volume_elements[i];
+                }
+              return;
+            }
 
-                return;
-              }
 
-            case mapping_piola:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
-                Assert(
-                  data.update_each & update_volume_elements,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_volume_elements"));
-                Assert(rank == 1, ExcMessage("Only for rank 1"));
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  {
-                    output[i] =
-                      apply_transformation(data.contravariant[i], input[i]);
-                    output[i] /= data.volume_elements[i];
-                  }
-                return;
-              }
+          // We still allow this operation as in the
+          // reference cell Derivatives are Tensor
+          // rather than DerivativeForm
+          case mapping_covariant:
+            {
+              Assert(data.update_each & update_contravariant_transformation,
+                     typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                       "update_contravariant_transformation"));
 
+              for (unsigned int i = 0; i < output.size(); ++i)
+                output[i] = apply_transformation(data.covariant[i], input[i]);
 
-            // 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"));
+              return;
+            }
 
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  output[i] = apply_transformation(data.covariant[i], input[i]);
+          default:
+            DEAL_II_NOT_IMPLEMENTED();
+        }
+    }
 
-                return;
-              }
 
-            default:
-              DEAL_II_NOT_IMPLEMENTED();
-          }
-      }
-
-
-      template <int dim, int spacedim, int rank, typename VectorType>
-      void
-      transform_differential_forms(
-        const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
-        const MappingKind                                        mapping_kind,
-        const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
-        const ArrayView<Tensor<rank + 1, spacedim>>             &output)
-      {
-        AssertDimension(input.size(), output.size());
-        Assert((dynamic_cast<
-                  const typename dealii::
-                    MappingFEField<dim, spacedim, VectorType>::InternalData *>(
-                  &mapping_data) != nullptr),
-               ExcInternalError());
-        const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-          InternalData &data = static_cast<
-            const typename dealii::MappingFEField<dim, spacedim, VectorType>::
-              InternalData &>(mapping_data);
-
-        switch (mapping_kind)
-          {
-            case mapping_covariant:
-              {
-                Assert(
-                  data.update_each & update_contravariant_transformation,
-                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
-                    "update_contravariant_transformation"));
+    template <int dim, int spacedim, int rank, typename VectorType>
+    void
+    transform_differential_forms(
+      const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
+      const MappingKind                                           mapping_kind,
+      const typename Mapping<dim, spacedim>::InternalDataBase    &mapping_data,
+      const ArrayView<Tensor<rank + 1, spacedim>>                &output)
+    {
+      AssertDimension(input.size(), output.size());
+      Assert((dynamic_cast<
+                const typename dealii::
+                  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
+                &mapping_data) != nullptr),
+             ExcInternalError());
+      const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+        InternalData &data = static_cast<
+          const typename dealii::MappingFEField<dim, spacedim, VectorType>::
+            InternalData &>(mapping_data);
+
+      switch (mapping_kind)
+        {
+          case mapping_covariant:
+            {
+              Assert(data.update_each & update_contravariant_transformation,
+                     typename FEValuesBase<dim>::ExcAccessToUninitializedField(
+                       "update_contravariant_transformation"));
 
-                for (unsigned int i = 0; i < output.size(); ++i)
-                  output[i] = apply_transformation(data.covariant[i], input[i]);
+              for (unsigned int i = 0; i < output.size(); ++i)
+                output[i] = apply_transformation(data.covariant[i], input[i]);
 
-                return;
-              }
-            default:
-              DEAL_II_NOT_IMPLEMENTED();
-          }
-      }
-    } // namespace
-  }   // namespace MappingFEFieldImplementation
+              return;
+            }
+          default:
+            DEAL_II_NOT_IMPLEMENTED();
+        }
+    }
+  } // namespace MappingFEFieldImplementation
 } // namespace internal
 
 

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.