Introduce ShapeInfo::subface_interpolation_matrix

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index 51d61df5cf82279301015afd2ca9fd8ef6491245..ac641152f8d2b66e2e7503bb76763bba9067ac9e 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -261,6 +261,12 @@ namespace internal
        */
       std::array<AlignedVector<Number>, 2> hessians_within_subface;
 
+      /**
+       * A 1D subface interpolation matrix to the first quadrant. This data
+       * structure is only set up for FE_Q for dim > 1.
+       */
+      AlignedVector<Number> subface_interpolation_matrix;
+
       /**
        * 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 2a9b360878ceb56b4607e73420a1d0bde843612d..3c716a4fe48436bb78a792e1f3d6a3b662a2aacc 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -30,9 +30,11 @@
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_q_dg0.h>
 #include <deal.II/fe/fe_simplex_p.h>
 #include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/fe_wedge_p.h>
 
 #include <deal.II/grid/reference_cell.h>
@@ -171,6 +173,29 @@ namespace internal
 
 
 
+    template <int dim_to, int dim, int spacedim>
+    std::unique_ptr<FiniteElement<dim_to, dim_to>>
+    create_fe(const FiniteElement<dim, spacedim> &fe)
+    {
+      std::string fe_name = fe.get_name();
+
+      Assert(
+        fe_name.find("FESystem") == std::string::npos,
+        ExcMessage(
+          "This function can not accept FESystem but only base elements."));
+
+      {
+        const std::size_t template_starts = fe_name.find_first_of('<');
+        Assert(fe_name[template_starts + 1] ==
+                 (dim == 1 ? '1' : (dim == 2 ? '2' : '3')),
+               ExcInternalError());
+        fe_name[template_starts + 1] = std::to_string(dim_to).c_str()[0];
+      }
+      return FETools::get_fe_by_name<dim_to, dim_to>(fe_name);
+    }
+
+
+
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo()
       : element_type(tensor_general)
@@ -600,6 +625,43 @@ namespace internal
             }
         }
 
+      if (dim > 1 && dynamic_cast<const FE_Q<dim> *>(&fe))
+        {
+          auto &subface_interpolation_matrix =
+            univariate_shape_data.subface_interpolation_matrix;
+
+          const auto fe_1d = create_fe<1>(fe);
+          const auto fe_2d = create_fe<2>(fe);
+
+          FullMatrix<double> interpolation_matrix(fe_2d->n_dofs_per_face(0),
+                                                  fe_2d->n_dofs_per_face(0));
+
+          fe_2d->get_subface_interpolation_matrix(*fe_2d,
+                                                  0,
+                                                  interpolation_matrix,
+                                                  0);
+
+          ElementType               element_type;
+          std::vector<unsigned int> scalar_lexicographic;
+          std::vector<unsigned int> lexicographic_numbering;
+
+          get_element_type_specific_information(*fe_1d,
+                                                *fe_1d,
+                                                0,
+                                                element_type,
+                                                scalar_lexicographic,
+                                                lexicographic_numbering);
+
+          subface_interpolation_matrix.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)
+              subface_interpolation_matrix[c] =
+                interpolation_matrix(scalar_lexicographic[i],
+                                     scalar_lexicographic[j]);
+        }
+
       // get gradient and Hessian transformation matrix for the polynomial
       // space associated with the quadrature rule (collocation space). We
       // need to avoid the case with more than a few hundreds of quadrature