diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index 66d3179d91e34e670ef8adfd6db90b3fe19d4459..aea4a7e4bd43e84c34ebefc5c2023bd2244a6c01 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -14,7 +14,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-from argparse import ArgumentParser
+from argparse import ArgumentParser, RawDescriptionHelpFormatter
 from dataclasses import dataclass, field
 from typing import Callable, Dict, List, Tuple, Union
 import os
@@ -38,6 +38,7 @@ from hog.function_space import (
     P2PlusBubbleSpace,
     TrialSpace,
     TestSpace,
+    TensorialVectorFunctionSpace,
 )
 from hog.forms import (
     mass,
@@ -114,6 +115,9 @@ class FormInfo:
     test_degree: int
     trial_family: str = "Lagrange"
     test_family: str = "Lagrange"
+    supported_geometry_options: List[str] = field(
+        default_factory=lambda: ["triangle", "tetrahedron"]
+    )
     quad_schemes: Dict[int, Union[int, str]] = field(
         default_factory=lambda: {2: 2, 3: 2}
     )
@@ -138,24 +142,24 @@ class FormInfo:
             descr_string = "n1e1"
         elif self.trial_family == "P2 enhanced with Bubble":
             if self.test_family == "P2 enhanced with Bubble":
-                descr_string =  "p2_plus_bubble"
+                descr_string = "p2_plus_bubble"
             elif self.test_family == "DG":
-                descr_string =  f"p2_plus_bubble_to_dg{self.test_degree}"
+                descr_string = f"p2_plus_bubble_to_dg{self.test_degree}"
         elif self.trial_family == "Lagrange" and self.test_family == "Lagrange":
             if self.trial_degree == self.test_degree:
-                descr_string =  f"p{self.trial_degree}"
+                descr_string = f"p{self.trial_degree}"
             else:
-                descr_string =  f"p{self.trial_degree}_to_p{self.test_degree}"
+                descr_string = f"p{self.trial_degree}_to_p{self.test_degree}"
         elif self.trial_family == "DG":
             if self.test_family == "DG":
                 if self.trial_degree == self.test_degree:
-                    descr_string =  f"dg{self.trial_degree}"
+                    descr_string = f"dg{self.trial_degree}"
                 else:
-                    descr_string =  f"dg{self.trial_degree}_to_dg{self.test_degree}"
+                    descr_string = f"dg{self.trial_degree}_to_dg{self.test_degree}"
             elif self.test_family == "Lagrange":
-                descr_string =  f"dg{self.trial_degree}_to_p{self.test_degree}"
+                descr_string = f"dg{self.trial_degree}_to_p{self.test_degree}"
             elif self.test_family == "P2 enhanced with Bubble":
-                descr_string =  f"dg{self.trial_degree}_to_p2_plus_bubble"
+                descr_string = f"dg{self.trial_degree}_to_p2_plus_bubble"
             else:
                 raise HOGException(
                     f"Do not know how to name combination of DGFunctionSpace with {self.test_family}."
@@ -225,6 +229,15 @@ class FormInfo:
             f"quadrature schemes/degree (dimension): {self.quad_schemes}, blending: {self.blending}"
         )
 
+    def supports_geometry(self, geometry: str) -> bool:
+        """Check if form supports a certain type of geometric element."""
+        if geometry == "triangle+tetrahedron":
+            return ("triangle" in self.supported_geometry_options) and (
+                "tetrahedron" in self.supported_geometry_options
+            )
+        else:
+            return geometry in self.supported_geometry_options
+
     def __repr__(self):
         return str(self)
 
@@ -277,6 +290,7 @@ form_infos = [
         test_degree=2,
         trial_family="P2 enhanced with Bubble",
         test_family="P2 enhanced with Bubble",
+        supported_geometry_options=["triangle"],
         quad_schemes={2: "exact"},
         integrate_rows=[],
     ),
@@ -312,6 +326,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 2},
         integrate_rows=[],
     ),
