diff --git a/hog/operator_generation/__init__.py b/hog/operator_generation/__init__.py
index 7de6b0e232ff9b89bccc8a525c7841951ce1c2ca..c4cab3c872f55c3312a8a9216b435f2daaf541dd 100644
--- a/hog/operator_generation/__init__.py
+++ b/hog/operator_generation/__init__.py
@@ -15,7 +15,7 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 __all__ = [
-    "function_space_impls",
+    "function_space_implementations",
     "indexing",
     "kernel_types",
     "loop_strategies",
diff --git a/hog/operator_generation/function_space_implementations/__init__.py b/hog/operator_generation/function_space_implementations/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..83f757d5d24ff92f53857294feaba132c63d2ecd
--- /dev/null
+++ b/hog/operator_generation/function_space_implementations/__init__.py
@@ -0,0 +1,24 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024  HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# 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/>.
+
+__all__ = [
+    "function_space_impl_base",
+    "function_space_impl_factor",
+    "p0_function_space_impl",
+    "p1_function_space_impl",
+    "p2_function_space_impl",
+    "n1e1_function_space_impl",
+]
diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_implementations/function_space_impl_base.py
similarity index 71%
rename from hog/operator_generation/function_space_impls.py
rename to hog/operator_generation/function_space_implementations/function_space_impl_base.py
index 67ebd430437f5a79015af75b9cc93cf837edf75c..a618d138db3c57c189c733a5858b675ab0e91930 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_implementations/function_space_impl_base.py
@@ -78,64 +78,6 @@ class FunctionSpaceImpl(ABC):
         self.type_descriptor = type_descriptor
         self.is_pointer = is_pointer
 
-    @staticmethod
-    def create_impl(  # type: ignore[no-untyped-def] # returns Self but Self is kind of new
-        func_space: FunctionSpace,
-        name: str,
-        type_descriptor: HOGType,
-        is_pointer: bool = False,
-    ):
-        """Takes a mathematical function space and produces the corresponding function space implementation.
-
-        :param func_space:      The mathematical function space.
-        :param name:            The C++ variable name/identifier of this function.
-        :param type_descriptor: The value type of this function.
-        :param is_pointer:      Whether the C++ variable is of a pointer type. Used
-                                to print member accesses.
-        """
-        impl_class: type
-
-        import hog.operator_generation.function_space_implementations.p0_space_impl
-        import hog.operator_generation.function_space_implementations.p1_space_impl
-        import hog.operator_generation.function_space_implementations.p2_space_impl
-        import hog.operator_generation.function_space_implementations.n1e1_space_impl
-
-        if isinstance(func_space, LagrangianFunctionSpace):
-            if func_space.degree == 2:
-                impl_class = P2FunctionSpaceImpl
-            elif func_space.degree == 1:
-                impl_class = P1FunctionSpaceImpl
-            elif func_space.degree == 0:
-                impl_class = P0FunctionSpaceImpl
-            else:
-                raise HOGException("Lagrangian function space must be of order 1 or 2.")
-        elif isinstance(func_space, TensorialVectorFunctionSpace):
-            if isinstance(func_space.component_function_space, LagrangianFunctionSpace):
-                if func_space.component_function_space.degree == 1:
-                    if func_space.single_component is None:
-                        impl_class = P1VectorFunctionSpaceImpl
-                    else:
-                        impl_class = P1FunctionSpaceImpl
-                elif func_space.component_function_space.degree == 2:
-                    if func_space.single_component is None:
-                        impl_class = P2VectorFunctionSpaceImpl
-                    else:
-                        impl_class = P2FunctionSpaceImpl
-                else:
-                    raise HOGException(
-                        "TensorialVectorFunctionSpaces not supported for the chosen components."
-                    )
-            else:
-                raise HOGException(
-                    "TensorialVectorFunctionSpaces are only supported with Lagrangian component spaces."
-                )
-        elif isinstance(func_space, N1E1Space):
-            impl_class = N1E1FunctionSpaceImpl
-        else:
-            raise HOGException("Unexpected function space")
-
-        return impl_class(func_space, name, type_descriptor, is_pointer)
-
     def _create_field(self, name: str) -> Field:
         """Creates a pystencils field with a given name and stride 1."""
         f = Field.create_generic(
diff --git a/hog/operator_generation/function_space_implementations/function_space_impl_factory.py b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py
new file mode 100644
index 0000000000000000000000000000000000000000..e0d0b1628a185adaaff094613399a7dc146f85f5
--- /dev/null
+++ b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py
@@ -0,0 +1,92 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024  HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# 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/>.
+
+from hog.function_space import (
+    FunctionSpace,
+    LagrangianFunctionSpace,
+    TensorialVectorFunctionSpace,
+    N1E1Space,
+)
+
+from hog.operator_generation.function_space_implementations.p0_space_impl import (
+    P0FunctionSpaceImpl,
+)
+from hog.operator_generation.function_space_implementations.p1_space_impl import (
+    P1FunctionSpaceImpl,
+    P1VectorFunctionSpaceImpl,
+)
+from hog.operator_generation.function_space_implementations.p2_space_impl import (
+    P2FunctionSpaceImpl,
+    P2VectorFunctionSpaceImpl,
+)
+from hog.operator_generation.function_space_implementations.n1e1_space_impl import (
+    N1E1FunctionSpaceImpl,
+)
+
+from hog.operator_generation.types import HOGType
+
+
+def create_impl(
+    func_space: FunctionSpace,
+    name: str,
+    type_descriptor: HOGType,
+    is_pointer: bool = False,
+):
+    """Takes a mathematical function space and produces the corresponding function space implementation.
+
+    :param func_space:      The mathematical function space.
+    :param name:            The C++ variable name/identifier of this function.
+    :param type_descriptor: The value type of this function.
+    :param is_pointer:      Whether the C++ variable is of a pointer type. Used
+                            to print member accesses.
+    """
+    impl_class: type
+
+    if isinstance(func_space, LagrangianFunctionSpace):
+        if func_space.degree == 2:
+            impl_class = P2FunctionSpaceImpl
+        elif func_space.degree == 1:
+            impl_class = P1FunctionSpaceImpl
+        elif func_space.degree == 0:
+            impl_class = P0FunctionSpaceImpl
+        else:
+            raise HOGException("Lagrangian function space must be of order 1 or 2.")
+    elif isinstance(func_space, TensorialVectorFunctionSpace):
+        if isinstance(func_space.component_function_space, LagrangianFunctionSpace):
+            if func_space.component_function_space.degree == 1:
+                if func_space.single_component is None:
+                    impl_class = P1VectorFunctionSpaceImpl
+                else:
+                    impl_class = P1FunctionSpaceImpl
+            elif func_space.component_function_space.degree == 2:
+                if func_space.single_component is None:
+                    impl_class = P2VectorFunctionSpaceImpl
+                else:
+                    impl_class = P2FunctionSpaceImpl
+            else:
+                raise HOGException(
+                    "TensorialVectorFunctionSpaces not supported for the chosen components."
+                )
+        else:
+            raise HOGException(
+                "TensorialVectorFunctionSpaces are only supported with Lagrangian component spaces."
+            )
+    elif isinstance(func_space, N1E1Space):
+        impl_class = N1E1FunctionSpaceImpl
+    else:
+        raise HOGException("Unexpected function space")
+
+    return impl_class(func_space, name, type_descriptor, is_pointer)
diff --git a/hog/operator_generation/function_space_implementations/n1e1_space_impl.py b/hog/operator_generation/function_space_implementations/n1e1_space_impl.py
index c0cdd5a18b25c335f0641dc444a2461b8a200347..16909cd489a0d04d10837b66c2f63b5faba0788b 100644
--- a/hog/operator_generation/function_space_implementations/n1e1_space_impl.py
+++ b/hog/operator_generation/function_space_implementations/n1e1_space_impl.py
@@ -14,8 +14,36 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-
-import hog.operator_generation.function_space_impls
+import itertools
+import numpy as np
+import sympy as sp
+from typing import List, Tuple, Set, Union, Dict
+
+import pystencils as ps
+from pystencils import Field, FieldType
+from pystencils.backends.cbackend import CustomCodeNode
+
+from hog.element_geometry import ElementGeometry
+from hog.exception import HOGException
+from hog.function_space import (
+    FunctionSpace,
+    LagrangianFunctionSpace,
+    TensorialVectorFunctionSpace,
+    N1E1Space,
+)
+from hog.operator_generation.indexing import (
+    CellType,
+    FaceType,
+    VolumeDoFMemoryLayout,
+    micro_element_to_vertex_indices,
+    micro_vertex_to_edge_indices,
+    micro_element_to_volume_indices,
+    IndexingInfo,
+)
+from hog.operator_generation.types import HOGType
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
 
 
 class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
diff --git a/hog/operator_generation/function_space_implementations/p0_space_impl.py b/hog/operator_generation/function_space_implementations/p0_space_impl.py
index 15180b1fcff9f42f2fce45285bff772f5315632c..48987bb8bddb0fdb2ea2b4589c98750b7ebb3b14 100644
--- a/hog/operator_generation/function_space_implementations/p0_space_impl.py
+++ b/hog/operator_generation/function_space_implementations/p0_space_impl.py
@@ -14,7 +14,36 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import hog.operator_generation.function_space_impls
+import itertools
+import numpy as np
+import sympy as sp
+from typing import List, Tuple, Set, Union, Dict
+
+import pystencils as ps
+from pystencils import Field, FieldType
+from pystencils.backends.cbackend import CustomCodeNode
+
+from hog.element_geometry import ElementGeometry
+from hog.exception import HOGException
+from hog.function_space import (
+    FunctionSpace,
+    LagrangianFunctionSpace,
+    TensorialVectorFunctionSpace,
+    N1E1Space,
+)
+from hog.operator_generation.indexing import (
+    CellType,
+    FaceType,
+    VolumeDoFMemoryLayout,
+    micro_element_to_vertex_indices,
+    micro_vertex_to_edge_indices,
+    micro_element_to_volume_indices,
+    IndexingInfo,
+)
+from hog.operator_generation.types import HOGType
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
 
 
 class P0FunctionSpaceImpl(FunctionSpaceImpl):
diff --git a/hog/operator_generation/function_space_implementations/p1_space_impl.py b/hog/operator_generation/function_space_implementations/p1_space_impl.py
index cb851e65873831dd647d949273a203465dbfff5d..0f993dce5f4820d6ab40de9f9bce4ea9d0d56ea9 100644
--- a/hog/operator_generation/function_space_implementations/p1_space_impl.py
+++ b/hog/operator_generation/function_space_implementations/p1_space_impl.py
@@ -14,7 +14,36 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-import hog.operator_generation.function_space_impls
+import itertools
+import numpy as np
+import sympy as sp
+from typing import List, Tuple, Set, Union, Dict
+
+import pystencils as ps
+from pystencils import Field, FieldType
+from pystencils.backends.cbackend import CustomCodeNode
+
+from hog.element_geometry import ElementGeometry
+from hog.exception import HOGException
+from hog.function_space import (
+    FunctionSpace,
+    LagrangianFunctionSpace,
+    TensorialVectorFunctionSpace,
+    N1E1Space,
+)
+from hog.operator_generation.indexing import (
+    CellType,
+    FaceType,
+    VolumeDoFMemoryLayout,
+    micro_element_to_vertex_indices,
+    micro_vertex_to_edge_indices,
+    micro_element_to_volume_indices,
+    IndexingInfo,
+)
+from hog.operator_generation.types import HOGType
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
 
 
 class P1FunctionSpaceImpl(FunctionSpaceImpl):
diff --git a/hog/operator_generation/function_space_implementations/p2_space_impl.py b/hog/operator_generation/function_space_implementations/p2_space_impl.py
index f56bca5a391838b0372a02a610f7d6e7a2791b43..01eec61b8cca34a84bb54bf04f87cba9f34d5219 100644
--- a/hog/operator_generation/function_space_implementations/p2_space_impl.py
+++ b/hog/operator_generation/function_space_implementations/p2_space_impl.py
@@ -14,8 +14,36 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
-
-import hog.operator_generation.function_space_impls
+import itertools
+import numpy as np
+import sympy as sp
+from typing import List, Tuple, Set, Union, Dict
+
+import pystencils as ps
+from pystencils import Field, FieldType
+from pystencils.backends.cbackend import CustomCodeNode
+
+from hog.element_geometry import ElementGeometry
+from hog.exception import HOGException
+from hog.function_space import (
+    FunctionSpace,
+    LagrangianFunctionSpace,
+    TensorialVectorFunctionSpace,
+    N1E1Space,
+)
+from hog.operator_generation.indexing import (
+    CellType,
+    FaceType,
+    VolumeDoFMemoryLayout,
+    micro_element_to_vertex_indices,
+    micro_vertex_to_edge_indices,
+    micro_element_to_volume_indices,
+    IndexingInfo,
+)
+from hog.operator_generation.types import HOGType
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
 
 
 class P2FunctionSpaceImpl(FunctionSpaceImpl):
diff --git a/hog/operator_generation/kernel_types.py b/hog/operator_generation/kernel_types.py
index f0ba9f13523100cb3c2862dd1bcce32faea0644f..712dc0365b750c6bc6b95088e39133bcf71fde6b 100644
--- a/hog/operator_generation/kernel_types.py
+++ b/hog/operator_generation/kernel_types.py
@@ -38,7 +38,12 @@ from hog.cpp_printing import (
 from hog.element_geometry import ElementGeometry
 from hog.exception import HOGException
 from hog.function_space import FunctionSpace, TrialSpace, TestSpace
-from hog.operator_generation.function_space_impls import FunctionSpaceImpl
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
+from hog.operator_generation.function_space_implementations.function_space_impl_factory import (
+    create_impl,
+)
 from hog.operator_generation.indexing import FaceType, CellType
 from hog.operator_generation.pystencils_extensions import create_generic_fields
 from hog.operator_generation.types import HOGPrecision, HOGType, hyteg_type
@@ -377,12 +382,8 @@ class ApplyWrapper(KernelWrapperType):
         dims: List[int] = [2, 3],
     ):
         self.name = "apply"
-        self.src: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
-            src_space, "src", type_descriptor
-        )
-        self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
-            dst_space, "dst", type_descriptor
-        )
+        self.src: FunctionSpaceImpl = create_impl(src_space, "src", type_descriptor)
+        self.dst: FunctionSpaceImpl = create_impl(dst_space, "dst", type_descriptor)
         self.src_fields = [self.src]
         self.dst_fields = [self.dst]
         self.dims = dims
