diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index 1b8c0fbd7edba7272c8b4c864658964338cee3cd..ac239881a143796dfcf45ca387702118139bfdf4 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -39,6 +39,7 @@ from hog.function_space import (
     P2PlusBubbleSpace,
     TrialSpace,
     TestSpace,
+    TensorialVectorFunctionSpace,
 )
 from hog.forms import (
     mass,
@@ -918,8 +919,8 @@ def form_func(
             raise HOGException("Invalid call to epsilon form.")
         # the input parameters for the epsilon operators are intended to be switched below (col ~ trial component, row ~ test component)
         return epsilon(
-            trial,
-            test,
+            TensorialVectorFunctionSpace(trial, single_component=col),
+            TensorialVectorFunctionSpace(test, single_component=row),
             geometry,
             symbolizer,
             blending=blending,
@@ -937,8 +938,8 @@ def form_func(
         if row not in [0, 1, 2] or col not in [0, 1, 2]:
             raise HOGException("Invalid call to epsilon form.")
         return full_stokes(
-            trial,
-            test,
+            TensorialVectorFunctionSpace(trial, single_component=col),
+            TensorialVectorFunctionSpace(test, single_component=row),
             geometry,
             symbolizer,
             blending=blending,
@@ -1263,7 +1264,7 @@ def main():
     if args.list:
         logger.info("Available forms:")
         for form_info in filtered_form_infos:
-            logger.info(f"- {form_info.full_form_name()}")
+            logger.info(f"- {form_info.full_form_name().sort()}")
         logger.info("Bye.")
         return
 
diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py
index 56e0fae22664f0dc375c988dc26ce7305d6bd337..83dbdc02f87264cd1f64435b32ce25014e93de0a 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -31,7 +31,12 @@ from hog.element_geometry import (
     LineElement,
 )
 from hog.exception import HOGException
-from hog.function_space import FunctionSpace, TrialSpace, TestSpace
+from hog.function_space import (
+    FunctionSpace,
+    TrialSpace,
+    TestSpace,
+    TensorialVectorFunctionSpace,
+)
 from hog.math_helpers import inv, det
 from hog.multi_assignment import MultiAssignment
 from hog.symbolizer import Symbolizer
@@ -49,7 +54,9 @@ from hog.dof_symbol import DoFSymbol
 
 
 def create_empty_element_matrix(
-    trial: TrialSpace, test: TestSpace, geometry: ElementGeometry
+    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
+    test: Union[TestSpace, TensorialVectorFunctionSpace],
+    geometry: ElementGeometry,
 ) -> sp.Matrix:
     """
     Returns a sympy matrix of the required size corresponding to the trial and test spaces, initialized with zeros.
@@ -72,7 +79,9 @@ class ElementMatrixData:
 
 
 def element_matrix_iterator(
-    trial: TrialSpace, test: TestSpace, geometry: ElementGeometry
+    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
+    test: Union[TestSpace, TensorialVectorFunctionSpace],
+    geometry: ElementGeometry,
 ) -> Iterator[ElementMatrixData]:
     """Call this to create a generator to conveniently fill the element matrix."""
     for row, (psi, grad_psi, hessian_psi) in enumerate(
diff --git a/hog/forms.py b/hog/forms.py
index 5ed4b0eeea7aeecc044ad81ba042bb9ff4ad7e0f..9e1d697ba35ea147fbf1df65e7e3e5c259065f3e 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -35,11 +35,12 @@ from hog.fem_helpers import (
 )
 from hog.function_space import (
     FunctionSpace,
+    LagrangianFunctionSpace,
     N1E1Space,
     P2PlusBubbleSpace,
-    TrialSpace,
     TestSpace,
-    LagrangianFunctionSpace,
+    TrialSpace,
+    TensorialVectorFunctionSpace,
 )
 from hog.math_helpers import dot, inv, abs, det, double_contraction
 from hog.quadrature import Quadrature, Tabulation
@@ -319,8 +320,8 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
 
 
 def epsilon(
-    trial: TrialSpace,
-    test: TestSpace,
+    trial: TensorialVectorFunctionSpace,
+    test: TensorialVectorFunctionSpace,
     geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),
@@ -631,8 +632,8 @@ def gradient(
 
 
 def full_stokes(
-    trial: TrialSpace,
-    test: TestSpace,
+    trial: TensorialVectorFunctionSpace,
+    test: TensorialVectorFunctionSpace,
     geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),
diff --git a/hog/integrand.py b/hog/integrand.py
index a95b19358ff6be220e022621941201cea9ac5149..78c999dcdbaa8b805cf23f65b8b5dd9f40447172 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -226,8 +226,8 @@ class IntegrandSymbols:
 
 def process_integrand(
     integrand: Callable[..., Any],
-    trial: TrialSpace,
-    test: TestSpace,
+    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
+    test: Union[TestSpace, TensorialVectorFunctionSpace],
     volume_geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),