diff --git a/generate_all_operators.py b/generate_all_operators.py
index 6ce4a12dcab8ddddd60fb2e05107591d8415a84a..25ca992b12f638af87fac7eb43e8be74affb0405 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -29,7 +29,13 @@ import sympy as sp
 from sympy.core.cache import clear_cache
 from tabulate import tabulate
 
-from hog.blending import GeometryMap, IdentityMap, AnnulusMap, IcosahedralShellMap, ParametricMap
+from hog.blending import (
+    GeometryMap,
+    IdentityMap,
+    AnnulusMap,
+    IcosahedralShellMap,
+    ParametricMap,
+)
 from hog.cse import CseImplementation
 from hog.element_geometry import (
     ElementGeometry,
@@ -46,6 +52,7 @@ from hog.forms import (
     shear_heating,
     epsilon,
     full_stokes,
+    mass,
     nonlinear_diffusion,
     nonlinear_diffusion_newton_galerkin,
     supg_diffusion,
@@ -64,6 +71,7 @@ from hog.function_space import (
     LagrangianFunctionSpace,
     N1E1Space,
     TensorialVectorFunctionSpace,
+    P2PlusBubbleSpace,
     TrialSpace,
     TestSpace,
 )
@@ -584,6 +592,7 @@ def all_operators(
     P1Vector = TensorialVectorFunctionSpace(P1)
     P2 = LagrangianFunctionSpace(2, symbolizer)
     P2Vector = TensorialVectorFunctionSpace(P2)
+    P2PlusBubble = P2PlusBubbleSpace(symbolizer)
     N1E1 = N1E1Space(symbolizer)
 
     two_d = list(geometries & {TriangleElement()})
@@ -593,6 +602,7 @@ def all_operators(
 
     # fmt: off
     # TODO switch to manual specification of opts for now/developement, later use default factory
+    
     ops.append(OperatorInfo("N1E1", "CurlCurl", TrialSpace(N1E1), TestSpace(N1E1), form=curl_curl,
                             type_descriptor=type_descriptor, geometries=three_d, opts=opts, blending=blending))
     ops.append(OperatorInfo("N1E1", "Mass", TrialSpace(N1E1), TestSpace(N1E1), form=mass_n1e1,
@@ -647,6 +657,12 @@ def all_operators(
                             form_boundary=mass_boundary, type_descriptor=type_descriptor, geometries=list(geometries),
                             opts=opts, blending=blending))
 
+    ops.append(OperatorInfo("P2PlusBubble", "Mass", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=mass,
+                            type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
+
+    ops.append(OperatorInfo("P2PlusBubble", "Diffusion", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=diffusion,
+                            type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
+
     # fmt: on
 
     p2vec_epsilon = partial(
diff --git a/hog/operator_generation/indexing.py b/hog/operator_generation/indexing.py
index 896bcae8e6c4ff5ea2d04c69af5b447a2fd9b70d..1680e3942d3a72132392150fc99234f1c487b9b7 100644
--- a/hog/operator_generation/indexing.py
+++ b/hog/operator_generation/indexing.py
@@ -265,7 +265,14 @@ def facedof_index(
     num_microedges_per_edge: sp.Symbol,
 ) -> int:
     """Indexes triangles/faces. Used to compute offsets in volume dof indexing in 2D and AoS layout."""
-    x, y, _ = index
+
+    # Ugly hack; why do we receive an index with only two entries,
+    # and why do we expect one with three?
+    if len(index) == 3:
+        x, y, _ = index
+    elif len(index) == 2:
+        x, y = index
+
     # width = num_faces_per_row_by_type(level, faceType)
     if faceType == FaceType.GRAY:
         return linear_macro_face_index(num_microedges_per_edge, x, y)