@@ -321,6 +336,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: "exact"},
         integrate_rows=[],
     ),
@@ -330,6 +346,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 2},
         blending=ExternalMap(),
         integrate_rows=[],
@@ -340,6 +357,7 @@ form_infos = [
         test_degree=2,
         trial_family="P2 enhanced with Bubble",
         test_family="P2 enhanced with Bubble",
+        supported_geometry_options=["triangle"],
         quad_schemes={2: "exact"},
         integrate_rows=[],
     ),
@@ -349,6 +367,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 0},
         integrate_rows=[],
     ),
@@ -358,6 +377,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 2},
         blending=ExternalMap(),
         integrate_rows=[],
@@ -473,6 +493,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 6},
         description="Implements a linear form of type: (k(x), psi) where psi a test function and k = k(x) a vectorial, external function.",
         integrate_rows=[],
@@ -483,6 +504,7 @@ form_infos = [
         test_degree=1,
         trial_family="N1E1",
         test_family="N1E1",
+        supported_geometry_options=["tetrahedron"],
         quad_schemes={3: 6},
         blending=ExternalMap(),
         description="Implements a linear form of type: (k(x), psi) where psi a test function and k = k(x) a vectorial, external function.",
@@ -506,27 +528,36 @@ form_infos = [
         "laplace_beltrami",
         trial_degree=1,
         test_degree=1,
+        supported_geometry_options=["embedded_triangle"],
         quad_schemes={3: 1},
     ),
     FormInfo(
         "laplace_beltrami",
         trial_degree=1,
         test_degree=1,
+        supported_geometry_options=["embedded_triangle"],
         quad_schemes={3: 3},
         blending=ExternalMap(),
     ),
-    FormInfo("manifold_mass", trial_degree=1, test_degree=1, quad_schemes={2: 1}),
     FormInfo(
         "manifold_mass",
         trial_degree=1,
         test_degree=1,
+        supported_geometry_options=["embedded_triangle"],
+        quad_schemes={2: 1},
+    ),
+    FormInfo(
+        "manifold_mass",
+        trial_degree=1,
+        test_degree=1,
+        supported_geometry_options=["embedded_triangle"],
         quad_schemes={2: 3},
         blending=ExternalMap(),
     ),
 ]
 
 for d in [1, 2]:
