From 566d222d055f4192deda3bfa950e251c586dc372 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 20 Aug 2024 08:49:25 -0600
Subject: [PATCH] Fix FiniteElement::has_support_points() for FE_Nothing.

---
 include/deal.II/fe/fe.h | 31 ++++++++++++++++++++-----------
 source/fe/fe.cc         | 13 ++++++++++++-
 2 files changed, 32 insertions(+), 12 deletions(-)

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index 744bfa0147..2d5dfe2bc0 100644
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -2023,20 +2023,29 @@ public:
 
   /**
    * Return whether a finite element has defined support points. If the result
-   * is true, then a call to the get_unit_support_points() yields a non-empty
-   * array.
-   *
-   * The result may be false if an element is not defined by interpolating
-   * shape functions, for example by P-elements on quadrilaterals. It will
-   * usually only be true if the element constructs its shape functions by the
-   * requirement that they be one at a certain point and zero at all the
-   * points associated with the other shape functions.
+   * is true, then a call to the get_unit_support_points() yields an
+   * array with `dofs_per_cell` entries. Note that the function's name is
+   * poorly chosen: The function does not return *whether* an element has
+   * support points, but whether its implementation is able to *provide*
+   * a list of support points via get_unit_support_points().
+   *
+   * The result may be false if a finite element defines itself not by
+   * interpolating shape functions, but by other means. A typical example are
+   * discontinuous $P$-type elements on quadrilaterals (rather than the common
+   * $Q$-type elements on quadrilaterals). Elements will generally only
+   * `true` if they construct their shape functions by the
+   * requirement that they be nonzero at a certain point and zero at all the
+   * points associated with the other shape functions. In other words,
+   * if the "node functionals" that form the dual space of the finite
+   * element space and are used to define the shape functions are point
+   * evaluations.
    *
    * In composed elements (i.e. for the FESystem class), the result will be
    * true if all the base elements have defined support points. FE_Nothing
-   * is a special case in FESystems, because it has 0 support points and
-   * has_support_points() is false, but an FESystem containing an FE_Nothing
-   * among other elements will return true.
+   * is a special case since it does not have support point (it has no
+   * shape functions, after all), and so returns an empty array from its
+   * get_unit_support_points() function. Nonetheless, because that array has
+   * the right size, it reports `true` in the current function.
    */
   bool
   has_support_points() const;
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index f9b3978cc7..8bae7157dd 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -1044,7 +1044,18 @@ template <int dim, int spacedim>
 bool
 FiniteElement<dim, spacedim>::has_support_points() const
 {
-  return (unit_support_points.size() != 0);
+  if (this->dofs_per_cell > 0)
+    return (unit_support_points.size() != 0);
+  else
+    {
+      // If the FE has no DoFs, we shouldn't expect the array
+      // size to be anything other than zero:
+      AssertDimension(unit_support_points.size(), 0);
+
+      // A finite element without DoFs *has* support points
+      // (which is then an empty array)
+      return true;
+    }
 }
 
 
-- 
2.39.5