Merge pull request #14953 from peterrum/si_is_supported

Work on ShapeInfo::is_supported()
dealii · Mar 27, 2023 · 7174f4d · 7174f4d
2 parents b8135fa + 2c0f378
commit 7174f4d
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 16 deletions.
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
@@ -383,7 +383,20 @@ namespace internal
              const unsigned int                  base_element = 0);
 
       /**
-       * Return which kinds of elements are supported by MatrixFree.
+       * Return whether an element is supported by MatrixFree.
+       *
+       * The following scalar elements are supported:
+       * - FENothing, FE_DGP, and FE_Q_DG0
+       * - polynomial tensor-product elements based on
+       *   Polynomials::Polynomial (FE_Q, FE_DG_Q, FE_DGQArbitraryNodes,
+       *   FE_DGQHermite, FE_DGQLegendre) or Polynomials::PiecewisePolynomial
+       *   (FE_Q_iso_Q1).
+       * - elements for simplex, pyramids, and wedges (FE_SimplexP,
+       *   FE_SimplexDGP, FE_PyramidP, FE_PyramidDGP, FE_WedgeP, FE_WedgeDGP)
+       *
+       * In the case of vectorial elements, FE_RaviartThomasNodal
+       * and FESystem with base elements from the scalar elements
+       * listed above are supported.
        */
       template <int dim, int spacedim>
       static bool

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
@@ -30,6 +30,7 @@
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_dgp.h>
 #include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_pyramid_p.h>
 #include <deal.II/fe/fe_q.h>
@@ -397,14 +398,12 @@ namespace internal
           return;
         }
       else if (quad_in.is_tensor_product() == false ||
-               dynamic_cast<const FE_SimplexP<dim> *>(
-                 &fe_in.base_element(base_element_number)) ||
-               dynamic_cast<const FE_SimplexDGP<dim> *>(
-                 &fe_in.base_element(base_element_number)) ||
-               dynamic_cast<const FE_WedgeP<dim> *>(
-                 &fe_in.base_element(base_element_number)) ||
-               dynamic_cast<const FE_PyramidP<dim> *>(
-                 &fe_in.base_element(base_element_number)))
+               dynamic_cast<const FE_SimplexPoly<dim, dim> *>(
+                 &fe_in.base_element(base_element_number)) != nullptr ||
+               dynamic_cast<const FE_WedgePoly<dim, dim> *>(
+                 &fe_in.base_element(base_element_number)) != nullptr ||
+               dynamic_cast<const FE_PyramidPoly<dim, dim> *>(
+                 &fe_in.base_element(base_element_number)) != nullptr)
         {
           // specialization for arbitrary finite elements and quadrature rules
           // as needed in the context, e.g., of simplices
@@ -1238,6 +1237,9 @@ namespace internal
       if (dim != spacedim)
         return false;
 
+      if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe))
+        return true;
+
       for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
         {
           const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(base));
@@ -1251,13 +1253,17 @@ namespace internal
                 dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr);
               // Simplices are a special case since the polynomial family is not
               // indicative of their support
-              if (dynamic_cast<const FE_SimplexP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_SimplexDGP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_WedgeP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_PyramidP<dim> *>(fe_poly_ptr))
+              if (dynamic_cast<const FE_SimplexPoly<dim, dim> *>(fe_poly_ptr) !=
+                    nullptr ||
+                  dynamic_cast<const FE_WedgePoly<dim, dim> *>(fe_poly_ptr) !=
+                    nullptr ||
+                  dynamic_cast<const FE_PyramidPoly<dim, dim> *>(fe_poly_ptr) !=
+                    nullptr)
                 return true;
 
-              if (dynamic_cast<const TensorProductPolynomials<dim> *>(
+              if (dynamic_cast<const TensorProductPolynomials<
+                      dim,
+                      Polynomials::Polynomial<double>> *>(
                     &fe_poly_ptr->get_poly_space()) == nullptr &&
                   dynamic_cast<const TensorProductPolynomials<
                       dim,
@@ -1269,6 +1275,9 @@ namespace internal
                     nullptr)
                 return false;
             }
+          else if (dynamic_cast<const FE_Nothing<dim, spacedim> *>(fe_ptr) !=
+                   nullptr)
+            return true;
           else
             return false;
         }

diff --git a/tests/matrix_free/is_supported_01.output b/tests/matrix_free/is_supported_01.output
@@ -16,7 +16,7 @@ DEAL::FE_FaceP<2>(1) supported by MatrixFree: false
 DEAL::FE_FaceQ<2>(1) supported by MatrixFree: false
 DEAL::FE_P1NC supported by MatrixFree: false
 DEAL::FE_Nedelec<2>(1) supported by MatrixFree: false
-DEAL::FE_Nothing<2>() supported by MatrixFree: false
+DEAL::FE_Nothing<2>() supported by MatrixFree: true
 DEAL::FE_Q<1,2>(1) supported by MatrixFree: false
 DEAL::FE_Q<2>(1) supported by MatrixFree: true
 DEAL::FE_Q_Bubbles<2>(1) supported by MatrixFree: false
@@ -25,5 +25,5 @@ DEAL::FE_Q_Hierarchical<2>(1) supported by MatrixFree: true
 DEAL::FE_Q_iso_Q1<2>(1) supported by MatrixFree: true
 DEAL::FE_RannacherTurek<2>(0, 2) supported by MatrixFree: false
 DEAL::FE_RaviartThomas<2>(1) supported by MatrixFree: false
-DEAL::FE_RaviartThomasNodal<2>(1) supported by MatrixFree: false
+DEAL::FE_RaviartThomasNodal<2>(1) supported by MatrixFree: true
 DEAL::FE_TraceQ<2>(1) supported by MatrixFree: false