From fae90c1fe9ff4454dbdded312bb5c5149c4aff23 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 17 Nov 2020 19:50:16 +0100
Subject: [PATCH] Actually use ShapeInfo for Simplex::FE_P and Simplex::QGauss

---
 include/deal.II/matrix_free/mapping_info.h    |  8 ++
 .../matrix_free/mapping_info.templates.h      | 88 +++++++++++++++----
 2 files changed, 77 insertions(+), 19 deletions(-)

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index e4852ea9eb..d717e3fc17 100644
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -125,6 +125,14 @@ namespace internal
          */
         QuadratureDescriptor();
 
+        /**
+         * Set up the lengths in the various members of this struct.
+         */
+        template <int dim_q>
+        void
+        initialize(const Quadrature<dim_q> &quadrature,
+                   const UpdateFlags update_flags_inner_faces = update_default);
+
         /**
          * Set up the lengths in the various members of this struct.
          */
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 1d89993cc5..e81fcb0fe8 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -51,6 +51,34 @@ namespace internal
 
 
 
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    template <int dim_q>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      QuadratureDescriptor::initialize(
+        const Quadrature<dim_q> &quadrature,
+        const UpdateFlags        update_flags_inner_faces)
+    {
+      Assert(structdim + 1 <= spacedim ||
+               update_flags_inner_faces == update_default,
+             ExcMessage("Volume cells do not allow for setting inner faces"));
+      this->quadrature = quadrature;
+      n_q_points       = quadrature.size();
+      quadrature_weights.resize(n_q_points);
+      for (unsigned int i = 0; i < n_q_points; ++i)
+        quadrature_weights[i] = quadrature.weight(i);
+
+      // note: quadrature_1 and tensor_quadrature_weights are not set up
+
+      // TODO: set up face_orientations
+      (void)update_flags_inner_faces;
+    }
+
+
+
     template <int structdim,
               int spacedim,
               typename Number,
@@ -338,27 +366,41 @@ namespace internal
           cell_data[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int q = 0; q < n_hp_quads; ++q)
             {
-              Assert(quad[my_q][q].is_tensor_product(), ExcNotImplemented());
-              for (unsigned int i = 1; i < dim; ++i)
-                {
-                  Assert(quad[my_q][q].get_tensor_basis()[0] ==
-                           quad[my_q][q].get_tensor_basis()[i],
-                         ExcNotImplemented());
-                }
+              bool flag = quad[my_q][q].is_tensor_product();
 
-              cell_data[my_q].descriptor[q].initialize(
-                quad[my_q][q].get_tensor_basis()[0], update_default);
+              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];
+                  }
+#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
+              if (flag == false)
+                cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
+                                                         update_default);
+              else
+#endif
+                cell_data[my_q].descriptor[q].initialize(
+                  quad[my_q][q].get_tensor_basis()[0], update_default);
             }
 
           face_data[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
             {
-              Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
-              for (unsigned int i = 1; i < dim; ++i)
+              bool flag = quad[my_q][hpq].is_tensor_product();
+
+              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];
+
+              if (flag == false)
                 {
-                  Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
-                           quad[my_q][hpq].get_tensor_basis()[i],
-                         ExcNotImplemented());
+#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
+                  continue; // TODO: face data is not set up
+#else
+                  Assert(false, ExcNotImplemented());
+#endif
                 }
 
               face_data[my_q].descriptor[hpq].initialize(
@@ -369,12 +411,20 @@ namespace internal
           face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
             {
-              Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
-              for (unsigned int i = 1; i < dim; ++i)
+              bool flag = quad[my_q][hpq].is_tensor_product();
+
+              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];
+
+              if (flag == false)
                 {
-                  Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
-                           quad[my_q][hpq].get_tensor_basis()[i],
-                         ExcNotImplemented());
+#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
+                  continue; // TODO: face data is not set up
+#else
+                  Assert(false, ExcNotImplemented());
+#endif
                 }
 
               face_data_by_cells[my_q].descriptor[hpq].initialize(
-- 
2.39.5