@@ -562,7 +563,7 @@ class AssembleDiagonalWrapper(KernelWrapperType):
         dims: List[int] = [2, 3],
     ):
         self.name = "computeInverseDiagonalOperatorValues"
-        self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
+        self.dst: FunctionSpaceImpl = create_impl(
             fe_space, dst_field, type_descriptor, is_pointer=True
         )
         self.src_fields = []
@@ -690,12 +691,8 @@ class AssembleWrapper(KernelWrapperType):
     ):
         idx_t = HOGType("idx_t", np.int64)
         self.name = "toMatrix"
-        self.src: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
-            src_space, "src", idx_t
-        )
-        self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
-            dst_space, "dst", idx_t
-        )
+        self.src: FunctionSpaceImpl = create_impl(src_space, "src", idx_t)
+        self.dst: FunctionSpaceImpl = create_impl(dst_space, "dst", idx_t)
 
         # Treating both src and dst as dst_fields because src_fields are loaded
         # explicitly from memory prior to the kernel operation but we do not
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index b137074b28e9c55a3edc791822939857d48b306a..2e039cd10973ec92076d1ef3ea1488f02b43eebd 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -45,7 +45,12 @@ from hog.hyteg_code_generation import (
     GCC_WARNING_WORKAROUND,
 )
 from hog.logger import TimedLogger
