diff --git a/generate_all_operators.py b/generate_all_operators.py
index f582363e815ce1073ea0811a30ad8731f28076f3..9842f60ff3f633dce74d6200492082c3e1b9a9c8 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -526,6 +526,7 @@ def all_operators(
     P1 = LagrangianFunctionSpace(1, symbolizer)
     P1Vector = TensorialVectorFunctionSpace(P1)
     P2 = LagrangianFunctionSpace(2, symbolizer)
+    P2Vector = TensorialVectorFunctionSpace(P2)
     N1E1 = N1E1Space(symbolizer)
 
     two_d = [TriangleElement()]
@@ -567,6 +568,24 @@ def all_operators(
 
     # fmt: on
 
+    p2vec_epsilon = partial(
+        epsilon,
+        variable_viscosity=True,
+        coefficient_function_space=P2,
+    )
+    ops.append(
+        OperatorInfo(
+            mapping=f"P2Vector",
+            name=f"Epsilon",
+            trial_space=P2Vector,
+            test_space=P2Vector,
+            form=p2vec_epsilon,
+            type_descriptor=type_descriptor,
+            opts=opts,
+            blending=blending,
+        )
+    )
+
     for c in [0, 1, 2]:
         # fmt: off
         div_geometries = [TriangleElement(), TetrahedronElement()]