From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 1 Nov 2022 21:22:37 +0000 (+0100)
Subject: Fast hanging-node algorithm for FE_Q_iso_Q1
X-Git-Tag: v9.5.0-rc1~880^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14385%2Fhead;p=dealii.git

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 b5cc0cf6bb..3942abed73 100644
--- 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 b9a7ba4419..3c73f22759 100644
--- 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];
               }