diff --git a/python/dune/codegen/sumfact/geometry.py b/python/dune/codegen/sumfact/geometry.py
index cef07354adb7dbfbbb85d42e72fc5e7e14d57a08..3a9a5a5ad9a6ea240bfc4ffc2534f610e12eb97d 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/python/dune/codegen/sumfact/geometry.py
@@ -56,6 +56,9 @@ import loopy as lp
 
 @geometry_mixin("sumfact_multilinear")
 class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
+    def nonsumfact_fallback(self):
+        return "generic"
+
     def _jacobian_inverse(self):
         return "jit"
 
@@ -243,6 +246,9 @@ class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
 
 @geometry_mixin("sumfact_axiparallel")
 class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
+    def nonsumfact_fallback(self):
+        return "axiparallel"
+
     def facet_normal(self, o):
         i, = self.indices
         self.indices = None
@@ -260,6 +266,9 @@ class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
 
 @geometry_mixin("sumfact_equidistant")
 class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParallelGeometryMixin):
+    def nonsumfact_fallback(self):
+        return "equidistant"
+
     def facet_jacobian_determinant(self, o):
         name = "fdetjac"
         self.define_facet_jacobian_determinant(name)
diff --git a/python/dune/codegen/sumfact/switch.py b/python/dune/codegen/sumfact/switch.py
index 67cb9df358112e456b4ef9250c3602f1f9aacdec..a031790dd5321841eb8c100195cde675651fe1d3 100644
--- a/python/dune/codegen/sumfact/switch.py
+++ b/python/dune/codegen/sumfact/switch.py
@@ -2,7 +2,8 @@
 
 import csv
 
-from dune.codegen.generation import (get_global_context_value,
+from dune.codegen.generation import (construct_from_mixins,
+                                     get_global_context_value,
                                      global_context,
                                      )
 from dune.codegen.pdelab.geometry import world_dimension
@@ -26,8 +27,18 @@ def sumfact_generate_kernels_per_integral(integrals):
         # Maybe skip sum factorization on boundary integrals
         if not get_form_option("sumfact_on_boundary"):
             set_form_option("sumfact", False)
+
+            # Try to find a fallback for sum factorized geometry mixins
+            geometry_backup = get_form_option("geometry_mixins")
+            mixin = construct_from_mixins(mixins=[geometry_backup])()
+            if hasattr(mixin, "nonsumfact_fallback"):
+                set_form_option("geometry_mixins", mixin.nonsumfact_fallback())
+
             for k in generate_kernels_per_integral(integrals):
                 yield k
+
+            # Reset state
+            set_form_option("geometry_mixins", geometry_backup)
             set_form_option("sumfact", True)
             return