diff --git a/hog/forms.py b/hog/forms.py
index 113f035fcb8f8a2bf7f29e23a192bbd7a1efb5f4..b634e3aacc93de873b1fe4e5656926c67ea4f3a7 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -30,6 +30,8 @@ from hog.fem_helpers import (
     element_matrix_iterator,
     scalar_space_dependent_coefficient,
     vector_space_dependent_coefficient,
+    fem_function_on_element,
+    fem_function_gradient_on_element,
 )
 from hog.function_space import FunctionSpace, N1E1Space
 from hog.math_helpers import dot, inv, abs, det, double_contraction
@@ -923,6 +925,7 @@ Weak formulation
         },
     )
 
+
 def grad_rho_by_rho_dot_u(
     trial: FunctionSpace,
     test: FunctionSpace,
@@ -945,7 +948,9 @@ Weak formulation
     ∫ ((∇ρ / ρ) · u) v
 """
 
-    with TimedLogger("assembling grad_rho_by_rho_dot_u-mass matrix", level=logging.DEBUG):
+    with TimedLogger(
+        "assembling grad_rho_by_rho_dot_u-mass matrix", level=logging.DEBUG
+    ):
         tabulation = Tabulation(symbolizer)
 
         jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
@@ -956,7 +961,9 @@ Weak formulation
         else:
             affine_coords = trafo_ref_to_affine(geometry, symbolizer)
             # jac_blending = blending.jacobian(affine_coords)
-            jac_blending_inv = symbolizer.jac_affine_to_blending_inv(geometry.dimensions)
+            jac_blending_inv = symbolizer.jac_affine_to_blending_inv(
+                geometry.dimensions
+            )
             jac_blending_det = symbolizer.abs_det_jac_affine_to_blending()
 
         # jac_blending_det = abs(det(jac_blending))
@@ -996,13 +1003,29 @@ Weak formulation
                 #     sp.Matrix([phi * psi * jac_affine_det]),
                 # )[0]
                 if blending == IdentityMap():
-                    form = dot(((jac_affine_inv.T * grad_rho) / rho[0]), phi) * psi * jac_affine_det
+                    form = (
+                        dot(((jac_affine_inv.T * grad_rho) / rho[0]), phi)
+                        * psi
+                        * jac_affine_det
+                    )
                 else:
-                    form = dot(((jac_blending_inv.T * jac_affine_inv.T * grad_rho) / rho[0]), phi) * psi * jac_affine_det * jac_blending_det
+                    form = (
+                        dot(
+                            (
+                                (jac_blending_inv.T * jac_affine_inv.T * grad_rho)
+                                / rho[0]
+                            ),
+                            phi,
+                        )
+                        * psi
+                        * jac_affine_det
+                        * jac_blending_det
+                    )
                 mat[data.row, data.col] = form
 
     return Form(mat, tabulation, symmetric=trial == test, docstring=docstring)
 
+
 def zero_form(
     trial: FunctionSpace,
     test: FunctionSpace,
diff --git a/hog/integrand.py b/hog/integrand.py
index e7c52bfec9c322b4033cce757a00bec0885bdff8..ce25acc178e99084e6982d308f7a9ed12eff0e83 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -79,7 +79,7 @@ class Form:
     integrand: sp.MatrixBase
     tabulation: Tabulation
     symmetric: bool
-    free_symbols: List[sp.Symbol]
+    free_symbols: List[sp.Symbol] = field(default_factory=lambda: list())
     docstring: str = ""
 
     def integrate(self, quad: Quadrature, symbolizer: Symbolizer) -> sp.Matrix:
@@ -415,12 +415,12 @@ def process_integrand(
 
         mat[data.row, data.col] = integrand(**asdict(s))
 
-    free_symbols = sorted(list(free_symbols), key=lambda x: str(x))
+    free_symbols_sorted = sorted(list(free_symbols), key=lambda x: str(x))
 
     return Form(
         mat,
         tabulation,
         symmetric=is_symmetric,
-        free_symbols=free_symbols,
+        free_symbols=free_symbols_sorted,
         docstring=docstring,
     )
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index 0dec0085b599f874bc146dab9ac6836558f2fce3..e3b1e0dc6ad73855ce221884234d2f78e2b17381 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -665,18 +665,20 @@ class HyTeGElementwiseOperator:
                     )
 
         # Free symbols that shall be settable through the ctor.
-        free_symbol_vars = set()
+        free_symbol_vars_set = set()
         for integration_infos in self.integration_infos.values():
             for integration_info in integration_infos:
                 for fs in integration_info.free_symbols:
-                    free_symbol_vars.add(f"{str(fs)}_{integration_info.integrand_name}")
+                    free_symbol_vars_set.add(
+                        f"{str(fs)}_{integration_info.integrand_name}"
+                    )
 
         free_symbol_vars = [
             CppVariable(
                 name=fs,
                 type=str(self._type_descriptor.pystencils_type),
             )
-            for fs in free_symbol_vars
+            for fs in free_symbol_vars_set
         ]
 
         free_symbol_vars = sorted(free_symbol_vars, key=lambda x: x.name)