-    for epstype in ["cc", "var"]:
+    for epstype in ["var"]:
         form_infos.append(
             FormInfo(
                 f"epsilon{epstype}",
@@ -540,7 +571,7 @@ for d in [1, 2]:
         )
 
 for d in [1, 2]:
-    for epstype in ["cc", "var"]:
+    for epstype in ["var"]:
         form_infos.append(
             FormInfo(
                 f"epsilon{epstype}",
@@ -668,6 +699,7 @@ for trial_deg, test_deg, transpose in [
                     test_degree=test_deg,
                     trial_family="P2 enhanced with Bubble",
                     test_family="DG",
+                    supported_geometry_options=["triangle"],
                     quad_schemes={
                         2: 3 + test_deg - 1,
                         3: 4 + test_deg - 1,
@@ -686,6 +718,7 @@ for trial_deg, test_deg, transpose in [
                     test_degree=test_deg,
                     trial_family="DG",
                     test_family="P2 enhanced with Bubble",
+                    supported_geometry_options=["triangle"],
                     quad_schemes={
                         2: trial_deg + 3 - 1,
                         3: trial_deg + 4 - 1,
@@ -703,6 +736,7 @@ for blending in [IdentityMap(), ExternalMap()]:
             "manifold_vector_mass",
             trial_degree=2,
             test_degree=2,
+            supported_geometry_options=["embedded_triangle"],
             quad_schemes={2: 3},
             row_dim=3,
             col_dim=3,
@@ -716,6 +750,7 @@ for blending in [IdentityMap(), ExternalMap()]:
             "manifold_normal_penalty",
             trial_degree=2,
             test_degree=2,
+            supported_geometry_options=["embedded_triangle"],
             quad_schemes={2: 3},
             row_dim=3,
             col_dim=3,
@@ -729,6 +764,7 @@ for blending in [IdentityMap(), ExternalMap()]:
             "manifold_epsilon",
             trial_degree=2,
             test_degree=2,
+            supported_geometry_options=["embedded_triangle"],
             quad_schemes={2: 3},
             row_dim=3,
             col_dim=3,
@@ -742,6 +778,7 @@ for blending in [IdentityMap(), ExternalMap()]:
             "manifold_epsilon",
             trial_degree=2,
             test_degree=2,
+            supported_geometry_options=["embedded_triangle"],
             quad_schemes={2: 6},
             row_dim=3,
             col_dim=3,
@@ -755,6 +792,7 @@ for blending in [IdentityMap(), ExternalMap()]:
             "vector_laplace_beltrami",
             trial_degree=2,
             test_degree=2,
+            supported_geometry_options=["embedded_triangle"],
             quad_schemes={2: 3},
             row_dim=3,
             col_dim=3,
@@ -771,6 +809,7 @@ for trial_deg, test_deg, transpose in [(1, 2, True), (2, 1, False)]:
                     "manifold_div",
                     trial_degree=trial_deg,
                     test_degree=test_deg,
+                    supported_geometry_options=["embedded_triangle"],
                     quad_schemes={2: 3},
                     row_dim=1,
                     col_dim=3,
@@ -784,6 +823,7 @@ for trial_deg, test_deg, transpose in [(1, 2, True), (2, 1, False)]:
                     "manifold_divt",
                     trial_degree=trial_deg,
                     test_degree=test_deg,
+                    supported_geometry_options=["embedded_triangle"],
                     quad_schemes={2: 3},
                     row_dim=3,
                     col_dim=1,
@@ -805,6 +845,7 @@ for trial_deg, test_deg, transpose in [
                     "manifold_vector_div",
                     trial_degree=trial_deg,
                     test_degree=test_deg,
+                    supported_geometry_options=["embedded_triangle"],
                     quad_schemes={2: 3},
                     row_dim=1,
                     col_dim=3,
@@ -818,6 +859,7 @@ for trial_deg, test_deg, transpose in [
                     "manifold_vector_divt",
                     trial_degree=trial_deg,
                     test_degree=test_deg,
+                    supported_geometry_options=["embedded_triangle"],
                     quad_schemes={2: 3},
                     row_dim=3,
                     col_dim=1,
@@ -869,8 +911,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,
+            TrialSpace(TensorialVectorFunctionSpace(trial, single_component=col)),
+            TestSpace(TensorialVectorFunctionSpace(test, single_component=row)),
             geometry,
             symbolizer,
             blending=blending,
@@ -888,8 +930,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,
+            TrialSpace(TensorialVectorFunctionSpace(trial, single_component=col)),
+            TestSpace(TensorialVectorFunctionSpace(test, single_component=row)),
             geometry,
             symbolizer,
             blending=blending,
@@ -1028,7 +1070,18 @@ def form_func(
 
 
 def parse_arguments():
-    parser = ArgumentParser(description="Generates all HyTeG forms.")
+    parser = ArgumentParser(
+        formatter_class=RawDescriptionHelpFormatter,
+        description="""Generates (selected) HyTeG forms
+
+Selection of forms to generate is based on the --filter and --geometry arguments:
+
+* if neither --geometry, nor --filter is given, we generate all defined forms
+* if only --filter is given, all forms matching the filter expression are generated
+* if only --geometry is given, all forms that support the given geometric element type are generated
+* if both --filter and --geometry are given, the intersection of both sets is generated""",
+    )
+
     parser.add_argument(
         "hyteg_base_path",
         type=str,
@@ -1041,23 +1094,28 @@ def parse_arguments():
         "-f",
         "--filter",
         type=str,
-        help="only generate forms that include the passed string (which can be a Python regular expression) (works in combination with --list to show all filtered forms and abort)",
+        help="only generate forms that include the passed string (which can be a Python regular expression)",
         default="",
     )
     parser.add_argument(
         "-g",
         "--geometry",
         type=str,
-        help="build form(s) only for triangle (2D) or tetrahedron (3D) elements; if not speficied we do both",
+        help="build form(s) only for triangle (2D), tetrahedron (3D) elements or embedded triangles (manifolds)",
         nargs="?",
-        choices=["triangle", "tetrahedron", "embedded_triangle", "both"],
-        default="both",
+        choices=[
+            "triangle",
+            "tetrahedron",
+            "triangle+tetrahedron",
+            "embedded_triangle",
+        ],
+        default="",
     )
     parser.add_argument(
         "-l",
         "--list",
         action="store_true",
-        help="list all available forms by name and abort (works in combination with --filter to show all filtered forms and abort)",
+        help="list all forms that would be generated for given values of --filter and --geometry; but do not generate anything",
     )
     parser.add_argument(
         "-o",
@@ -1092,6 +1150,78 @@ def valid_base_dir(hyteg_base_path):
     return True
 
 
+def assemble_list_of_forms_to_generate(
+    form_infos: List[FormInfo], logger: logging.Logger
+) -> List[FormInfo]:
+    """From the list of all defined forms extract those to generate.
+
+    The extraction works as follows, depending on the command-line
+    options given:
+
+    - if neither --geometry, nor --filter was given, we return all defined forms
+    - if only --filter was given, it is applied to the list of all defined forms
+    - if only --geometry was given, we extract all forms that support the given geometry element type
+    - if both --filter and --geometry was given, determine the intersection of both sets
+    """
+
+    form_list: List[FormInfo] = []
+
+    parser, args = parse_arguments()
+
+    if args.geometry == "" and args.filter == "":
+        form_list = form_infos
+        logger.info(
+            "No '--filter' and no '--geometry' given: selecting all forms available"
+        )
+
+    elif args.geometry != "" and args.filter == "":
+        logger.info(
+            f"No '--filter' given, extracting all forms supporting '--geometry {args.geometry}'"
+        )
+        for fi in form_infos:
+            if fi.supports_geometry(args.geometry):
+                form_list.append(fi)
+
+    elif args.geometry == "" and args.filter != "":
+        logger.info(f"Extracting forms based on '--filter {args.filter}'")
+        form_list = [
+            fi for fi in form_infos if re.search(args.filter, fi.full_form_name())
+        ]
+
+    else:
+        logger.info(f"Extracting forms based on '--filter {args.filter}'")
+        form_list = [
+            fi for fi in form_infos if re.search(args.filter, fi.full_form_name())
+        ]
+        logger.info(f"Found {len(form_list)} matching forms")
+        for fi in form_list:
+            logger.info(f"* {fi.full_form_name()}")
+
+        logger.info(f"Checking forms against '--geometry {args.geometry}'")
+
+        aux_list = form_list.copy()
+        for fi in form_list:
+            if not fi.supports_geometry(args.geometry):
+                logger.info(f"* deselecting '{fi.full_form_name()}' again")
+                aux_list.remove(fi)
+        form_list = aux_list
+
+    # sort alphabetically by full form name
+    form_list.sort(key=FormInfo.full_form_name)
+
+    return form_list
+
+
+def geometry_string_to_list(geometry_string: str) -> List[str]:
+    glist: List[str] = []
+    if geometry_string == "triangle+tetrahedron":
+        glist.append("triangle")
+        glist.append("tetrahedron")
+    else:
+        glist.append(geometry_string)
+    return glist
+
+
 def main():
     clear_cache()
 
@@ -1121,9 +1251,7 @@ def main():
     logger.info("### HyTeG Operator Generator ###")
     logger.info("################################")
 
-    filtered_form_infos = [
-        fi for fi in form_infos if re.search(args.filter, fi.full_form_name())
-    ]
+    filtered_form_infos = assemble_list_of_forms_to_generate(form_infos, logger)
 
     if args.list:
         logger.info("Available forms:")
@@ -1140,24 +1268,16 @@ def main():
 
     symbolizer = Symbolizer()
 
-    # determine geometries to use
-    geometries: List[ElementGeometry]
-    if args.geometry == "triangle":
-        logger.info(f"- selected geometry: triangle")
-        geometries = [TriangleElement()]
-    elif args.geometry == "tetrahedron":
-        logger.info(f"- selected geometry: tetrahedron")
-        geometries = [TetrahedronElement()]
-    elif args.geometry == "embedded_triangle":
-        logger.info(f"- selected geometry: embedded triangle")
-        geometries = [TriangleElement(space_dimension=3)]
-    else:
-        logger.info(f"- selected geometries: triangle, tetrahedron")
-        geometries = [TriangleElement(), TetrahedronElement()]
+    # no forms -> nothing to do
+    if len(filtered_form_infos) == 0:
+        logger.info(f"Found no matching forms to generate.")
+        logger.info(f"Bye.")
+        quit()
 
     logger.info(
         f"Generating {len(filtered_form_infos)} form{'s' if len(filtered_form_infos) > 1 else ''}:"
     )
+
     for form_info in filtered_form_infos:
         logger.info(f"- {form_info.full_form_name()}")
 
@@ -1185,6 +1305,30 @@ def main():
 
         form_classes = []
 
+        # determine geometries to use for this form
+        geometries: List[ElementGeometry] = []
+
+        target_geometries = []
+        if args.geometry == "":
+            target_geometries = form_info.supported_geometry_options
+        else:
+            target_geometries = geometry_string_to_list(args.geometry)
+
+        if (
+            "triangle" in target_geometries
+            or "triangle+tetrahedron" in target_geometries
+        ):
+            geometries.append(TriangleElement())
+
+        if (
+            "tetrahedron" in target_geometries
+            or "triangle+tetrahedron" in target_geometries
+        ):
+            geometries.append(TetrahedronElement())
+
+        if "embedded_triangle" in target_geometries:
+            geometries.append(TriangleElement(space_dimension=3))
+
         for row in range(0, form_info.row_dim):
             for col in range(0, form_info.col_dim):
                 form_codes = []
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 56e0fae22664f0dc375c988dc26ce7305d6bd337..638e9e31daf90c145038d48fda636b5dfa4cb895 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -31,7 +31,11 @@ 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,
+)
 from hog.math_helpers import inv, det
 from hog.multi_assignment import MultiAssignment
 from hog.symbolizer import Symbolizer
@@ -49,7 +53,9 @@ from hog.dof_symbol import DoFSymbol
 
 
 def create_empty_element_matrix(
-    trial: TrialSpace, test: TestSpace, geometry: ElementGeometry
+    trial: TrialSpace,
+    test: TestSpace,
+    geometry: ElementGeometry,
 ) -> sp.Matrix:
     """
     Returns a sympy matrix of the required size corresponding to the trial and test spaces, initialized with zeros.
@@ -72,7 +78,9 @@ class ElementMatrixData:
 
 
 def element_matrix_iterator(
-    trial: TrialSpace, test: TestSpace, geometry: ElementGeometry
+    trial: TrialSpace,
+    test: TestSpace,
+    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..009a7bd296247a1e29a5c33bc2052f22bf650bbb 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -35,11 +35,11 @@ from hog.fem_helpers import (
 )
 from hog.function_space import (
     FunctionSpace,
+    LagrangianFunctionSpace,
     N1E1Space,
     P2PlusBubbleSpace,
-    TrialSpace,
     TestSpace,
-    LagrangianFunctionSpace,
+    TrialSpace,
 )
 from hog.math_helpers import dot, inv, abs, det, double_contraction
 from hog.quadrature import Quadrature, Tabulation
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)