Work on get_face_quadrature_collection() and get_unique_face_quadratures()

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 898c7a192b96dfdb01a022bddafc30005057fecd..2a9b360878ceb56b4607e73420a1d0bde843612d 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -328,9 +328,12 @@ namespace internal
               {
                 this->n_q_points_face = quad_face[0].size();
 
-                n_q_points_faces.resize(quad_face.size());
-                for (unsigned int i = 0; i < quad_face.size(); ++i)
-                  n_q_points_faces[i] = quad_face[i].size();
+                const unsigned int n_faces = temp.first.n_faces();
+
+                n_q_points_faces.resize(n_faces);
+                for (unsigned int i = 0; i < n_faces; ++i)
+                  n_q_points_faces[i] =
+                    quad_face[quad_face.size() == 1 ? 0 : i].size();
 
                 unsigned int n_q_points_face_max = 0;
 
@@ -340,8 +343,7 @@ namespace internal
 
                 unsigned int n_max_vertices = 0;
 
-                for (unsigned int face_no = 0; face_no < quad_face.size();
-                     ++face_no)
+                for (unsigned int face_no = 0; face_no < n_faces; ++face_no)
                   n_max_vertices = std::max(
                     n_max_vertices,
                     reference_cell.face_reference_cell(face_no).n_vertices());
@@ -353,16 +355,16 @@ namespace internal
                 const unsigned int n_max_face_orientations =
                   dim == 2 ? 2 : (2 * n_max_vertices);
 
-                shape_values_face.reinit({quad_face.size(),
+                shape_values_face.reinit({n_faces,
                                           n_max_face_orientations,
                                           n_dofs * n_q_points_face_max});
 
-                shape_gradients_face.reinit({quad_face.size(),
+                shape_gradients_face.reinit({n_faces,
                                              n_max_face_orientations,
                                              dim,
                                              n_dofs * n_q_points_face_max});
 
-                for (unsigned int f = 0; f < quad_face.size(); ++f)
+                for (unsigned int f = 0; f < n_faces; ++f)
                   {
                     const unsigned int n_face_orientations =
                       dim == 2 ?
@@ -370,7 +372,8 @@ namespace internal
                         (2 *
                          reference_cell.face_reference_cell(f).n_vertices());
 
-                    const unsigned int n_q_points_face = quad_face[f].size();
+                    const unsigned int n_q_points_face =
+                      quad_face[quad_face.size() == 1 ? 0 : f].size();
 
                     for (unsigned int o = 0; o < n_face_orientations; ++o)
                       {

diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
index 6b43808cfd5a345813100ca7aedda6c8a615af2a..323f2093e0486ddbda5fd2f2a2c6c8e14d2dee0a 100644 (file)
--- a/include/deal.II/matrix_free/util.h
+++ b/include/deal.II/matrix_free/util.h
@@ -34,6 +34,11 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
+    /**
+     * Given a quadrature rule @p quad defined on a cell, return the type of
+     * the cell and a collection of lower-dimensional quadrature rules that
+     * are defined on each face.
+     */
     template <int dim>
     inline std::pair<dealii::ReferenceCell, dealii::hp::QCollection<dim - 1>>
     get_face_quadrature_collection(const Quadrature<dim> &quad,
@@ -43,29 +48,15 @@ namespace internal
         {
           for (unsigned int i = 1; i <= 4; ++i)
             if (quad == QGaussSimplex<dim>(i))
-              {
-                QGaussSimplex<dim - 1> tri(i);
-
-                if (dim == 2)
-                  return {ReferenceCells::Triangle,
-                          dealii::hp::QCollection<dim - 1>(tri, tri, tri)};
-                else
-                  return {ReferenceCells::Tetrahedron,
-                          dealii::hp::QCollection<dim - 1>(tri, tri, tri, tri)};
-              }
+              return {ReferenceCells::get_simplex<dim>(),
+                      dealii::hp::QCollection<dim - 1>(
+                        QGaussSimplex<dim - 1>(i))};
 
           for (unsigned int i = 1; i <= 5; ++i)
             if (quad == QWitherdenVincentSimplex<dim>(i))
-              {
-                QWitherdenVincentSimplex<dim - 1> tri(i);
-
-                if (dim == 2)
-                  return {ReferenceCells::Triangle,
-                          dealii::hp::QCollection<dim - 1>(tri, tri, tri)};
-                else
-                  return {ReferenceCells::Tetrahedron,
-                          dealii::hp::QCollection<dim - 1>(tri, tri, tri, tri)};
-              }
+              return {ReferenceCells::get_simplex<dim>(),
+                      dealii::hp::QCollection<dim - 1>(
+                        QWitherdenVincentSimplex<dim - 1>(i))};
         }
 
       if (dim == 3)
@@ -92,6 +83,13 @@ namespace internal
                 dealii::hp::QCollection<dim - 1>(quad, tri, tri, tri, tri)};
             }
 
+      // note: handle hypercubes last since normally this function is not
+      // called for hypercubes
+      for (unsigned int i = 1; i <= 5; ++i)
+        if (quad == QGauss<dim>(i))
+          return {ReferenceCells::get_hypercube<dim>(),
+                  dealii::hp::QCollection<dim - 1>(QGauss<dim - 1>(i))};
+
       if (do_assert)
         AssertThrow(false, ExcNotImplemented());
 
@@ -100,19 +98,41 @@ namespace internal
 
 
 
+    /**
+     * Return face quadrature rules. In contrast to
+     * get_face_quadrature_collection(), it does not return the quadrature
+     * face for each face but returns one of each type. The first entry
+     * of the returned pair might contain a quadrature rule defined on lines and
+     * quadrilaterals, while the second entry might contain a quadrature rule
+     * defined on a triangle.
+     */
     template <int dim>
     inline std::pair<Quadrature<dim - 1>, Quadrature<dim - 1>>
     get_unique_face_quadratures(const Quadrature<dim> &quad)
     {
       if (dim == 2 || dim == 3)
-        for (unsigned int i = 1; i <= 3; ++i)
-          if (quad == QGaussSimplex<dim>(i))
-            {
-              if (dim == 2)
-                return {QGaussSimplex<dim - 1>(i), Quadrature<dim - 1>()};
-              else
-                return {Quadrature<dim - 1>(), QGaussSimplex<dim - 1>(i)};
-            }
+        {
+          for (unsigned int i = 1; i <= 4; ++i)
+            if (quad == QGaussSimplex<dim>(i))
+              {
+                if (dim == 2)
+                  return {QGaussSimplex<dim - 1>(i), // line!
+                          Quadrature<dim - 1>()};
+                else
+                  return {Quadrature<dim - 1>(), QGaussSimplex<dim - 1>(i)};
+              }
+
+          for (unsigned int i = 1; i <= 5; ++i)
+            if (quad == QWitherdenVincentSimplex<dim>(i))
+              {
+                if (dim == 2)
+                  return {QWitherdenVincentSimplex<dim - 1>(i), // line!
+                          Quadrature<dim - 1>()};
+                else
+                  return {Quadrature<dim - 1>(),
+                          QWitherdenVincentSimplex<dim - 1>(i)};
+              }
+        }
 
       if (dim == 3)
         for (unsigned int i = 1; i <= 3; ++i)
@@ -124,9 +144,15 @@ namespace internal
           if (quad == QGaussPyramid<dim>(i))
             return {QGauss<dim - 1>(i), QGaussSimplex<dim - 1>(i)};
 
+      // note: handle hypercubes last since normally this function is not
+      // called for hypercubes
+      for (unsigned int i = 1; i <= 5; ++i)
+        if (quad == QGauss<dim>(i))
+          return {QGauss<dim - 1>(i), Quadrature<dim - 1>()};
+
       AssertThrow(false, ExcNotImplemented());
 
-      return {QGauss<dim - 1>(1), QGauss<dim - 1>(1)};
+      return {Quadrature<dim - 1>(), Quadrature<dim - 1>()};
     }
   } // end of namespace MatrixFreeFunctions
 } // end of namespace internal