Introduce internal::VectorizedArrayTrait and use it

diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h
index 93432f5a9e1a883941be33a7629c5add80971935..369863610f15d869eb855ba0319af465f02cbc0a 100644 (file)
--- a/include/deal.II/base/vectorization.h
+++ b/include/deal.II/base/vectorization.h
@@ -5317,6 +5317,22 @@ compare_and_apply_mask(const VectorizedArray<double, 2> &left,
 #endif // DOXYGEN
 
 
+namespace internal
+{
+  template <typename T>
+  struct VectorizedArrayTrait
+  {
+    using value_type = T;
+  };
+
+  template <typename T, std::size_t width>
+  struct VectorizedArrayTrait<VectorizedArray<T, width>>
+  {
+    using value_type = T;
+  };
+} // namespace internal
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 /**

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index bab18e73a4988bd1cf47186b9e484ace9576d414..59fe6a9bc16d7bcff29ab406f2f0cc1093bc2675 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -23,6 +23,7 @@
 #include <deal.II/base/aligned_vector.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/vectorization.h>
 
 #include <deal.II/fe/fe.h>
 
@@ -267,6 +268,15 @@ namespace internal
        */
       std::array<AlignedVector<Number>, 2> subface_interpolation_matrices;
 
+      /**
+       * Same as above but stored in a scalar format independent of the type of
+       * Number
+       */
+      std::array<AlignedVector<typename dealii::internal::VectorizedArrayTrait<
+                   Number>::value_type>,
+                 2>
+        subface_interpolation_matrices_scalar;
+
       /**
        * We store a copy of the one-dimensional quadrature formula
        * used for initialization.

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index e1643d51daa26b011be6e526bd6ce0b60213b942..c8a482d752c7b83b36891e2c9e0a07e6d1c3f342 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -631,6 +631,10 @@ namespace internal
             univariate_shape_data.subface_interpolation_matrices[0];
           auto &subface_interpolation_matrix_1 =
             univariate_shape_data.subface_interpolation_matrices[1];
+          auto &subface_interpolation_matrix_scalar_0 =
+            univariate_shape_data.subface_interpolation_matrices_scalar[0];
+          auto &subface_interpolation_matrix_scalar_1 =
+            univariate_shape_data.subface_interpolation_matrices_scalar[1];
 
           const auto fe_1d = create_fe<1>(fe);
           const auto fe_2d = create_fe<2>(fe);
@@ -666,6 +670,11 @@ namespace internal
           subface_interpolation_matrix_1.resize(fe_1d->n_dofs_per_cell() *
                                                 fe_1d->n_dofs_per_cell());
 
+          subface_interpolation_matrix_scalar_0.resize(
+            fe_1d->n_dofs_per_cell() * fe_1d->n_dofs_per_cell());
+          subface_interpolation_matrix_scalar_1.resize(
+            fe_1d->n_dofs_per_cell() * fe_1d->n_dofs_per_cell());
+
           for (unsigned int i = 0, c = 0; i < fe_1d->n_dofs_per_cell(); ++i)
             for (unsigned int j = 0; j < fe_1d->n_dofs_per_cell(); ++j, ++c)
               {
@@ -675,6 +684,13 @@ namespace internal
                 subface_interpolation_matrix_1[c] =
                   interpolation_matrix_1(scalar_lexicographic[i],
                                          scalar_lexicographic[j]);
+
+                subface_interpolation_matrix_scalar_0[c] =
+                  interpolation_matrix_0(scalar_lexicographic[i],
+                                         scalar_lexicographic[j]);
+                subface_interpolation_matrix_scalar_1[c] =
+                  interpolation_matrix_1(scalar_lexicographic[i],
+                                         scalar_lexicographic[j]);
               }
         }