diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index d96ca5defea2c8af6a16ac3616d77dc6eb526d6f..cf26ddb885243c1cf9de51fde9a73bebdf3b1995 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -34,6 +34,7 @@ from hog.element_geometry import (
 from hog.function_space import (
     LagrangianFunctionSpace,
     N1E1Space,
+    P2PlusBubbleSpace,
     TrialSpace,
     TestSpace,
 )
@@ -133,6 +134,8 @@ class FormInfo:
         """A compact representation of the function spaces."""
         if self.trial_family == "N1E1":
             return "n1e1"
+        elif self.trial_family == "P2 enhanced with Bubble":
+            return "p2_plus_bubble"
         elif self.trial_degree == self.test_degree:
             return f"p{self.trial_degree}"
         else:
@@ -163,6 +166,8 @@ class FormInfo:
         sub_dir = "p1"
         if self.trial_family == "N1E1":
             sub_dir = "n1e1"
+        elif self.trial_family == "P2 enhanced with Bubble":
+            return "p2_plus_bubble"
         elif self.trial_degree == self.test_degree:
             sub_dir = f"p{self.trial_degree}"
         else:
@@ -242,6 +247,15 @@ form_infos = [
         quad_schemes={2: 3, 3: 3},
         blending=ExternalMap(),
     ),
+    FormInfo(
+        "diffusion",
+        trial_degree=2,
+        test_degree=2,
+        trial_family="P2 enhanced with Bubble",
+        test_family="P2 enhanced with Bubble",
+        quad_schemes={2: "exact", 3: "exact"},
+        integrate_rows=[],
+    ),
     FormInfo(
         "mass",
         trial_degree=1,
@@ -296,6 +310,15 @@ form_infos = [
         blending=ExternalMap(),
         integrate_rows=[],
     ),
+    FormInfo(
+        "mass",
+        trial_degree=2,
+        test_degree=2,
+        trial_family="P2 enhanced with Bubble",
+        test_family="P2 enhanced with Bubble",
+        quad_schemes={2: "exact", 3: "exact"},
+        integrate_rows=[],
+    ),
     FormInfo(
         "curl_curl",
         trial_degree=1,
@@ -1077,6 +1100,8 @@ def main():
         trial: TrialSpace
         if form_info.trial_family == "N1E1":
             trial = TrialSpace(N1E1Space(symbolizer))
+        elif form_info.trial_family == "P2 enhanced with Bubble":
+            trial = TrialSpace(P2PlusBubbleSpace(symbolizer))
         else:
             trial = TrialSpace(
                 LagrangianFunctionSpace(form_info.trial_degree, symbolizer)
@@ -1085,6 +1110,8 @@ def main():
         test: TestSpace
         if form_info.test_family == "N1E1":
             test = TestSpace(N1E1Space(symbolizer))
+        elif form_info.test_family == "P2 enhanced with Bubble":
+            test = TestSpace(P2PlusBubbleSpace(symbolizer))
         else:
             test = TestSpace(LagrangianFunctionSpace(form_info.test_degree, symbolizer))
 
diff --git a/hog/forms.py b/hog/forms.py
index f67f0781d698489310f163fcd996231aecee66e0..ec630ce41999ca80a8396f1d0de0838f231fe435 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -36,6 +36,7 @@ from hog.fem_helpers import (
 from hog.function_space import (
     FunctionSpace,
     N1E1Space,
+    P2PlusBubbleSpace,
     TrialSpace,
     TestSpace,
     LagrangianFunctionSpace,
@@ -1097,7 +1098,7 @@ Weak formulation
 
     ∫ ((∇ρ / ρ) · u) v
 """
-    
+
     from hog.recipes.integrands.volume.frozen_velocity import integrand
 
     return process_integrand(
@@ -1114,6 +1115,7 @@ Weak formulation
         },
     )
 
+
 def zero_form(
     trial: TrialSpace,
     test: TestSpace,
diff --git a/hog/hyteg_form_template.py b/hog/hyteg_form_template.py
index f0d43087871aab1c97978286f5a0a4776f6eafbb..b1c7c34124dd97cb339213ac35a8b6fc11cc2b32 100644
--- a/hog/hyteg_form_template.py
+++ b/hog/hyteg_form_template.py
@@ -21,7 +21,7 @@ from hog.ast import Assignment, CodeBlock
 from hog.exception import HOGException
 from hog.quadrature import Quadrature
 from hog.symbolizer import Symbolizer
-from hog.function_space import N1E1Space, TrialSpace, TestSpace
+from hog.function_space import N1E1Space, P2PlusBubbleSpace, TrialSpace, TestSpace
 from hog.element_geometry import ElementGeometry
 from hog.code_generation import code_block_from_element_matrix
 from hog.multi_assignment import Member
@@ -292,6 +292,8 @@ class HyTeGFormClass:
 
         if isinstance(self.trial, N1E1Space):
             super_class = f"n1e1::N1E1Form"
+        elif isinstance(self.trial, P2PlusBubbleSpace):
+            super_class = f"P2PlusBubbleForm"
         elif self.trial.degree == self.test.degree:
             super_class = f"P{self.trial.degree}FormHyTeG"
         else:
@@ -358,6 +360,8 @@ class HyTeGForm:
 
         if isinstance(self.trial, N1E1Space):
             super_class = f"N1E1Form"
+        elif isinstance(self.trial, P2PlusBubbleSpace):
+            super_class = f"P2PlusBubbleForm"
         elif self.trial.degree == self.test.degree:
             super_class = f"form_hyteg_base/P{self.trial.degree}FormHyTeG"
         else: