From bcd0c3ad97f575df2fab3947825eed92b7c76c48 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 12 Apr 2018 13:37:27 +0200
Subject: [PATCH] Add methods to compute the face data in MappingInfo.

---
 include/deal.II/matrix_free/mapping_info.h    |  30 +-
 .../matrix_free/mapping_info.templates.h      | 780 +++++++++++++++++-
 .../matrix_free/matrix_free.templates.h       |  16 +-
 3 files changed, 820 insertions(+), 6 deletions(-)

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index ed4f4ad8d1..37bc22c860 100644
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -25,6 +25,7 @@
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/mapping.h>
 #include <deal.II/matrix_free/helper_functions.h>
+#include <deal.II/matrix_free/face_info.h>
 
 #include <memory>
 
@@ -289,10 +290,14 @@ namespace internal
        */
       void initialize (const dealii::Triangulation<dim>                &tria,
                        const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
+                       const FaceInfo<VectorizedArray<Number>::n_array_elements> &faces,
                        const std::vector<unsigned int>         &active_fe_index,
                        const Mapping<dim>                      &mapping,
                        const std::vector<dealii::hp::QCollection<1> >  &quad,
-                       const UpdateFlags                        update_flags_cells);
+                       const UpdateFlags                        update_flags_cells,
+                       const UpdateFlags                        update_flags_boundary_faces,
+                       const UpdateFlags                        update_flags_inner_faces,
+                       const UpdateFlags                        update_flags_faces_by_cells);
 
       /**
        * Return the type of a given cell as detected during initialization.
@@ -361,6 +366,29 @@ namespace internal
                              const std::vector<dealii::hp::QCollection<1> >  &quad,
                              const UpdateFlags                        update_flags_cells);
 
+      /**
+       * Computes the information in the given faces, called within
+       * initialize.
+       */
+      void initialize_faces (const dealii::Triangulation<dim>        &tria,
+                             const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+                             const std::vector<FaceToCellTopology<VectorizedArray<Number>::n_array_elements> > &faces,
+                             const Mapping<dim>                      &mapping,
+                             const std::vector<dealii::hp::QCollection<1> >  &quad,
+                             const UpdateFlags                        update_flags_boundary_faces,
+                             const UpdateFlags                        update_flags_inner_faces);
+
+      /**
+       * Computes the information in the given faces, called within
+       * initialize.
+       */
+      void initialize_faces_by_cells
+      (const dealii::Triangulation<dim>        &tria,
+       const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+       const Mapping<dim>                      &mapping,
+       const std::vector<dealii::hp::QCollection<1> > &quad,
+       const UpdateFlags                        update_flags_faces_by_cells);
+
       /**
        * Helper function to determine which update flags must be set in the
        * internal functions to initialize all data as requested by the user.
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 8bd1f58e5f..796e2e9c87 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -261,16 +261,23 @@ namespace internal
     MappingInfo<dim,Number>::initialize
     (const dealii::Triangulation<dim>                          &tria,
      const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
+     const FaceInfo<VectorizedArray<Number>::n_array_elements> &face_info,
      const std::vector<unsigned int>                           &active_fe_index,
      const Mapping<dim>                                        &mapping,
      const std::vector<dealii::hp::QCollection<1> >            &quad,
-     const UpdateFlags                                          update_flags_cells)
+     const UpdateFlags                                          update_flags_cells,
+     const UpdateFlags                                          update_flags_boundary_faces,
+     const UpdateFlags                                          update_flags_inner_faces,
+     const UpdateFlags                                          update_flags_faces_by_cells)
     {
       clear();
 
       // Could call these functions in parallel, but not useful because the
       // work inside is nicely split up already
       initialize_cells(tria, cells, active_fe_index, mapping, quad, update_flags_cells);
+      initialize_faces(tria, cells, face_info.faces, mapping, quad,
+                       update_flags_boundary_faces, update_flags_inner_faces);
+      initialize_faces_by_cells(tria, cells, mapping, quad, update_flags_faces_by_cells);
     }
 
 
@@ -1069,6 +1076,777 @@ namespace internal
 
 
 
+    /* ------------------------- initialization of faces ------------------- */
+
+    // Anonymous namespace with implementation of extraction of values on cell
+    // range
+    namespace
+    {
+      template <int dim, typename Number>
+      struct CompressedFaceData
+      {
+        // Constructor. As a scaling factor for the FPArrayComparator, we
+        // select the inverse of the Jacobian (not the Jacobian as in the
+        // CompressedCellData) and add another factor of 512 to account for
+        // some roundoff effects.
+        CompressedFaceData(const Number jacobian_size)
+          :
+          data(FPArrayComparator<Number>(512./jacobian_size)),
+          jacobian_size(jacobian_size)
+        {}
+
+        // Store the Jacobians on both sides of a face (2*(dim*dim) entries),
+        // the normal vector (dim entries), and the Jacobian determinant on
+        // the face (first tensor) for each entry in the vectorized array
+        // (inner tensor). We cannot choose a VectorizedArray type directly
+        // because std::map does not provide the necessary alignment upon
+        // memory allocation.
+        std::map<Tensor<1,2*dim *dim+dim+1,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> >
+        , unsigned int, FPArrayComparator<Number> > data;
+
+        // Store the scaling factor
+        const Number jacobian_size;
+      };
+
+
+
+      // We always put the derivative normal to the face in the last slot for
+      // simpler unit cell gradient computations. This function reorders the
+      // indices of the Jacobian appropriately.
+      template <int dim>
+      unsigned int
+      reorder_face_derivative_indices(const unsigned int face_no,
+                                      const unsigned int index)
+      {
+        Assert(index < dim, ExcInternalError());
+        if (dim == 3)
+          {
+            unsigned int table [3][3] = {{1, 2, 0}, {2, 0, 1}, {0, 1, 2}};
+            return table[face_no/2][index];
+          }
+        else if (dim == 2)
+          {
+            unsigned int table [2][2] = {{1, 0}, {0, 1}};
+            return table[face_no/2][index];
+          }
+        else if (dim == 1)
+          return 0;
+        else
+          Assert(false, ExcNotImplemented("Not possible in dim=" + std::to_string(dim)));
+      }
+
+
+
+      template <int dim, typename Number>
+      void
+      initialize_face_range
+      (const std::pair<unsigned int,unsigned int>                 face_range,
+       const dealii::Triangulation<dim>                          &tria,
+       const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
+       const std::vector<FaceToCellTopology<VectorizedArray<Number>::n_array_elements> > &faces,
+       const Mapping<dim>                                        &mapping,
+       const UpdateFlags                                          update_flags_boundary,
+       const UpdateFlags                                          update_flags_inner,
+       MappingInfo<dim,Number>                                   &mapping_info,
+       std::pair<std::vector<MappingInfoStorage<dim-1,dim,Number> >,
+       CompressedFaceData<dim,Number> >                          &data)
+      {
+        FE_Nothing<dim> dummy_fe;
+
+        std::vector<std::vector<std::shared_ptr<FEFaceValues<dim> > > >
+        fe_face_values_container(mapping_info.face_data.size());
+        for (unsigned int my_q=0; my_q<mapping_info.face_data.size(); ++my_q)
+          fe_face_values_container[my_q].resize(mapping_info.face_data[my_q].
+                                                descriptor.size());
+
+        std::vector<std::vector<std::shared_ptr<FEFaceValues<dim> > > >
+        fe_boundary_face_values_container(mapping_info.face_data.size());
+        for (unsigned int my_q=0; my_q<mapping_info.face_data.size(); ++my_q)
+          fe_boundary_face_values_container[my_q].resize(mapping_info.face_data[my_q].
+                                                         descriptor.size());
+
+        std::vector<std::vector<std::shared_ptr<FESubfaceValues<dim> > > >
+        fe_subface_values_container(mapping_info.face_data.size());
+        for (unsigned int my_q=0; my_q<mapping_info.face_data.size(); ++my_q)
+          fe_subface_values_container[my_q].resize(mapping_info.face_data[my_q].
+                                                   descriptor.size());
+
+        LocalData<dim,Number> face_data (get_jacobian_size(tria));
+
+        const unsigned int end_face =
+          std::min(std::size_t(face_range.second), faces.size());
+        for (unsigned int face=face_range.first; face<end_face; ++face)
+          for (unsigned int my_q=0; my_q<mapping_info.face_data.size(); ++my_q)
+            {
+              // currently only non-hp case...
+              Assert(mapping_info.face_data[my_q].descriptor.size() == 1,
+                     ExcNotImplemented());
+              const Quadrature<dim-1> &quadrature =
+                mapping_info.face_data[my_q].descriptor[0].quadrature;
+
+              const bool is_boundary_face =
+                faces[face].cells_exterior[0] == numbers::invalid_unsigned_int;
+
+              if (is_boundary_face &&
+                  fe_boundary_face_values_container[my_q][0].get() == 0)
+                fe_boundary_face_values_container[my_q][0].reset
+                (new FEFaceValues<dim> (mapping, dummy_fe, quadrature,
+                                        update_flags_boundary));
+              else if (fe_face_values_container[my_q][0].get() == 0)
+                fe_face_values_container[my_q][0].reset
+                (new FEFaceValues<dim> (mapping, dummy_fe, quadrature,
+                                        update_flags_inner));
+
+              FEFaceValues<dim> &fe_face_values = is_boundary_face ?
+                                                  *fe_boundary_face_values_container[my_q][0] :
+                                                  *fe_face_values_container[my_q][0];
+              const unsigned int n_q_points = fe_face_values.n_quadrature_points;
+              face_data.resize(n_q_points);
+
+              bool normal_is_similar = true;
+              bool JxW_is_similar = true;
+              bool cell_is_cartesian = true;
+              for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                {
+                  Tensor<2,dim> jacobian_0;
+                  double compare_norm_jac = 1.;
+                  if (faces[face].cells_interior[v] != numbers::invalid_unsigned_int)
+                    {
+                      typename dealii::Triangulation<dim>::cell_iterator
+                      cell_it (&tria, cells[faces[face].cells_interior[v]].first,
+                               cells[faces[face].cells_interior[v]].second);
+
+                      fe_face_values.reinit(cell_it, faces[face].interior_face_no);
+
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        {
+                          if (std::abs(fe_face_values.JxW(q)*quadrature.weight(0)-
+                                       fe_face_values.JxW(0)*quadrature.weight(q))
+                              > 2048. * std::numeric_limits<double>::epsilon() *
+                              fe_face_values.JxW(0)*quadrature.weight(q))
+                            JxW_is_similar = false;
+                          face_data.JxW_values[q][v] = fe_face_values.JxW(q);
+
+                          DerivativeForm<1,dim,dim> inv_jac = fe_face_values.jacobian(q).covariant_form();
+                          Tensor<1,dim> normal = fe_face_values.normal_vector(q);
+
+                          // Filter out very small values in normal. No need
+                          // to re-normalize because these values cannot enter
+                          // the norm significantly: Total size is 1 but 1e-13
+                          // squared is 1e-26.
+                          for (unsigned int d=0; d<dim; ++d)
+                            if (std::abs(normal[d]) < 1024. *
+                                std::numeric_limits<double>::epsilon())
+                              normal[d] = 0.;
+
+                          if (q==0)
+                            {
+                              jacobian_0 = inv_jac;
+                              compare_norm_jac = jacobian_0.norm();
+                            }
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              {
+                                if (std::abs(inv_jac[d][e]-jacobian_0[d][e])
+                                    > 2048. * std::numeric_limits<double>::epsilon() *
+                                    compare_norm_jac)
+                                  JxW_is_similar = false;
+                                const unsigned int ee=
+                                  reorder_face_derivative_indices<dim>(faces[face].interior_face_no,e);
+                                face_data.general_jac[q][d][e][v] = inv_jac[d][ee];
+                              }
+
+                          for (unsigned int d=0; d<dim; ++d)
+                            {
+                              face_data.normal_vectors[q][d][v] = normal[d];
+                              if (std::abs(normal[d]-
+                                           fe_face_values.normal_vector(0)[d]) > 1024. *
+                                  std::numeric_limits<double>::epsilon())
+                                normal_is_similar = false;
+                            }
+                          if (std::abs(std::abs(normal[faces[face].interior_face_no/2])-
+                                       1.) > 1024. * std::numeric_limits<double>::epsilon())
+                            cell_is_cartesian = false;
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              if (e != d &&
+                                  std::abs(inv_jac[d][e]) >
+                                  2048*std::numeric_limits<double>::epsilon() *
+                                  compare_norm_jac)
+                                cell_is_cartesian = false;
+
+                          if (fe_face_values.get_update_flags() & update_quadrature_points)
+                            for (unsigned int d=0; d<dim; ++d)
+                              face_data.quadrature_points[q][d][v] = fe_face_values.quadrature_point(q)[d];
+                        }
+                    }
+                  // Fill up with data of the zeroth component to avoid
+                  // false negatives when checking for similarities
+                  else
+                    {
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        {
+                          face_data.JxW_values[q][v] = face_data.JxW_values[q][0];
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              {
+                                face_data.general_jac[q][d][e][v] =
+                                  face_data.general_jac[q][d][e][0];
+                                face_data.normal_vectors[q][d][v] =
+                                  face_data.normal_vectors[q][d][0];
+                              }
+                        }
+                    }
+                  if (is_boundary_face == false &&
+                      faces[face].cells_exterior[v] != numbers::invalid_unsigned_int)
+                    {
+                      typename dealii::Triangulation<dim>::cell_iterator
+                      cell_it (&tria, cells[faces[face].cells_exterior[v]].first,
+                               cells[faces[face].cells_exterior[v]].second);
+
+                      const FEValuesBase<dim> *actual_fe_face_values = 0;
+                      if (faces[face].subface_index >=
+                          GeometryInfo<dim>::max_children_per_cell)
+                        {
+                          fe_face_values.reinit(cell_it, faces[face].exterior_face_no);
+                          actual_fe_face_values = &fe_face_values;
+                        }
+                      else
+                        {
+                          if (fe_subface_values_container[my_q][0].get() == 0)
+                            fe_subface_values_container[my_q][0].reset
+                            (new FESubfaceValues<dim> (mapping, dummy_fe, quadrature,
+                                                       update_flags_inner));
+                          fe_subface_values_container[my_q][0]->
+                          reinit(cell_it, faces[face].exterior_face_no,
+                                 faces[face].subface_index);
+                          actual_fe_face_values = fe_subface_values_container[my_q][0].get();
+                        }
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        {
+                          DerivativeForm<1,dim,dim> inv_jac = actual_fe_face_values->jacobian(q).covariant_form();
+                          if (q==0)
+                            {
+                              jacobian_0 = inv_jac;
+                              compare_norm_jac = jacobian_0.norm();
+                            }
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              {
+                                if (std::abs(inv_jac[d][e]-jacobian_0[d][e])
+                                    > 2048. * std::numeric_limits<double>::epsilon() *
+                                    compare_norm_jac)
+                                  JxW_is_similar = false;
+                                const unsigned int ee=
+                                  reorder_face_derivative_indices<dim>(faces[face].exterior_face_no,e);
+                                face_data.general_jac[n_q_points+q][d][e][v] = inv_jac[d][ee];
+                              }
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              if (e != d &&
+                                  std::abs(inv_jac[d][e]) >
+                                  2048*std::numeric_limits<double>::epsilon() *
+                                  compare_norm_jac)
+                                cell_is_cartesian = false;
+                        }
+                    }
+                  // Fill up with 'known' values
+                  else if (is_boundary_face == false)
+                    {
+                      Assert(faces[face].cells_exterior[0] != numbers::invalid_unsigned_int,
+                             ExcInternalError());
+                      for (unsigned int q=0; q<n_q_points; ++q)
+                        for (unsigned int d=0; d<dim; ++d)
+                          for (unsigned int e=0; e<dim; ++e)
+                            face_data.general_jac[n_q_points+q][d][e][v] = face_data.general_jac[n_q_points+q][d][e][0];
+                    }
+                  // If boundary face, simply set the data to zero (will not
+                  // be used). Note that faces over periodic boundary
+                  // conditions will be treated as interior ones in this setup.
+                  else
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      for (unsigned int d=0; d<dim; ++d)
+                        for (unsigned int e=0; e<dim; ++e)
+                          face_data.general_jac[n_q_points+q][d][e][v] = 0.;
+                }
+
+              // check if face is affine or at least if it is flat
+              // (i.e., all normal vectors are the same)
+              if (my_q == 0)
+                {
+                  GeometryType face_type = affine;
+                  if (JxW_is_similar == false)
+                    face_type = flat_faces;
+                  if (normal_is_similar == false)
+                    face_type = general;
+                  if (face_type == affine && cell_is_cartesian)
+                    face_type = cartesian;
+                  mapping_info.face_type[face] = face_type;
+                }
+
+              // Fill in quadrature points
+              if (fe_face_values.get_update_flags() & update_quadrature_points)
+                {
+                  data.first[my_q].quadrature_point_offsets.
+                  push_back(data.first[my_q].quadrature_points.size());
+                  if (fe_face_values.get_update_flags() & update_quadrature_points)
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      data.first[my_q].quadrature_points.push_back(face_data.quadrature_points[q]);
+                }
+
+              typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
+              unsigned int insert_position = data.first[my_q].JxW_values.size();
+
+              // Fill in JxW values, apply compression
+              if (mapping_info.face_type[face] <= affine)
+                {
+                  if (my_q == 0)
+                    {
+                      // find out if we already had the same JxW values before
+                      std::pair<Tensor<1,2*dim *dim+dim+1,VEC_ARRAY>,unsigned int> new_entry;
+                      new_entry.second = data.second.data.size();
+                      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                        new_entry.first[2*dim*dim+dim][v] = face_data.JxW_values[0][v]/quadrature.weight(0)/
+                                                            Utilities::fixed_power<dim>(face_data.jac_size);
+
+                      new_entry.second = data.second.data.size();
+                      for (unsigned int d=0; d<dim; ++d)
+                        for (unsigned int e=0; e<dim; ++e)
+                          for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                            new_entry.first[d*dim+e][v] = face_data.general_jac[0][d][e][v];
+                      if (is_boundary_face == false)
+                        for (unsigned int d=0; d<dim; ++d)
+                          for (unsigned int e=0; e<dim; ++e)
+                            for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                              new_entry.first[dim*dim+d*dim+e][v] = face_data.general_jac[n_q_points][d][e][v];
+                      // we need to add the normal vector here because we
+                      // store both the inverse jacobian and the normal vector
+                      // times the jacobian; of course, there will be
+                      // different values in their product for normal vectors
+                      // oriented in different ways (the memory saving is
+                      // still significant); we need to divide by the jacobian
+                      // size to get the right scaling
+                      for (unsigned int d=0; d<dim; ++d)
+                        for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                          new_entry.first[2*dim*dim+d][v] = face_data.normal_vectors[0][d][v] / face_data.jac_size;
+
+                      insert_position = data.second.data.insert(new_entry).first->second;
+                    }
+                  else
+                    insert_position = data.first[0].data_index_offsets[face];
+                }
+              data.first[my_q].data_index_offsets.push_back(insert_position);
+              if (mapping_info.face_type[face] > affine)
+                {
+                  for (unsigned int q=0; q<n_q_points; ++q)
+                    data.first[my_q].JxW_values.push_back(face_data.JxW_values[q]);
+                  for (unsigned int q=0; q<n_q_points; ++q)
+                    data.first[my_q].normal_vectors.push_back(face_data.normal_vectors[q]);
+                  for (unsigned int q=0; q<n_q_points; ++q)
+                    data.first[my_q].jacobians[0].push_back(face_data.general_jac[q]);
+                  if (is_boundary_face == false)
+                    for (unsigned int q=0; q<n_q_points; ++q)
+                      data.first[my_q].jacobians[1].push_back(face_data.general_jac[n_q_points+q]);
+                }
+            }
+      }
+
+      template <int dim, typename Number>
+      void
+      compute_normal_times_jacobian (const unsigned int           first_face,
+                                     const unsigned int           last_face,
+                                     const std::vector<GeometryType> &face_type,
+                                     const std::vector<FaceToCellTopology<VectorizedArray<Number>::n_array_elements> > &faces,
+                                     MappingInfoStorage<dim-1,dim,Number> &data_faces)
+      {
+        for (unsigned int face = first_face; face < last_face; ++face)
+          {
+            const bool is_boundary_face = faces[face].cells_exterior[0] == numbers::invalid_unsigned_int;
+            const unsigned int n_q_points_work =
+              face_type[face] > affine ? data_faces.descriptor[0].n_q_points : 1;
+            const unsigned int offset = data_faces.data_index_offsets[face];
+
+            for (unsigned int q=0; q<n_q_points_work; ++q)
+              {
+                data_faces.normals_times_jacobians[0][offset+q] =
+                  data_faces.normal_vectors[offset+q] *
+                  data_faces.jacobians[0][offset+q];
+                if (is_boundary_face == false)
+                  data_faces.normals_times_jacobians[1][offset+q] =
+                    data_faces.normal_vectors[offset+q] *
+                    data_faces.jacobians[1][offset+q];
+              }
+          }
+      }
+
+    } // end of anonymous namespace
+
+
+
+    template <int dim, typename Number>
+    void
+    MappingInfo<dim,Number>::initialize_faces
+    (const dealii::Triangulation<dim>                         &tria,
+     const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+     const std::vector<FaceToCellTopology<VectorizedArray<Number>::n_array_elements> > &faces,
+     const Mapping<dim>                                       &mapping,
+     const std::vector<dealii::hp::QCollection<1> >           &quad,
+     const UpdateFlags                                         update_flags_boundary_faces,
+     const UpdateFlags                                         update_flags_inner_faces)
+    {
+      face_type.resize(faces.size(), general);
+      face_data.resize (quad.size());
+
+      // We currently always set the same flags on both inner and boundary
+      // faces to the same value. At some point, we might want to separate the
+      // two.
+      UpdateFlags update_flags_compute_boundary =
+        ((update_flags_inner_faces | update_flags_boundary_faces) & update_quadrature_points ?
+         update_quadrature_points :
+         update_default)
+        | update_normal_vectors | update_JxW_values | update_jacobians;
+      UpdateFlags update_flags_compute_inner =
+        ((update_flags_inner_faces | update_flags_boundary_faces) & update_quadrature_points ?
+         update_quadrature_points :
+         update_default)
+        | update_normal_vectors | update_JxW_values | update_jacobians;
+      UpdateFlags update_flags_common = update_flags_inner_faces | update_flags_boundary_faces;
+
+      for (unsigned int my_q=0; my_q<quad.size(); ++my_q)
+        {
+          face_data[my_q].descriptor.resize(quad[my_q].size());
+          for (unsigned int hpq=0; hpq<quad[my_q].size(); ++hpq)
+            face_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
+                                                       update_flags_compute_inner);
+        }
+
+      if (faces.size() == 0)
+        return;
+
+      // Create as many chunks of cells as we have threads and spawn the work
+      unsigned int work_per_chunk =
+        std::max(std::size_t(8), (faces.size() + MultithreadInfo::n_threads() - 1) /
+                 MultithreadInfo::n_threads());
+
+      std::vector<std::pair<std::vector<MappingInfoStorage<dim-1,dim,Number> >,
+          CompressedFaceData<dim,Number> > > data_faces_local;
+      // Reserve enough space to avoid re-allocation (which would destroy the
+      // references passed to the tasks!)
+      data_faces_local.reserve(MultithreadInfo::n_threads());
+
+      {
+        Threads::TaskGroup<> tasks;
+        std::pair<unsigned int,unsigned int> face_range(0U, work_per_chunk);
+        while (face_range.first < faces.size())
+          {
+            data_faces_local.push_back
+            (std::make_pair (std::vector<MappingInfoStorage<dim-1,dim,Number> >(quad.size()),
+                             CompressedFaceData<dim,Number>(get_jacobian_size(tria))));
+            tasks += Threads::new_task(&initialize_face_range<dim,Number>,
+                                       face_range, tria, cells, faces, mapping,
+                                       update_flags_compute_boundary,
+                                       update_flags_compute_inner,
+                                       *this, data_faces_local.back());
+            face_range.first = face_range.second;
+            face_range.second += work_per_chunk;
+          }
+
+        tasks.join_all();
+      }
+
+      // Fill in each thread's constant data (normals, JxW, normals x
+      // Jacobian) into the data of the zeroth chunk in serial
+      std::vector<std::vector<unsigned int> > indices_compressed(data_faces_local.size());
+      for (unsigned int i=0; i<data_faces_local.size(); ++i)
+        merge_compressed_data(data_faces_local[i].second.data,
+                              data_faces_local[0].second.data,
+                              indices_compressed[i]);
+
+      // Collect all data in the final data fields.
+      // First allocate the memory
+      const unsigned int n_constant_data = data_faces_local[0].second.data.size();
+      for (unsigned int my_q=0; my_q<face_data.size(); ++my_q)
+        {
+          face_data[my_q].data_index_offsets.resize(face_type.size());
+          std::vector<std::array<std::size_t,2> > shift(data_faces_local.size());
+          shift[0][0] = n_constant_data;
+          shift[0][1] = 0;
+          for (unsigned int i=1; i<data_faces_local.size(); ++i)
+            {
+              shift[i][0] = shift[i-1][0] + data_faces_local[i-1].first[my_q].JxW_values.size();
+              shift[i][1] = shift[i-1][1] + data_faces_local[i-1].first[my_q].quadrature_points.size();
+            }
+          face_data[my_q].JxW_values.
+          resize_fast(shift.back()[0] + data_faces_local.back().first[my_q].
+                      JxW_values.size());
+          face_data[my_q].normal_vectors.resize_fast(face_data[my_q].JxW_values.size());
+          face_data[my_q].jacobians[0].resize_fast(face_data[my_q].JxW_values.size());
+          face_data[my_q].jacobians[1].resize_fast(face_data[my_q].JxW_values.size());
+          if (update_flags_common & update_jacobian_grads)
+            {
+              face_data[my_q].jacobian_gradients[0].resize_fast(face_data[my_q].JxW_values.size());
+              face_data[my_q].jacobian_gradients[1].resize_fast(face_data[my_q].JxW_values.size());
+            }
+          face_data[my_q].normals_times_jacobians[0].resize_fast(face_data[my_q].JxW_values.size());
+          face_data[my_q].normals_times_jacobians[1].resize_fast(face_data[my_q].JxW_values.size());
+          if (update_flags_common & update_quadrature_points)
+            {
+              face_data[my_q].quadrature_point_offsets.resize(face_type.size());
+              face_data[my_q].quadrature_points.
+              resize_fast(shift.back()[1] + data_faces_local.back().first[my_q].
+                          quadrature_points.size());
+            }
+
+          // start the tasks to gather the data in parallel
+          Threads::TaskGroup<> tasks;
+          for (unsigned int i=0; i<data_faces_local.size(); ++i)
+            tasks += Threads::new_task(&copy_data<dim-1,dim,Number>,
+                                       work_per_chunk * i, shift[i],
+                                       indices_compressed[i], face_type,
+                                       data_faces_local[i].first[my_q],
+                                       face_data[my_q]);
+
+          // fill the constant data fields (in parallel to the loop above)
+          if (my_q == 0)
+            {
+              const Number jac_size = get_jacobian_size(tria);
+              for (auto &it : data_faces_local[0].second.data)
+                {
+                  // JxW values; invert previously applied scaling
+                  for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                    face_data[my_q].JxW_values[it.second][v]
+                      = it.first[2*dim*dim+dim][v]*
+                        Utilities::fixed_power<dim>(jac_size);
+
+                  // inverse Jacobians
+                  for (unsigned int i=0; i<2; ++i)
+                    for (unsigned int d=0; d<dim; ++d)
+                      for (unsigned int e=0; e<dim; ++e)
+                        for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                          face_data[my_q].jacobians[i][it.second][d][e][v]
+                            = it.first[i*dim*dim+d*dim+e][v];
+
+                  // normal vectors; invert previously applied scaling
+                  for (unsigned int d=0; d<dim; ++d)
+                    for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+                      face_data[my_q].normal_vectors[it.second][d][v]
+                        = it.first[2*dim*dim+d][v] * jac_size;
+                }
+            }
+          else
+            {
+              for (unsigned int i=0; i<n_constant_data; ++i)
+                {
+                  face_data[my_q].JxW_values[i] = face_data[0].JxW_values[i];
+                  face_data[my_q].normal_vectors[i] = face_data[0].normal_vectors[i];
+                  for (unsigned k=0; k<2; ++k)
+                    face_data[my_q].jacobians[k][i] = face_data[0].jacobians[k][i];
+                }
+            }
+
+          // ... wait for the parallel work to finish
+          tasks.join_all();
+
+          // finally compute the normal times the jacobian
+          for (unsigned int i=0; i<data_faces_local.size(); ++i)
+            tasks += Threads::new_task(&compute_normal_times_jacobian<dim,Number>,
+                                       work_per_chunk * i,
+                                       std::min(work_per_chunk * (i+1),
+                                                (unsigned int)faces.size()),
+                                       face_type, faces, face_data[my_q]);
+          tasks.join_all();
+        }
+    }
+
+
+
+    template <int dim, typename Number>
+    void
+    MappingInfo<dim,Number>::initialize_faces_by_cells
+    (const dealii::Triangulation<dim>                         &tria,
+     const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+     const Mapping<dim>                                       &mapping,
+     const std::vector<dealii::hp::QCollection<1> >           &quad,
+     const UpdateFlags                                         update_flags_faces_by_cells)
+    {
+      if (update_flags_faces_by_cells == update_default)
+        return;
+
+      face_data_by_cells.resize(quad.size());
+      const unsigned int n_quads = quad.size();
+      const unsigned int vectorization_width =
+        VectorizedArray<Number>::n_array_elements;
+      UpdateFlags update_flags =
+        (update_flags_faces_by_cells & update_quadrature_points ?
+         update_quadrature_points :
+         update_default)
+        | update_normal_vectors | update_JxW_values | update_jacobians;
+
+      for (unsigned int my_q=0; my_q<n_quads; ++my_q)
+        {
+          const unsigned int n_hp_quads = quad[my_q].size();
+          AssertIndexRange (0, n_hp_quads);
+          face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
+          for (unsigned int q=0; q<n_hp_quads; ++q)
+            face_data_by_cells[my_q].descriptor[q].initialize(quad[my_q][q],
+                                                              update_default);
+
+          // since we already know the cell type, we can pre-allocate the right
+          // amount of data straight away and we just need to do some basic counting
+          AssertDimension(cell_type.size(), cells.size()/vectorization_width);
+          face_data_by_cells[my_q].data_index_offsets.
+          resize(cell_type.size()*GeometryInfo<dim>::faces_per_cell);
+          if (update_flags & update_quadrature_points)
+            face_data_by_cells[my_q].quadrature_point_offsets.
+            resize(cell_type.size()*GeometryInfo<dim>::faces_per_cell);
+          std::size_t storage_length = 0;
+          for (unsigned int i=0; i<cell_type.size(); ++i)
+            for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+              {
+                if (cell_type[i] <= affine)
+                  {
+                    face_data_by_cells[my_q].data_index_offsets
+                    [i*GeometryInfo<dim>::faces_per_cell+face] = storage_length;
+                    ++storage_length;
+                  }
+                else
+                  {
+                    face_data_by_cells[my_q].data_index_offsets
+                    [i*GeometryInfo<dim>::faces_per_cell+face] = storage_length;
+                    storage_length += face_data_by_cells[my_q].descriptor[0].n_q_points;
+                  }
+                if (update_flags & update_quadrature_points)
+                  face_data_by_cells[my_q].quadrature_point_offsets
+                  [i*GeometryInfo<dim>::faces_per_cell+face] =
+                    (i*GeometryInfo<dim>::faces_per_cell+face)*
+                    face_data_by_cells[my_q].descriptor[0].n_q_points;
+              }
+          face_data_by_cells[my_q].JxW_values.resize_fast(storage_length *
+                                                          GeometryInfo<dim>::faces_per_cell);
+          face_data_by_cells[my_q].jacobians[0].resize_fast(storage_length *
+                                                            GeometryInfo<dim>::faces_per_cell);
+          if (update_flags & update_normal_vectors)
+            face_data_by_cells[my_q].normal_vectors.resize_fast(storage_length *
+                                                                GeometryInfo<dim>::faces_per_cell);
+          if (update_flags & update_normal_vectors &&
+              update_flags & update_jacobians)
+            face_data_by_cells[my_q].normals_times_jacobians[0].resize_fast(storage_length *
+                GeometryInfo<dim>::faces_per_cell);
+          if (update_flags & update_jacobian_grads)
+            face_data_by_cells[my_q].jacobian_gradients[0].resize_fast(storage_length *
+                                                                       GeometryInfo<dim>::faces_per_cell);
+
+          if (update_flags & update_quadrature_points)
+            face_data_by_cells[my_q].quadrature_points.
+            resize_fast(cell_type.size()*GeometryInfo<dim>::faces_per_cell*
+                        face_data_by_cells[my_q].descriptor[0].n_q_points);
+        }
+
+      FE_Nothing<dim> dummy_fe;
+      // currently no hp-indices implemented
+      const unsigned int fe_index = 0;
+      std::vector<std::vector<std::shared_ptr<dealii::FEFaceValues<dim> > > >
+      fe_face_values (face_data_by_cells.size());
+      for (unsigned int i=0; i<fe_face_values.size(); ++i)
+        fe_face_values[i].resize(face_data_by_cells[i].descriptor.size());
+      for (unsigned int cell=0; cell<cell_type.size(); ++cell)
+        for (unsigned int my_q=0; my_q<face_data_by_cells.size(); ++my_q)
+          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+            {
+              if (fe_face_values[my_q][fe_index].get() == nullptr)
+                fe_face_values[my_q][fe_index].reset
+                (new dealii::FEFaceValues<dim> (mapping, dummy_fe, face_data_by_cells[my_q].
+                                                descriptor[fe_index].quadrature,
+                                                update_flags));
+              dealii::FEFaceValues<dim> &fe_val = *fe_face_values[my_q][fe_index];
+              const unsigned int offset =
+                face_data_by_cells[my_q].data_index_offsets
+                [cell*GeometryInfo<dim>::faces_per_cell+face];
+
+              for (unsigned int v=0; v<vectorization_width; ++v)
+                {
+                  typename dealii::Triangulation<dim>::cell_iterator
+                  cell_it (&tria, cells[cell*vectorization_width+v].first,
+                           cells[cell*vectorization_width+v].second);
+                  fe_val.reinit(cell_it, face);
+
+                  // copy data for affine data type
+                  if (cell_type[cell] <= affine)
+                    {
+                      if (update_flags & update_JxW_values)
+                        face_data_by_cells[my_q].JxW_values[offset][v]
+                          = fe_val.JxW(0) / face_data_by_cells[my_q].
+                            descriptor[fe_index].quadrature.weight(0);
+                      if (update_flags & update_jacobians)
+                        {
+                          DerivativeForm<1,dim,dim> inv_jac = fe_val.jacobian(0).covariant_form();
+                          for (unsigned int d=0; d<dim; ++d)
+                            for (unsigned int e=0; e<dim; ++e)
+                              {
+                                const unsigned int ee=
+                                  reorder_face_derivative_indices<dim>(face,e);
+                                face_data_by_cells[my_q].jacobians[0][offset][d][e][v]
+                                  = inv_jac[d][ee];
+                              }
+                        }
+                      if (update_flags & update_jacobian_grads)
+                        {
+                          Assert(false, ExcNotImplemented());
+                        }
+                      if (update_flags & update_normal_vectors)
+                        for (unsigned int d=0; d<dim; ++d)
+                          face_data_by_cells[my_q].normal_vectors[offset][d][v]
+                            = fe_val.normal_vector(0)[d];
+                    }
+                  // copy data for general data type
+                  else
+                    {
+                      if (update_flags & update_JxW_values)
+                        for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+                          face_data_by_cells[my_q].JxW_values[offset+q][v]
+                            = fe_val.JxW(q);
+                      if (update_flags & update_jacobians)
+                        for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+                          {
+                            DerivativeForm<1,dim,dim> inv_jac = fe_val.jacobian(q).covariant_form();
+                            for (unsigned int d=0; d<dim; ++d)
+                              for (unsigned int e=0; e<dim; ++e)
+                                {
+                                  const unsigned int ee=
+                                    reorder_face_derivative_indices<dim>(face,e);
+                                  face_data_by_cells[my_q].jacobians[0][offset+q][d][e][v]
+                                    = inv_jac[d][ee];
+                                }
+                          }
+                      if (update_flags & update_jacobian_grads)
+                        {
+                          Assert(false, ExcNotImplemented());
+                        }
+                      if (update_flags & update_normal_vectors)
+                        for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+                          for (unsigned int d=0; d<dim; ++d)
+                            face_data_by_cells[my_q].normal_vectors[offset+q][d][v]
+                              = fe_val.normal_vector(q)[d];
+                    }
+                  if (update_flags & update_quadrature_points)
+                    for (unsigned int q=0; q<fe_val.n_quadrature_points; ++q)
+                      for (unsigned int d=0; d<dim; ++d)
+                        face_data_by_cells[my_q].quadrature_points
+                        [face_data_by_cells[my_q].quadrature_point_offsets
+                         [cell*GeometryInfo<dim>::faces_per_cell+face]+q][d][v] = fe_val.quadrature_point(q)[d];
+                }
+              if (update_flags & update_normal_vectors &&
+                  update_flags & update_jacobians)
+                for (unsigned int q=0; q<(cell_type[cell]<=affine ? 1 :
+                                          fe_val.n_quadrature_points); ++q)
+                  face_data_by_cells[my_q].normals_times_jacobians[0][offset+q]
+                    = face_data_by_cells[my_q].normal_vectors[offset+q] *
+                      face_data_by_cells[my_q].jacobians[0][offset+q];
+            }
+    }
+
+
+
     template <int dim, typename Number>
     std::size_t MappingInfo<dim,Number>::memory_consumption() const
     {
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 390f0143c7..475386de79 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -31,6 +31,7 @@
 #include <deal.II/matrix_free/shape_info.templates.h>
 #include <deal.II/matrix_free/mapping_info.templates.h>
 #include <deal.II/matrix_free/dof_info.templates.h>
+#include <deal.II/matrix_free/face_info.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -188,9 +189,13 @@ internal_reinit(const Mapping<dim>                          &mapping,
   // general case?
   if (additional_data.initialize_mapping == true)
     {
+      std::vector<unsigned int> dummy;
       mapping_info.initialize (dof_handler[0]->get_triangulation(), cell_level_index,
-                               dof_info[0].cell_active_fe_index, mapping, quad,
-                               additional_data.mapping_update_flags);
+                               internal::MatrixFreeFunctions:: FaceInfo
+                               <VectorizedArray<Number>::n_array_elements>(),
+                               dummy, mapping,
+                               quad, additional_data.mapping_update_flags,
+                               update_default, update_default, update_default);
 
       mapping_is_initialized = true;
     }
@@ -316,8 +321,11 @@ internal_reinit(const Mapping<dim>                            &mapping,
   if (additional_data.initialize_mapping == true)
     {
       mapping_info.initialize (dof_handler[0]->get_triangulation(), cell_level_index,
-                               dof_info[0].cell_active_fe_index, mapping, quad,
-                               additional_data.mapping_update_flags);
+                               internal::MatrixFreeFunctions::FaceInfo
+                               <VectorizedArray<Number>::n_array_elements>(),
+                               dof_info[0].cell_active_fe_index, mapping,
+                               quad, additional_data.mapping_update_flags,
+                               update_default, update_default, update_default);
 
       mapping_is_initialized = true;
     }
-- 
2.39.5