diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index 200938907e74fa18bce3bf6ada33abf23d677611..f875202c65756e4c8a80c9d15998a2487d7e44e9 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -919,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(
-            TensorialVectorFunctionSpace(trial, single_component=col),
-            TensorialVectorFunctionSpace(test, single_component=row),
+            TrialSpace(TensorialVectorFunctionSpace(trial, single_component=col)),
+            TestSpace(TensorialVectorFunctionSpace(test, single_component=row)),
             geometry,
             symbolizer,
             blending=blending,
@@ -938,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(
-            TensorialVectorFunctionSpace(trial, single_component=col),
-            TensorialVectorFunctionSpace(test, single_component=row),
+            TrialSpace(TensorialVectorFunctionSpace(trial, single_component=col)),
+            TestSpace(TensorialVectorFunctionSpace(test, single_component=row)),
             geometry,
             symbolizer,
             blending=blending,
diff --git a/hog/code_generation.py b/hog/code_generation.py
index afa5aaea38a5a0a1a11b1ec72239c4cb747b9234..fdf086c01d966bcf6ac958d6042acdc2a8c0a5fd 100644
--- a/hog/code_generation.py
+++ b/hog/code_generation.py
@@ -160,9 +160,9 @@ def replace_multi_assignments(
         # Actually filling the dict.
         for ma in multi_assignments:
             replacement_symbol = replacement_symbols[ma.output_arg()]
-            multi_assignments_replacement_symbols[ma.unique_identifier] = (
-                replacement_symbol
-            )
+            multi_assignments_replacement_symbols[
+                ma.unique_identifier
+            ] = replacement_symbol
 
     if multi_assignments_replacement_symbols:
         with TimedLogger(
diff --git a/hog/cpp_printing.py b/hog/cpp_printing.py
index 8093439c49e7a5fa763cd2e8f264c64b2af400dd..a99b736674c095d842a7814bc47ebdb28398afea 100644
--- a/hog/cpp_printing.py
+++ b/hog/cpp_printing.py
@@ -746,7 +746,7 @@ def apply_clang_format(
         if not fail_if_no_binary:
             return
         else:
-            raise HOGException( f"Could not find clang-format binary '{binary}'.")
+            raise HOGException(f"Could not find clang-format binary '{binary}'.")
 
     cmd = [binary, "-i", cpp_file_path]
 
diff --git a/hog/element_geometry.py b/hog/element_geometry.py
index 9fe91fc1b7a7e0f886003a20beec765e4f0862c7..1a8aad207d5c307de8cc10f4a519c542f4258327 100644
--- a/hog/element_geometry.py
+++ b/hog/element_geometry.py
@@ -19,7 +19,6 @@ from hog.exception import HOGException
 
 class ElementGeometry:
     def __init__(self, dimensions: int, num_vertices: int, space_dimension: int):
-
         if space_dimension < dimensions:
             raise HOGException(
                 "The space dimension should be larger or equal to the dimension of the geometry."
diff --git a/hog/exception.py b/hog/exception.py
index 0f50d0e20a10ec80cd4ab16d6041ad4c54fbeadd..d1627d9401f2cefddec652b7738dc618e1e08253 100644
--- a/hog/exception.py
+++ b/hog/exception.py
@@ -14,5 +14,6 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+
 class HOGException(Exception):
-    pass
\ No newline at end of file
+    pass
diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py
index 83dbdc02f87264cd1f64435b32ce25014e93de0a..638e9e31daf90c145038d48fda636b5dfa4cb895 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -35,7 +35,6 @@ from hog.function_space import (
     FunctionSpace,
     TrialSpace,
     TestSpace,
-    TensorialVectorFunctionSpace,
 )
 from hog.math_helpers import inv, det
 from hog.multi_assignment import MultiAssignment
@@ -54,8 +53,8 @@ from hog.dof_symbol import DoFSymbol
 
 
 def create_empty_element_matrix(
-    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
-    test: Union[TestSpace, TensorialVectorFunctionSpace],
+    trial: TrialSpace,
+    test: TestSpace,
     geometry: ElementGeometry,
 ) -> sp.Matrix:
     """
@@ -79,8 +78,8 @@ class ElementMatrixData:
 
 
 def element_matrix_iterator(
-    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
-    test: Union[TestSpace, TensorialVectorFunctionSpace],
+    trial: TrialSpace,
+    test: TestSpace,
     geometry: ElementGeometry,
 ) -> Iterator[ElementMatrixData]:
     """Call this to create a generator to conveniently fill the element matrix."""
diff --git a/hog/forms.py b/hog/forms.py
index 9e1d697ba35ea147fbf1df65e7e3e5c259065f3e..009a7bd296247a1e29a5c33bc2052f22bf650bbb 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -40,7 +40,6 @@ from hog.function_space import (
     P2PlusBubbleSpace,
     TestSpace,
     TrialSpace,
-    TensorialVectorFunctionSpace,
 )
 from hog.math_helpers import dot, inv, abs, det, double_contraction
 from hog.quadrature import Quadrature, Tabulation
@@ -320,8 +319,8 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
 
 
 def epsilon(
-    trial: TensorialVectorFunctionSpace,
-    test: TensorialVectorFunctionSpace,
+    trial: TrialSpace,
+    test: TestSpace,
     geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),
@@ -632,8 +631,8 @@ def gradient(
 
 
 def full_stokes(
-    trial: TensorialVectorFunctionSpace,
-    test: TensorialVectorFunctionSpace,
+    trial: TrialSpace,
+    test: TestSpace,
     geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),
diff --git a/hog/integrand.py b/hog/integrand.py
index 78c999dcdbaa8b805cf23f65b8b5dd9f40447172..a95b19358ff6be220e022621941201cea9ac5149 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -226,8 +226,8 @@ class IntegrandSymbols:
 
 def process_integrand(
     integrand: Callable[..., Any],
-    trial: Union[TrialSpace, TensorialVectorFunctionSpace],
-    test: Union[TestSpace, TensorialVectorFunctionSpace],
+    trial: TrialSpace,
+    test: TestSpace,
     volume_geometry: ElementGeometry,
     symbolizer: Symbolizer,
     blending: GeometryMap = IdentityMap(),
diff --git a/hog/logger.py b/hog/logger.py
index 78497076b8d2259d22c8fc4b9c187a82b8c16146..956f9dc2ce00efd917628982f69b1c9626e9a94a 100644
--- a/hog/logger.py
+++ b/hog/logger.py
@@ -53,7 +53,6 @@ def get_logger(level=logging.INFO):
 
 
 class TimedLogger:
-
     LOG_LEVEL = logging.INFO
 
     _CURRENT_DEPTH = 0
diff --git a/hog/multi_assignment.py b/hog/multi_assignment.py
index 065090696352b1fb052b79eb9857e648590a98ed..75fb2b6b5e8b5e581f9428d443bc7544c2e155bf 100644
--- a/hog/multi_assignment.py
+++ b/hog/multi_assignment.py
@@ -78,6 +78,7 @@ class MultiAssignment(sp.Function):
 
         out_1_expr = MyFunc("alpha", out_idx, input_parameters)
     """
+
     unique_identifier: uuid.UUID
 
     @classmethod
@@ -101,7 +102,7 @@ class MultiAssignment(sp.Function):
     def implementation(self) -> str:
         """
         Should be overridden by subclass.
-        
+
         Returns the implementation (only code block) of the C-function.
         """
         raise HOGException("No implementation has been defined.")
@@ -126,15 +127,17 @@ class MultiAssignment(sp.Function):
 
     def variable_name(self) -> str:
         """
-        Returns the name of a specific instance of the function. 
-        If there are e.g. multiple scalar coefficients, both may have the same function_name() 
+        Returns the name of a specific instance of the function.
+        If there are e.g. multiple scalar coefficients, both may have the same function_name()
         but different variable_name() (e.g. 'alpha' and 'beta').
         """
         return self.args[0]
 
     def symbol_name(self, call_id: int, output_arg: int) -> str:
         """Returns a string that serves as a sympy symbol name."""
-        return f"{self.function_name()}_{self.variable_name()}_out{output_arg}_id{call_id}"
+        return (
+            f"{self.function_name()}_{self.variable_name()}_out{output_arg}_id{call_id}"
+        )
 
     def output_arg(self) -> int:
         return self.args[1]
@@ -160,7 +163,6 @@ class MultiAssignment(sp.Function):
         return self.unique_identifier.int
 
     def __new__(cls, *args):
-
         arg = args[0]
         if not isinstance(arg, sp.Symbol):
             raise HOGException("First argument of MultiAssignment must be a symbol.")
diff --git a/hog/operator_generation/loop_strategies.py b/hog/operator_generation/loop_strategies.py
index 0c370e42b7901a15b0fea556bf75365f608078cf..ee29c93e7555d923220e325094ab8418076c6307 100644
--- a/hog/operator_generation/loop_strategies.py
+++ b/hog/operator_generation/loop_strategies.py
@@ -339,7 +339,6 @@ class BOUNDARY(LoopStrategy):
         self.element_loops: Dict[Union[FaceType, CellType], LoopOverCoordinate] = dict()
 
     def create_loop(self, dim, element_index, micro_edges_per_macro_edge):
-
         if dim == 2:
             if self.facet_id not in [0, 1, 2]:
                 raise HOGException("Invalid facet ID for BOUNDARY loop strategy in 2D.")
diff --git a/hog/recipes/__init__.py b/hog/recipes/__init__.py
index 66608a37410f3c7a2e12498fa816fda9799ddac7..4f026e513de97c126ce27a765a931d2f730e1c39 100644
--- a/hog/recipes/__init__.py
+++ b/hog/recipes/__init__.py
@@ -12,4 +12,4 @@
 # GNU General Public License for more details.
 #
 # You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <https://www.gnu.org/licenses/>.
\ No newline at end of file
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
diff --git a/hog/recipes/integrands/boundary/freeslip_nitsche_divergence.py b/hog/recipes/integrands/boundary/freeslip_nitsche_divergence.py
index d85a4c96ad95c2c4abadc8f02f561f623123265f..f5616c6bb663246c880af71ac77b58a8ff4e1f31 100644
--- a/hog/recipes/integrands/boundary/freeslip_nitsche_divergence.py
+++ b/hog/recipes/integrands/boundary/freeslip_nitsche_divergence.py
@@ -18,7 +18,6 @@ from hog.recipes.common import *
 
 
 def integrand(*, u, v, x, jac_a_boundary, jac_b, matrix, **_):
-
     space_dim = len(x)
 
     A = matrix("A", space_dim, space_dim)
diff --git a/hog/recipes/integrands/boundary/freeslip_nitsche_gradient.py b/hog/recipes/integrands/boundary/freeslip_nitsche_gradient.py
index 569e2381a43721e11074ba4c32d1f444993c6ca4..67042401c79c23cd26b54f9c05dea74a6b3dd471 100644
--- a/hog/recipes/integrands/boundary/freeslip_nitsche_gradient.py
+++ b/hog/recipes/integrands/boundary/freeslip_nitsche_gradient.py
@@ -18,7 +18,6 @@ from hog.recipes.common import *
 
 
 def integrand(*, u, v, x, jac_a_boundary, jac_b, matrix, **_):
-
     space_dim = len(x)
 
     A = matrix("A", space_dim, space_dim)
diff --git a/hog/recipes/integrands/volume/divdiv.py b/hog/recipes/integrands/volume/divdiv.py
index 5a417a617fc132e47dfac2f416611bd6d2cf3c0f..a0544046f5467417ab26dc54a0d3a10e644a89d7 100644
--- a/hog/recipes/integrands/volume/divdiv.py
+++ b/hog/recipes/integrands/volume/divdiv.py
@@ -31,4 +31,10 @@ def integrand(
 ):
     div_u = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)).trace()
     div_v = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)).trace()
-    return (k["k"] if "k" in k else 1.0) * div_u * div_v * tabulate(jac_a_abs_det) * jac_b_abs_det
+    return (
+        (k["k"] if "k" in k else 1.0)
+        * div_u
+        * div_v
+        * tabulate(jac_a_abs_det)
+        * jac_b_abs_det
+    )
diff --git a/hog/recipes/integrands/volume/epsilon.py b/hog/recipes/integrands/volume/epsilon.py
index ab7888e56fa474658e72208089d050b8d4fe02b5..36a91c141c5b4a088c7c0b3dd9ac57b7500f8362 100644
--- a/hog/recipes/integrands/volume/epsilon.py
+++ b/hog/recipes/integrands/volume/epsilon.py
@@ -29,7 +29,6 @@ def integrand(
     tabulate,
     **_,
 ):
-
     grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
     grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)
 
diff --git a/hog/recipes/integrands/volume/epsilon_affine.py b/hog/recipes/integrands/volume/epsilon_affine.py
index 049f334a027fa363c38be1f18f59da494cb6c12c..b040f62fac8de2c10d144a2bac9efa7683add599 100644
--- a/hog/recipes/integrands/volume/epsilon_affine.py
+++ b/hog/recipes/integrands/volume/epsilon_affine.py
@@ -27,7 +27,6 @@ def integrand(
     tabulate,
     **_,
 ):
-
     grad_u_chain = jac_a_inv.T * grad_u
     grad_v_chain = jac_a_inv.T * grad_v
 
diff --git a/hog/recipes/integrands/volume/frozen_velocity.py b/hog/recipes/integrands/volume/frozen_velocity.py
index c0b24c1ca0907acf16f0ac4008f3bea9691939cc..dbe4bf0d356b5f040b0f03f0ebc4fa8b5aef2913 100644
--- a/hog/recipes/integrands/volume/frozen_velocity.py
+++ b/hog/recipes/integrands/volume/frozen_velocity.py
@@ -30,6 +30,9 @@ def integrand(
     tabulate,
     **_,
 ):
-    return dot(
-         (jac_b_inv.T * jac_a_inv.T * grad_k["rho"]) / k["rho"], u
-    ) * v * tabulate(jac_a_abs_det) * jac_b_abs_det
+    return (
+        dot((jac_b_inv.T * jac_a_inv.T * grad_k["rho"]) / k["rho"], u)
+        * v
+        * tabulate(jac_a_abs_det)
+        * jac_b_abs_det
+    )
diff --git a/hog/recipes/integrands/volume/full_stokes.py b/hog/recipes/integrands/volume/full_stokes.py
index d8e40a1d0e94c278278de8d0f340c28d0ca34c51..03d37e5aa9be4a7e97c2e575e0e651fa9387693d 100644
--- a/hog/recipes/integrands/volume/full_stokes.py
+++ b/hog/recipes/integrands/volume/full_stokes.py
@@ -29,7 +29,6 @@ def integrand(
     tabulate,
     **_,
 ):
-
     grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
     grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)