test supersmooth bcs

firedrake/bcs.py

@@ -137,14 +137,14 @@ def hermite_stride(bcnodes):
             fe = self._function_space.finat_element
             tdim = self._function_space.mesh().topological_dimension()
             if isinstance(fe, finat.Hermite) and tdim == 1:
-                return bcnodes[::2]  # every second dof is the vertex value
-            elif isinstance(fe, finat.AlfeldSorokina):
-                # Skip derivative nodes
-                deriv_nodes = [k for k, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
-                deriv_ids = self._function_space.cell_node_list[:, deriv_nodes]
-                return np.setdiff1d(bcnodes, deriv_ids)
-            else:
-                return bcnodes
+                bcnodes = bcnodes[::2]  # every second dof is the vertex value
+            elif self._function_space.ufl_element().sobolev_space == ufl.H1:
+                # Skip derivative nodes for supersmooth H1 functions
+                deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
+                if len(deriv_nodes) > 0:
+                    deriv_ids = self._function_space.cell_node_list[:, deriv_nodes]
+                    bcnodes = np.setdiff1d(bcnodes, deriv_ids)
+            return bcnodes
 
         sub_d = (self.sub_domain, ) if isinstance(self.sub_domain, str) else as_tuple(self.sub_domain)
         sub_d = [s if isinstance(s, str) else as_tuple(s) for s in sub_d]

tests/macro/test_macro_interp_project.py

@@ -123,6 +123,7 @@ def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op):
                          [(2, 'Alfeld-Sorokina', 2, 2),
                           (3, 'Alfeld-Sorokina', 2, 2),
                           (2, 'Reduced-Arnold-Qin', 2, 1),
+                          (3, 'Guzman-Neilan', 3, 1),
                           (3, 'Christiansen-Hu', 1, 1),
                           (2, 'Bernardi-Raugel', 2, 1),
                           (3, 'Bernardi-Raugel', 3, 1),
@@ -134,3 +135,45 @@ def test_piola_convergence(hierarchy, dim, el, deg, convrate, op):
     if op == proj:
         convrate += 1
     run_convergence(hierarchy[dim], el, deg, convrate, op)
+
+
+# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions
+@pytest.mark.parametrize('enriched', (False, True))
+def test_supersmooth_bcs(mesh, enriched):
+    tdim = mesh.topological_dimension()
+    if enriched:
+        cell = mesh.ufl_cell()
+        element = EnrichedElement(FiniteElement("AS", cell=cell, degree=2),
+                                  FiniteElement("GNB", cell=cell, degree=tdim))
+        V = FunctionSpace(mesh, element)
+    else:
+        V = FunctionSpace(mesh, "Alfeld-Sorokina", 2)
+
+    # check that V in H1
+    assert V.ufl_element().sobolev_space == H1
+
+    # check that V is supersmooth
+    nodes = V.finat_element.fiat_equivalent.dual.nodes
+    deriv_nodes = [i for i, node in enumerate(nodes) if len(node.deriv_dict)]
+    assert len(deriv_nodes) == tdim + 1
+
+    deriv_ids = V.cell_node_list[:, deriv_nodes]
+    u = Function(V)
+
+    CG = FunctionSpace(mesh, "Lagrange", 2)
+    RT = FunctionSpace(mesh, "RT", 1)
+    for sub in [1, (1, 2), "on_boundary"]:
+        bc = DirichletBC(V, 0, sub)
+
+        # check that we have the expected number of bc nodes
+        nnodes = len(bc.nodes)
+        expected = tdim * len(DirichletBC(CG, 0, sub).nodes)
+        if enriched:
+            expected += len(DirichletBC(RT, 0, sub).nodes)
+        assert nnodes == expected
+
+        # check that the bc does not set the derivative nodes
+        u.assign(111)
+        u.dat.data_wo[deriv_ids] = 42
+        bc.zero(u)
+        assert numpy.allclose(u.dat.data_ro[deriv_ids], 42)