-from hog.operator_generation.function_space_impls import FunctionSpaceImpl
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+    FunctionSpaceImpl,
+)
+from hog.operator_generation.function_space_implementations.function_space_impl_factory import (
+    create_impl,
+)
 from hog.operator_generation.pystencils_extensions import create_generic_fields
 
 import pystencils as ps
@@ -219,7 +224,6 @@ def micro_vertex_permutation_for_facet(
     """
 
     if volume_geometry == TriangleElement():
-
         if element_type == FaceType.BLUE:
             return [0, 1, 2]
 
@@ -232,7 +236,6 @@ def micro_vertex_permutation_for_facet(
         return shuffle_order_gray[facet_id]
 
     elif volume_geometry == TetrahedronElement():
-
         if element_type == CellType.WHITE_DOWN:
             return [0, 1, 2, 3]
 
@@ -557,7 +560,6 @@ class HyTeGElementwiseOperator:
         with TimedLogger(
             f"Generating kernels for operator {self.name}", level=logging.INFO
         ):
-
             # Generate each kernel type (apply, gemv, ...).
             self.generate_kernels()
 
@@ -602,7 +604,6 @@ class HyTeGElementwiseOperator:
         with TimedLogger("Generating C code from kernel AST(s)"):
             # Add all kernels to the class.
             for operator_method in self.operator_methods:
-
                 num_integrals = len(operator_method.integration_infos)
 
                 if num_integrals != len(
@@ -798,7 +799,6 @@ class HyTeGElementwiseOperator:
         )
         blending_includes = set()
         for dim, integration_infos in self.integration_infos.items():
-
             if not all(
                 [
                     integration_infos[0].blending.coupling_includes()
@@ -1194,7 +1194,6 @@ class HyTeGElementwiseOperator:
             element_types = list(integration_info.loop_strategy.element_loops.keys())
 
         for element_type in element_types:
-
             # Re-ordering micro-element vertices for the handling of domain boundary integrals.
             #
             # Boundary integrals are handled by looping over all (volume-)elements that have a facet at one of the
@@ -1283,7 +1282,7 @@ class HyTeGElementwiseOperator:
             coeffs = dict(
                 (
                     dof_symbol.function_id,
-                    FunctionSpaceImpl.create_impl(
+                    create_impl(
                         dof_symbol.function_space,
                         dof_symbol.function_id,
                         self._type_descriptor,
@@ -1530,7 +1529,6 @@ class HyTeGElementwiseOperator:
 
         for kernel_wrapper_type in self.kernel_wrapper_types:
             for dim, integration_infos in self.integration_infos.items():
-
                 kernel_functions = []
                 kernel_op_counts = []
                 platform_dep_kernels = []
@@ -1546,13 +1544,11 @@ class HyTeGElementwiseOperator:
                     )
 
                 for integration_info in integration_infos:
-
                     # generate AST of kernel loop
                     with TimedLogger(
                         f"Generating kernel {integration_info.name} ({kernel_wrapper_type.name}, {dim}D)",
                         logging.INFO,
                     ):
-
                         (
                             function_body,
                             kernel_op_count,
@@ -1655,7 +1651,6 @@ class HyTeGElementwiseOperator:
                 for kernel_function, integration_info in zip(
                     kernel_functions, integration_infos
                 ):
-
                     pre_call_code = ""
                     post_call_code = ""
 
@@ -1663,7 +1658,6 @@ class HyTeGElementwiseOperator:
                         integration_info.integration_domain
                         == MacroIntegrationDomain.DOMAIN_BOUNDARY
                     ):
-
                         if not isinstance(integration_info.loop_strategy, BOUNDARY):
                             raise HOGException(
                                 "The loop strategy should be BOUNDARY for boundary integrals."