Fast hanging-node algorithm for FE_Q_iso_Q1

diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h
index b5cc0cf6bb21a69c9712aa3028a9479c7eabb629..3942abed73d84df4917a598c18b7b5ae6344bd7a 100644 (file)
--- a/include/deal.II/matrix_free/hanging_nodes_internal.h
+++ b/include/deal.II/matrix_free/hanging_nodes_internal.h
@@ -24,6 +24,7 @@
 #include <deal.II/dofs/dof_accessor.h>
 
 #include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
 #include <deal.II/fe/fe_tools.h>
 
 #include <deal.II/hp/fe_collection.h>
@@ -491,9 +492,13 @@ namespace internal
             for (unsigned int c = 0;
                  c < fe_collection[i].element_multiplicity(base_element_index);
                  ++c, ++comp)
-              if (dim == 1 || dynamic_cast<const FE_Q<dim> *>(
-                                &fe_collection[i].base_element(
-                                  base_element_index)) == nullptr)
+              if (dim == 1 ||
+                  (dynamic_cast<const FE_Q<dim> *>(
+                     &fe_collection[i].base_element(base_element_index)) ==
+                     nullptr &&
+                   dynamic_cast<const FE_Q_iso_Q1<dim> *>(
+                     &fe_collection[i].base_element(base_element_index)) ==
+                     nullptr))
                 supported_components[i][comp] = false;
               else
                 supported_components[i][comp] = true;

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index b9a7ba44199f6d6dd30bbbf3bee82ac34ff9ba7c..3c73f22759dfd0c40027fcbf65f35aff48e95650 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -753,7 +753,8 @@ namespace internal
             }
         }
 
-      if (dim > 1 && dynamic_cast<const FE_Q<dim> *>(&fe))
+      if (dim > 1 && (dynamic_cast<const FE_Q<dim> *>(&fe) ||
+                      dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe)))
         {
           auto &subface_interpolation_matrix_0 =
             univariate_shape_data.subface_interpolation_matrices[0];
@@ -770,19 +771,33 @@ namespace internal
           subface_interpolation_matrix_scalar_0.resize(nn * nn);
           subface_interpolation_matrix_scalar_1.resize(nn * nn);
 
-          std::vector<Point<1>> fe_q_points = QGaussLobatto<1>(nn).get_points();
-          const std::vector<Polynomials::Polynomial<double>> poly =
+          const bool is_feq = dynamic_cast<const FE_Q<dim> *>(&fe) != nullptr;
+
+          std::vector<Point<1>> fe_q_points =
+            is_feq ? QGaussLobatto<1>(nn).get_points() :
+                     QIterated<1>(QTrapezoid<1>(), nn - 1).get_points();
+
+          const std::vector<Polynomials::Polynomial<double>> poly_feq =
             Polynomials::generate_complete_Lagrange_basis(fe_q_points);
 
+          const std::vector<Polynomials::PiecewisePolynomial<double>>
+            poly_feq_iso_q1 =
+              Polynomials::generate_complete_Lagrange_basis_on_subdivisions(nn -
+                                                                              1,
+                                                                            1);
+
           for (unsigned int i = 0, c = 0; i < nn; ++i)
             for (unsigned int j = 0; j < nn; ++j, ++c)
               {
                 subface_interpolation_matrix_scalar_0[c] =
-                  poly[j].value(0.5 * fe_q_points[i][0]);
+                  is_feq ? poly_feq[j].value(0.5 * fe_q_points[i][0]) :
+                           poly_feq_iso_q1[j].value(0.5 * fe_q_points[i][0]);
                 subface_interpolation_matrix_0[c] =
                   subface_interpolation_matrix_scalar_0[c];
                 subface_interpolation_matrix_scalar_1[c] =
-                  poly[j].value(0.5 + 0.5 * fe_q_points[i][0]);
+                  is_feq ?
+                    poly_feq[j].value(0.5 + 0.5 * fe_q_points[i][0]) :
+                    poly_feq_iso_q1[j].value(0.5 + 0.5 * fe_q_points[i][0]);
                 subface_interpolation_matrix_1[c] =
                   subface_interpolation_matrix_scalar_1[c];
               }