diff --git a/generate_all_operators.py b/generate_all_operators.py
index af0df0d527b10e42fa51c1b0e0ed256b134066b5..7f0c04a116b8d1fd9073b54a4c746274df4479ed 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -210,7 +210,7 @@ def parse_arguments():
         "--blending",
         default=str(IdentityMap()),
         help=f"Enter a blending type for the operator. Note some geometry maps only work in 2D or 3D respectively. "
-        f"Available geometry maps: {[str(m) for m in ALL_GEOMETRY_MAPS]}",
+             f"Available geometry maps: {[str(m) for m in ALL_GEOMETRY_MAPS]}",
         # TODO add choices list and type=string.lower
     )
 
@@ -257,7 +257,7 @@ def parse_arguments():
         "--clang-format-binary",
         default="clang-format",
         help=f"Allows to specify the name of the clang-format binary and/or optionally its full path."
-        f" By default 'clang-format' will be used.",
+             f" By default 'clang-format' will be used.",
     )
 
     args = parser.parse_args()
@@ -362,9 +362,9 @@ def main():
     re_form = re.compile(args.form)
 
     if (
-        not args.list_operators
-        and args.quad_degree == [-1, -1]
-        and args.quad_rule is None
+            not args.list_operators
+            and args.quad_degree == [-1, -1]
+            and args.quad_rule is None
     ):
         parser.error("Either --quad-degree or --quad-rule are required.")
     quad = {
@@ -392,7 +392,7 @@ def main():
     filtered_operators = list(
         filter(
             lambda o: re.fullmatch(args.space_mapping, o.mapping, re.IGNORECASE)
-            and re.fullmatch(args.form, o.name, re.IGNORECASE),
+                      and re.fullmatch(args.form, o.name, re.IGNORECASE),
             operators,
         )
     )
@@ -499,34 +499,35 @@ class OperatorInfo:
     trial_space: TrialSpace
     test_space: TestSpace
     form: (
-        Callable[
-            [
-                TrialSpace,
-                TestSpace,
-                ElementGeometry,
-                Symbolizer,
-                GeometryMap,
-            ],
-            Form,
-        ]
-        | None
+            Callable[
+                [
+                    TrialSpace,
+                    TestSpace,
+                    ElementGeometry,
+                    Symbolizer,
+                    GeometryMap,
+                ],
+                Form,
+            ]
+            | None
     )
     type_descriptor: HOGType
     form_boundary: (
-        Callable[
-            [
-                TrialSpace,
-                TestSpace,
-                ElementGeometry,
-                ElementGeometry,
-                Symbolizer,
-                GeometryMap,
-            ],
-            Form,
-        ]
-        | None
+            Callable[
+                [
+                    TrialSpace,
+                    TestSpace,
+                    ElementGeometry,
+                    ElementGeometry,
+                    Symbolizer,
+                    GeometryMap,
+                ],
+                Form,
+            ]
+            | None
     ) = None
-    kernel_types: List[KernelWrapperType] = None  # type: ignore[assignment] # will definitely be initialized in __post_init__
+    kernel_types: List[
+        KernelWrapperType] = None  # type: ignore[assignment] # will definitely be initialized in __post_init__
     geometries: Sequence[ElementGeometry] = field(
         default_factory=lambda: [TriangleElement(), TetrahedronElement()]
     )
@@ -575,11 +576,11 @@ class OperatorInfo:
 
 
 def all_operators(
-    symbolizer: Symbolizer,
-    opts: List[Tuple[Set[Opts], LoopStrategy, str]],
-    type_descriptor: HOGType,
-    blending: GeometryMap,
-    geometries: Set[Union[TriangleElement, TetrahedronElement]],
+        symbolizer: Symbolizer,
+        opts: List[Tuple[Set[Opts], LoopStrategy, str]],
+        type_descriptor: HOGType,
+        blending: GeometryMap,
+        geometries: Set[Union[TriangleElement, TetrahedronElement]],
 ) -> List[OperatorInfo]:
     P1 = LagrangianFunctionSpace(1, symbolizer)
     P1Vector = TensorialVectorFunctionSpace(P1)
@@ -622,35 +623,40 @@ def all_operators(
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
     ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkin", TrialSpace(P1),
                             TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin,
-                            coefficient_function_space=P1, only_newton_galerkin_part_of_form=False),
+                                                        coefficient_function_space=P1,
+                                                        only_newton_galerkin_part_of_form=False),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
     ops.append(OperatorInfo("P1Vector", "Diffusion", TrialSpace(P1Vector), TestSpace(P1Vector),
-                            form=diffusion, type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+                            form=diffusion, type_descriptor=type_descriptor, geometries=list(geometries), opts=opts,
+                            blending=blending))
 
