Use QCollection within MapptingInfo::initialize

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index 8c0d0c77d46e7badb001ecbca42b1b036788f719..4ba8c6b73c202ca3d233c66a5b50c7a2bcba0ab9 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -27,6 +27,8 @@
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/mapping.h>
 
+#include <deal.II/grid/reference_cell.h>
+
 #include <deal.II/hp/mapping_collection.h>
 #include <deal.II/hp/q_collection.h>
 
@@ -194,6 +196,14 @@ namespace internal
        */
       std::vector<QuadratureDescriptor> descriptor;
 
+      /**
+       * Collection of quadrature formulae applied on the given face.
+       *
+       * @note Only filled for faces, since faces might be quadrilateral or
+       *   triangle shaped.
+       */
+      std::vector<dealii::hp::QCollection<structdim>> q_collection;
+
       /**
        * Stores the index offset into the arrays @p jxw_values, @p jacobians,
        * @p normal_vectors and the second derivatives. Note that affine cells
@@ -454,6 +464,12 @@ namespace internal
        */
       SmartPointer<const Mapping<dim>> mapping;
 
+      /**
+       * Reference-cell type related to each quadrature and active quadrature
+       * index.
+       */
+      std::vector<std::vector<dealii::ReferenceCell>> reference_cell_types;
+
       /**
        * Internal function to compute the geometry for the case the mapping is
        * a MappingQ and a single quadrature formula per slot (non-hp-case) is

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 445e28dea7f26704c99ba788cac0f764308b469e..d9d2f84522f64804cc36ddcc142c121efed47373 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -362,37 +362,21 @@ namespace internal
       this->update_flags_inner_faces    = this->update_flags_boundary_faces;
       this->update_flags_faces_by_cells = update_flags_faces_by_cells;
 
+      reference_cell_types.resize(quad.size());
+
       for (unsigned int my_q = 0; my_q < quad.size(); ++my_q)
         {
           const unsigned int n_hp_quads = quad[my_q].size();
           AssertIndexRange(0, n_hp_quads);
-          cell_data[my_q].descriptor.resize(n_hp_quads);
-          for (unsigned int q = 0; q < n_hp_quads; ++q)
-            {
-              bool flag = quad[my_q][q].is_tensor_product();
 
-              if (flag)
-                for (unsigned int i = 1; i < dim; ++i)
-                  {
-                    flag &= quad[my_q][q].get_tensor_basis()[0] ==
-                            quad[my_q][q].get_tensor_basis()[i];
-                  }
+          const unsigned int scale = std::max<unsigned int>(1, dim - 1);
 
-              if (flag == false)
-                {
-#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
-                                                           update_default);
-#else
-                  Assert(false, ExcNotImplemented());
-#endif
-                }
-              else
-                cell_data[my_q].descriptor[q].initialize(
-                  quad[my_q][q].get_tensor_basis()[0], update_default);
-            }
+          cell_data[my_q].descriptor.resize(n_hp_quads);
+          face_data[my_q].descriptor.resize(n_hp_quads * scale);
+          face_data[my_q].q_collection.resize(n_hp_quads);
+          face_data_by_cells[my_q].descriptor.resize(n_hp_quads * scale);
+          reference_cell_types[my_q].resize(n_hp_quads);
 
-          face_data[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
             {
               bool flag = quad[my_q][hpq].is_tensor_product();
@@ -405,58 +389,57 @@ namespace internal
               if (flag == false)
                 {
 #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  try
+                  cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
+                                                             update_default);
+                  const auto quad_face =
+                    get_unique_face_quadratures(quad[my_q][hpq]);
+
+                  if (quad_face.first.size() > 0) // line, quad
                     {
-                      const auto quad_face =
-                        get_face_quadrature(quad[my_q][hpq]);
-                      face_data[my_q].descriptor[hpq].initialize(
-                        quad_face, update_default);
+                      face_data[my_q].descriptor[hpq * scale].initialize(
+                        quad_face.first, update_default);
+                      face_data_by_cells[my_q]
+                        .descriptor[hpq * scale]
+                        .initialize(quad_face.first, update_default);
                     }
-                  catch (...)
+
+                  if (quad_face.second.size() > 0) // triangle
                     {
-                      // TODO: nothing to do for now for wedges and pyramids.
+                      AssertDimension(dim, 3);
+                      face_data[my_q].descriptor[hpq * scale + 1].initialize(
+                        quad_face.second, update_default);
+                      face_data_by_cells[my_q]
+                        .descriptor[hpq * scale + 1]
+                        .initialize(quad_face.second, update_default);
                     }
-#else
-                  Assert(false, ExcNotImplemented());
-#endif
-                }
-              else
-                face_data[my_q].descriptor[hpq].initialize(
-                  quad[my_q][hpq].get_tensor_basis()[0],
-                  update_flags_boundary_faces);
-            }
 
-          face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
-          for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
-            {
-              bool flag = quad[my_q][hpq].is_tensor_product();
+                  const auto face_quadrature_collection =
+                    get_face_quadrature_collection(quad[my_q][hpq]);
 
-              if (flag)
-                for (unsigned int i = 1; i < dim; ++i)
-                  flag &= quad[my_q][hpq].get_tensor_basis()[0] ==
-                          quad[my_q][hpq].get_tensor_basis()[i];
+                  face_data[my_q].q_collection[hpq] =
+                    face_quadrature_collection.second;
 
-              if (flag == false)
-                {
-#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  try
-                    {
-                      const auto quad_face =
-                        get_face_quadrature(quad[my_q][hpq]);
-                      face_data_by_cells[my_q].descriptor[hpq].initialize(
-                        quad_face, update_default);
-                    }
-                  catch (...)
-                    {
-                      // TODO: nothing to do for now for wedges and pyramids.
-                    }
+                  reference_cell_types[my_q][hpq] =
+                    face_quadrature_collection.first;
 #else
                   Assert(false, ExcNotImplemented());
 #endif
                 }
               else
-                face_data_by_cells[my_q].descriptor[hpq].initialize(
-                  quad[my_q][hpq].get_tensor_basis()[0], update_default);
+                {
+                  cell_data[my_q].descriptor[hpq].initialize(
+                    quad[my_q][hpq].get_tensor_basis()[0], update_default);
+
+                  const auto quad_face = quad[my_q][hpq].get_tensor_basis()[0];
+                  face_data[my_q].descriptor[hpq * scale].initialize(
+                    quad_face, update_flags_boundary_faces);
+                  face_data[my_q].q_collection[hpq] =
+                    dealii::hp::QCollection<dim - 1>(quad_face);
+                  face_data_by_cells[my_q].descriptor[hpq * scale].initialize(
+                    quad_face, update_default);
+                  reference_cell_types[my_q][hpq] =
+                    dealii::ReferenceCells::get_hypercube<dim>();
+                }
             }
         }
 
@@ -1804,7 +1787,20 @@ namespace internal
             MappingInfoStorage<dim - 1, dim, Number, VectorizedArrayType>>,
           CompressedFaceData<dim, Number, VectorizedArrayType>> &data)
       {
-        FE_Nothing<dim> dummy_fe;
+        std::vector<std::vector<std::shared_ptr<FE_Nothing<dim>>>> dummy_fe(
+          mapping_info.reference_cell_types.size());
+        for (unsigned int my_q = 0;
+             my_q < mapping_info.reference_cell_types.size();
+             ++my_q)
+          {
+            const unsigned int n_hp_quads =
+              mapping_info.reference_cell_types[my_q].size();
+            dummy_fe[my_q].resize(n_hp_quads);
+
+            for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
+              dummy_fe[my_q][hpq] = std::make_shared<FE_Nothing<dim>>(
+                mapping_info.reference_cell_types[my_q][hpq]);
+          }
 
         const unsigned int max_active_fe_index =
           active_fe_index.size() > 0 ?
@@ -1860,11 +1856,9 @@ namespace internal
               const unsigned int hp_mapping_index =
                 mapping_in.size() == 1 ? 0 : fe_index;
 
-              const auto &               mapping = mapping_in[hp_mapping_index];
-              const Quadrature<dim - 1> &quadrature =
-                mapping_info.face_data[my_q]
-                  .descriptor[hp_quad_index]
-                  .quadrature;
+              const auto &mapping = mapping_in[hp_mapping_index];
+              const auto &quadrature =
+                mapping_info.face_data[my_q].q_collection[hp_quad_index];
 
               const bool is_boundary_face =
                 faces[face].cells_exterior[0] == numbers::invalid_unsigned_int;
@@ -1874,14 +1868,14 @@ namespace internal
                 fe_boundary_face_values_container[my_q][fe_index] =
                   std::make_shared<FEFaceValues<dim>>(
                     mapping,
-                    dummy_fe,
+                    *dummy_fe[my_q][hp_quad_index],
                     quadrature,
                     mapping_info.update_flags_boundary_faces);
               else if (fe_face_values_container[my_q][fe_index] == nullptr)
                 fe_face_values_container[my_q][fe_index] =
                   std::make_shared<FEFaceValues<dim>>(
                     mapping,
-                    dummy_fe,
+                    *dummy_fe[my_q][hp_quad_index],
                     quadrature,
                     mapping_info.update_flags_inner_faces);
 
@@ -1889,9 +1883,9 @@ namespace internal
                 is_boundary_face ?
                   *fe_boundary_face_values_container[my_q][fe_index] :
                   *fe_face_values_container[my_q][fe_index];
-              const unsigned int n_q_points =
-                fe_face_values.n_quadrature_points;
-              face_data.resize(n_q_points);
+
+              unsigned int n_q_points = 0; // will be override once FEFaceValues
+                                           // is set up
 
               bool normal_is_similar = true;
               bool JxW_is_similar    = true;
@@ -1911,13 +1905,22 @@ namespace internal
                       fe_face_values.reinit(cell_it,
                                             faces[face].interior_face_no);
 
+                      if (v == 0)
+                        {
+                          n_q_points = fe_face_values.n_quadrature_points;
+                          face_data.resize(n_q_points);
+                        }
+
                       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)) >
+                                fe_face_values.JxW(q) *
+                                  fe_face_values.get_quadrature().weight(0) -
+                                fe_face_values.JxW(0) *
+                                  fe_face_values.get_quadrature().weight(q)) >
                               2048. * std::numeric_limits<double>::epsilon() *
-                                fe_face_values.JxW(0) * quadrature.weight(q))
+                                fe_face_values.JxW(0) *
+                                fe_face_values.get_quadrature().weight(q))
                             JxW_is_similar = false;
                           face_data.JxW_values[q][v] = fe_face_values.JxW(q);
 
@@ -2043,17 +2046,16 @@ namespace internal
                             mapping_in.size() == 1 ? 0 : fe_index;
 
                           const auto &mapping = mapping_in[hp_mapping_index];
-                          const Quadrature<dim - 1> &quadrature =
+                          const auto &quadrature =
                             mapping_info.face_data[my_q]
-                              .descriptor[hp_quad_index]
-                              .quadrature;
+                              .q_collection[hp_quad_index];
 
                           if (fe_face_values_container[my_q][fe_index] ==
                               nullptr)
                             fe_face_values_container[my_q][fe_index] =
                               std::make_shared<FEFaceValues<dim>>(
                                 mapping,
-                                dummy_fe,
+                                *dummy_fe[my_q][hp_quad_index],
                                 quadrature,
                                 mapping_info.update_flags_boundary_faces);
 
@@ -2069,7 +2071,7 @@ namespace internal
                             fe_subface_values_container[my_q][0] =
                               std::make_shared<FESubfaceValues<dim>>(
                                 mapping,
-                                dummy_fe,
+                                *dummy_fe[my_q][hp_quad_index],
                                 quadrature,
                                 mapping_info.update_flags_inner_faces);
                           fe_subface_values_container[my_q][0]->reinit(
@@ -2178,7 +2180,8 @@ namespace internal
                       for (unsigned int v = 0; v < VectorizedArrayType::size();
                            ++v)
                         new_entry.first[2 * dim * dim + dim][v] =
-                          face_data.JxW_values[0][v] / quadrature.weight(0) /
+                          face_data.JxW_values[0][v] /
+                          fe_face_values.get_quadrature().weight(0) /
                           Utilities::fixed_power<dim>(face_data.jac_size);
 
                       new_entry.second = data.second.data.size();
@@ -3165,14 +3168,14 @@ namespace internal
                   std::make_shared<dealii::FEFaceValues<dim>>(
                     mapping,
                     dummy_fe,
-                    face_data_by_cells[my_q].descriptor[fe_index].quadrature,
+                    face_data[my_q].q_collection[fe_index],
                     update_flags);
               if (fe_face_values_neigh[my_q][fe_index].get() == nullptr)
                 fe_face_values_neigh[my_q][fe_index] =
                   std::make_shared<dealii::FEFaceValues<dim>>(
                     mapping,
                     dummy_fe,
-                    face_data_by_cells[my_q].descriptor[fe_index].quadrature,
+                    face_data[my_q].q_collection[fe_index],
                     update_flags);
               dealii::FEFaceValues<dim> &fe_val =
                 *fe_face_values[my_q][fe_index];
@@ -3215,9 +3218,7 @@ namespace internal
                     {
                       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);
+                          fe_val.JxW(0) / fe_val.get_quadrature().weight(0);
                       if (update_flags & update_jacobians)
                         {
                           DerivativeForm<1, dim, dim> inv_jac =

diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
index 164b759e275189deb0befa2726c54d77ea3e1271..2c47b96fab756855d3e1110a89ddd0fd9b490661 100644 (file)
--- a/include/deal.II/matrix_free/util.h
+++ b/include/deal.II/matrix_free/util.h
@@ -34,20 +34,6 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
-    template <int dim>
-    inline Quadrature<dim - 1>
-    get_face_quadrature(const Quadrature<dim> &quad)
-    {
-      if (dim == 2 || dim == 3)
-        for (unsigned int i = 1; i <= 3; ++i)
-          if (quad == QGaussSimplex<dim>(i))
-            return QGaussSimplex<dim - 1>(i);
-
-      AssertThrow(false, ExcNotImplemented());
-
-      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,
@@ -112,6 +98,36 @@ namespace internal
       return {ReferenceCells::Invalid, dealii::hp::QCollection<dim - 1>()};
     }
 
+
+
+    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)};
+            }
+
+      if (dim == 3)
+        for (unsigned int i = 1; i <= 3; ++i)
+          if (quad == QGaussWedge<dim>(i))
+            return {QGauss<dim - 1>(i), QGaussSimplex<dim - 1>(i)};
+
+      if (dim == 3)
+        for (unsigned int i = 1; i <= 2; ++i)
+          if (quad == QGaussPyramid<dim>(i))
+            return {QGauss<dim - 1>(i), QGaussSimplex<dim - 1>(i)};
+
+      AssertThrow(false, ExcNotImplemented());
+
+      return {QGauss<dim - 1>(1), QGauss<dim - 1>(1)};
+    }
   } // end of namespace MatrixFreeFunctions
 } // end of namespace internal
 

diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc
index 0cf57e3f635b14cc8523ec854eedf3a51ab30add..131750788cce05a0eda8f433841cbbc60cf1c35c 100644 (file)
--- a/tests/simplex/matrix_free_03.cc
+++ b/tests/simplex/matrix_free_03.cc
@@ -169,8 +169,10 @@ public:
           }
       },
       [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
-        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, true);
-        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval_neighbor(data, false);
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, face_range, true);
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval_neighbor(data,
+                                                                 face_range,
+                                                                 false);
 
         for (unsigned int face = face_range.first; face < face_range.second;
              face++)
@@ -214,7 +216,7 @@ public:
           }
       },
       [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
-        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, true);
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, face_range, true);
         for (unsigned int face = face_range.first; face < face_range.second;
              face++)
           {