diff --git a/firedrake/bcs.py b/firedrake/bcs.py
index ae65282629..ae20ab0e2c 100644
--- a/firedrake/bcs.py
+++ b/firedrake/bcs.py
@@ -137,14 +137,16 @@ 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 fe.complex.is_macrocell() and self._function_space.ufl_element().sobolev_space == ufl.H1:
+                # Skip derivative nodes for supersmooth H1 functions
+                nodes = fe.fiat_equivalent.dual_basis()
+                deriv_nodes = [i for i, node in enumerate(nodes)
+                               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]
diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py
index a982795b26..1ea27a6663 100644
--- a/firedrake/mg/embedded.py
+++ b/firedrake/mg/embedded.py
@@ -14,7 +14,8 @@
 
 native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement"])
 alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced-Hsieh-Clough-Tocher", "Johnson-Mercier",
-                             "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu"])
+                             "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu",
+                             "Guzman-Neilan", "Guzman-Neilan Bubble"])
 non_native_variants = frozenset(["integral", "fdm", "alfeld"])
 
 
diff --git a/tests/macro/test_macro_interp_project.py b/tests/macro/test_macro_interp_project.py
index d7518819a5..12d175fdb2 100644
--- a/tests/macro/test_macro_interp_project.py
+++ b/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)