Generalize ShapeInfo::reinit

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index 23d07de2bbd639b3ba95be07f55e6ae999ce18fa..a69022a7aa0767faeb1b973817c8b2c50f524584 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -420,6 +420,12 @@ namespace internal
        */
       unsigned int n_q_points_face;
 
+      /**
+       * Stores the number of quadrature points of a face in @p dim dimensions
+       * for simplex, wedge and pyramid reference cells.
+       */
+      std::vector<unsigned int> n_q_points_faces;
+
       /**
        * Stores the number of DoFs per face in @p dim dimensions.
        */

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 1ec957421eaad8764fcca8db11419996a23147c6..4e77bb7e333a631d3a70c85a61b3a60c713e055a 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -314,64 +314,97 @@ namespace internal
                                   q] = grad[d];
               }
 
-          try
-            {
-              const dealii::ReferenceCell reference_cell =
-                dealii::ReferenceCells::get_simplex<dim>();
+          {
+            const auto reference_cell_type = fe.reference_cell();
+
+            const auto  temp      = get_face_quadrature_collection(quad, false);
+            const auto &quad_face = temp.second;
+
+            if (reference_cell_type != temp.first)
+              {
+                // TODO: this might happen if the quadrature rule and the
+                // the FE do not match
+                this->n_q_points_face = 0;
+              }
+            else
+              {
+                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 auto quad_face  = get_face_quadrature(quad);
-              this->n_q_points_face = quad_face.size();
+                unsigned int n_q_points_face_max = 0;
 
-              const unsigned int n_face_orientations = dim == 2 ? 2 : 6;
-              const unsigned int n_faces =
-                dealii::internal::ReferenceCell::get_cell(reference_cell)
-                  .n_faces();
+                for (unsigned int i = 0; i < quad_face.size(); ++i)
+                  n_q_points_face_max =
+                    std::max(n_q_points_face_max, quad_face[i].size());
 
-              const auto projected_quad_face =
-                QProjector<dim>::project_to_all_faces(reference_cell,
-                                                      quad_face);
+                unsigned int n_max_vertices = 0;
 
-              shape_values_face.reinit(
-                {n_faces, n_face_orientations, n_dofs * n_q_points_face});
+                for (unsigned int face_no = 0; face_no < quad_face.size();
+                     ++face_no)
+                  n_max_vertices = std::max(n_max_vertices,
+                                            internal::ReferenceCell::get_face(
+                                              reference_cell_type, face_no)
+                                              .n_vertices());
 
-              shape_gradients_face.reinit(
-                {n_faces, n_face_orientations, dim, n_dofs * n_q_points_face});
+                const auto projected_quad_face =
+                  QProjector<dim>::project_to_all_faces(reference_cell_type,
+                                                        quad_face);
 
-              for (unsigned int f = 0; f < n_faces; ++f)
-                for (unsigned int o = 0; o < n_face_orientations; ++o)
+                const unsigned int n_max_face_orientations =
+                  dim == 2 ? 2 : (2 * n_max_vertices);
+
+                shape_values_face.reinit({quad_face.size(),
+                                          n_max_face_orientations,
+                                          n_dofs * n_q_points_face_max});
+
+                shape_gradients_face.reinit({quad_face.size(),
+                                             n_max_face_orientations,
+                                             dim,
+                                             n_dofs * n_q_points_face_max});
+
+                for (unsigned int f = 0; f < quad_face.size(); ++f)
                   {
-                    const auto offset =
-                      QProjector<dim>::DataSetDescriptor::face(
-                        reference_cell,
-                        f,
-                        (o ^ 1) & 1,  // face_orientation
-                        (o >> 1) & 1, // face_flip
-                        (o >> 2) & 1, // face_rotation
-                        n_q_points_face);
-
-                    for (unsigned int i = 0; i < n_dofs; ++i)
-                      for (unsigned int q = 0; q < n_q_points_face; ++q)
-                        {
-                          const auto point =
-                            projected_quad_face.point(q + offset);
-
-                          shape_values_face(f, o, i * n_q_points_face + q) =
-                            fe.shape_value(i, point);
-
-                          const auto grad = fe.shape_grad(i, point);
-
-                          for (int d = 0; d < dim; ++d)
-                            shape_gradients_face(
-                              f, o, d, i * n_q_points_face + q) = grad[d];
-                        }
+                    const unsigned int n_face_orientations =
+                      dim == 2 ? 2 :
+                                 (2 * internal::ReferenceCell::get_face(
+                                        reference_cell_type, f)
+                                        .n_vertices());
+
+                    const unsigned int n_q_points_face = quad_face[f].size();
+
+                    for (unsigned int o = 0; o < n_face_orientations; ++o)
+                      {
+                        const auto offset =
+                          QProjector<dim>::DataSetDescriptor::face(
+                            reference_cell_type,
+                            f,
+                            (o ^ 1) & 1,  // face_orientation
+                            (o >> 1) & 1, // face_flip
+                            (o >> 2) & 1, // face_rotation
+                            quad_face);
+
+                        for (unsigned int i = 0; i < n_dofs; ++i)
+                          for (unsigned int q = 0; q < n_q_points_face; ++q)
+                            {
+                              const auto point =
+                                projected_quad_face.point(q + offset);
+
+                              shape_values_face(f, o, i * n_q_points_face + q) =
+                                fe.shape_value(i, point);
+
+                              const auto grad = fe.shape_grad(i, point);
+
+                              for (int d = 0; d < dim; ++d)
+                                shape_gradients_face(
+                                  f, o, d, i * n_q_points_face + q) = grad[d];
+                            }
+                      }
                   }
-            }
-          catch (...)
-            {
-              // TODO: this might happen if the quadrature rule and the
-              // the FE do not match
-              this->n_q_points_face = 0;
-            }
+              }
+          }
 
           // TODO: also fill shape_hessians, inverse_shape_values,
           //   shape_data_on_face, quadrature_data_on_face,

diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
index b46f3c27f5da87aab6227a1f26dd81c84890b994..495c6a6e6dc7840c556ec697f6e0bb6ab62f33b8 100644 (file)
--- a/include/deal.II/matrix_free/util.h
+++ b/include/deal.II/matrix_free/util.h
@@ -22,6 +22,10 @@
 
 #include <deal.II/base/quadrature.h>
 
+#include <deal.II/grid/reference_cell.h>
+
+#include <deal.II/hp/q_collection.h>
+
 #include <deal.II/simplex/quadrature_lib.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -45,6 +49,70 @@ namespace internal
       return Quadrature<dim - 1>();
     }
 
+    template <int dim>
+    inline std::pair<dealii::ReferenceCell, dealii::hp::QCollection<dim - 1>>
+    get_face_quadrature_collection(const Quadrature<dim> &quad,
+                                   const bool             do_assert = true)
+    {
+      if (dim == 2 || dim == 3)
+        {
+          for (unsigned int i = 1; i <= 4; ++i)
+            if (quad == Simplex::QGauss<dim>(i))
+              {
+                Simplex::QGauss<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)};
+              }
+
+          for (unsigned int i = 1; i <= 5; ++i)
+            if (quad == Simplex::QWitherdenVincent<dim>(i))
+              {
+                Simplex::QWitherdenVincent<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)};
+              }
+        }
+
+      if (dim == 3)
+        for (unsigned int i = 1; i <= 3; ++i)
+          if (quad == Simplex::QGaussWedge<dim>(i))
+            {
+              QGauss<dim - 1>          quad(i);
+              Simplex::QGauss<dim - 1> tri(i);
+
+              return {
+                ReferenceCells::Wedge,
+                dealii::hp::QCollection<dim - 1>(tri, tri, quad, quad, quad)};
+            }
+
+      if (dim == 3)
+        for (unsigned int i = 1; i <= 2; ++i)
+          if (quad == Simplex::QGaussPyramid<dim>(i))
+            {
+              QGauss<dim - 1>          quad(i);
+              Simplex::QGauss<dim - 1> tri(i);
+
+              return {
+                ReferenceCells::Pyramid,
+                dealii::hp::QCollection<dim - 1>(quad, tri, tri, tri, tri)};
+            }
+
+      if (do_assert)
+        AssertThrow(false, ExcNotImplemented());
+
+      return {ReferenceCells::Invalid, dealii::hp::QCollection<dim - 1>()};
+    }
+
   } // end of namespace MatrixFreeFunctions
 } // end of namespace internal