-    ops.append(OperatorInfo("P2", "SUPGDiffusion", TrialSpace(P2), TestSpace(P2), 
-                            form=partial(supg_diffusion, velocity_function_space=P2, diffusivityXdelta_function_space=P2), 
+    ops.append(OperatorInfo("P2", "SUPGDiffusion", TrialSpace(P2), TestSpace(P2),
+                            form=partial(supg_diffusion, velocity_function_space=P2,
+                                         diffusivityXdelta_function_space=P2),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
-    
-    ops.append(OperatorInfo("P2", "AdvectionConstCp", TrialSpace(P2), TestSpace(P2), 
-                            form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2, constant_cp = True), 
+
+    ops.append(OperatorInfo("P2", "AdvectionConstCp", TrialSpace(P2), TestSpace(P2),
+                            form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2,
+                                         constant_cp=True),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
-    ops.append(OperatorInfo("P2", "AdvectionVarCp", TrialSpace(P2), TestSpace(P2), 
-                            form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2), 
+    ops.append(OperatorInfo("P2", "AdvectionVarCp", TrialSpace(P2), TestSpace(P2),
+                            form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
-    
-    ops.append(OperatorInfo("P2", "SUPGAdvection", TrialSpace(P2), TestSpace(P2), 
-                            form=partial(supg_advection, velocity_function_space=P2, coefficient_function_space=P2), 
+
+    ops.append(OperatorInfo("P2", "SUPGAdvection", TrialSpace(P2), TestSpace(P2),
+                            form=partial(supg_advection, velocity_function_space=P2, coefficient_function_space=P2),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
     ops.append(OperatorInfo("P2", "BoundaryMass", TrialSpace(P2), TestSpace(P2), form=None,
                             form_boundary=mass_boundary, type_descriptor=type_descriptor, geometries=list(geometries),
                             opts=opts, blending=blending))
 
-    ops.append(OperatorInfo("P1Vector", "LinearElasticity", TrialSpace(P1Vector), TestSpace(P1Vector), form=linear_elasticity
-                            , form_boundary=mass_boundary, type_descriptor=type_descriptor,
-                            geometries=list(geometries), opts=opts, blending=blending))
+    ops.append(OperatorInfo("P1Vector", "LinearElasticity", TrialSpace(P1Vector), TestSpace(P1Vector),
+                            form=partial(linear_elasticity, coefficient_function_space=P1)
+                            , type_descriptor=type_descriptor, geometries=list(geometries), opts=opts,
+                            blending=blending))
 
     # fmt: on
 
@@ -823,16 +829,16 @@ def all_operators(
 
 
 def generate_elementwise_op(
-    args: argparse.Namespace,
-    symbolizer: Symbolizer,
-    op_info: OperatorInfo,
-    name: str,
-    optimizations: Set[Opts],
-    loop_strategy: LoopStrategy,
-    blending: GeometryMap,
-    quad_info: dict[ElementGeometry, int | str],
-    quad_info_boundary: dict[ElementGeometry, int | str],
-    type_descriptor: HOGType,
+        args: argparse.Namespace,
+        symbolizer: Symbolizer,
+        op_info: OperatorInfo,
+        name: str,
+        optimizations: Set[Opts],
+        loop_strategy: LoopStrategy,
+        blending: GeometryMap,
+        quad_info: dict[ElementGeometry, int | str],
+        quad_info_boundary: dict[ElementGeometry, int | str],
+        type_descriptor: HOGType,
 ) -> None:
     """Generates a single operator and writes it to the HyTeG directory."""
 
diff --git a/hog/multi_assignment.py b/hog/multi_assignment.py
index 065090696352b1fb052b79eb9857e648590a98ed..64fdb48389d1ad64a5e7a594428e28329ab01974 100644
--- a/hog/multi_assignment.py
+++ b/hog/multi_assignment.py
@@ -163,7 +163,7 @@ class MultiAssignment(sp.Function):
 
         arg = args[0]
         if not isinstance(arg, sp.Symbol):
-            raise HOGException("First argument of MultiAssignment must be a symbol.")
+            raise HOGException("First argument of MultiAssignment must be a symbol, is " + str(arg))
 
         arg = args[1]
         if type(arg) != int or arg < 0:
diff --git a/hog/recipes/integrands/volume/elasticity_affine.py b/hog/recipes/integrands/volume/elasticity_affine.py
index 8d9c4edf9e1082b9b9dc448c2281b4f0009e14d1..e025dd7e83c72533f8b1e308f40d36eaace04e0e 100644
--- a/hog/recipes/integrands/volume/elasticity_affine.py
+++ b/hog/recipes/integrands/volume/elasticity_affine.py
@@ -39,6 +39,5 @@ def integrand(
     sigma_u = k["lambda"] * grad_u_chain.trace() * sp.eye(grad_u_chain.shape[0]) + k["mu"] * 2 * symm_grad_u
     symm_grad_v = symm_grad(grad_v_chain)
 
-    return tabulate(
-        double_contraction(sigma_u, symm_grad_v) * jac_a_abs_det
-    )
+    return double_contraction(sigma_u, grad_v_chain) * jac_a_abs_det
+