diff --git a/hog/forms.py b/hog/forms.py
index 8ae046393a7a8ad2b0330b12c960bb6eccfa7c1b..9e79d6a4f1247b74b2e532a2474b2ce5aa7d2205 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -321,6 +321,8 @@ def epsilon(
     component_test: int = 0,
     variable_viscosity: bool = True,
     coefficient_function_space: Optional[FunctionSpace] = None,
+    rotation_wrapper: bool = False,
+    rigid_body_mode: bool = False,
 ) -> Form:
     docstring = f"""
 "Epsilon" operator.
@@ -343,6 +345,9 @@ where
         raise HOGException("Constant viscosity currently not supported.")
 
     from hog.recipes.integrands.volume.epsilon import integrand
+    from hog.recipes.integrands.volume.epsilon_rigid import (
+        integrand as integrand_rigid,
+    )
     from hog.recipes.integrands.volume.epsilon_affine import (
         integrand as integrand_affine,
     )
@@ -352,6 +357,8 @@ where
     integr: Callable[..., Any] = integrand
     if blending == IdentityMap():
         integr = integrand_affine
+    elif rigid_body_mode:
+        integr = integrand_rigid
 
     return process_integrand(
         integr,
@@ -939,6 +946,49 @@ Weak formulation
         },
     )
 
+def supg_mass(
+    trial: TrialSpace,
+    test: TestSpace,
+    geometry: ElementGeometry,
+    symbolizer: Symbolizer,
+    velocity_function_space: FunctionSpace,
+    delta_function_space: FunctionSpace,
+    blending: GeometryMap = IdentityMap(),
+) -> Form:
+    docstring = f"""
+advection operator which needs to be used in combination with SUPG
+
+Geometry map:    {blending}
+
+Weak formulation
+
+    u: trial function (scalar space: {trial})
+    v: test function  (scalar space: {test})
+    𝛿: delta function (vectorial space: {delta_function_space})
+    w: velocity function (vectorial space: {velocity_function_space})
+
+    ∫ u 𝛿(w · ∇v)
+"""
+
+    from hog.recipes.integrands.volume.supg_mass import integrand
+
+    return process_integrand(
+        integrand,
+        trial,
+        test,
+        geometry,
+        symbolizer,
+        blending=blending,
+        is_symmetric=False,
+        docstring=docstring,
+        fe_coefficients={
+            "ux": velocity_function_space,
+            "uy": velocity_function_space,
+            "uz": velocity_function_space,
+            "delta": delta_function_space,
+        },
+    )
+
 
 def supg_diffusion(
     trial: TrialSpace,
diff --git a/hog/recipes/integrands/volume/supg_mass.py b/hog/recipes/integrands/volume/supg_mass.py
new file mode 100644
index 0000000000000000000000000000000000000000..91efcc60df4228dd88fd2b0f02190f70fbd55a81
--- /dev/null
+++ b/hog/recipes/integrands/volume/supg_mass.py
@@ -0,0 +1,44 @@
+# 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.recipes.common import *
+
+
+def integrand(
+    *,
+    u,
+    jac_a_inv,
+    jac_a_abs_det,
+    jac_b_inv,
+    jac_b_abs_det,
+    grad_v,
+    k,
+    volume_geometry,
+    tabulate,
+    **_,
+):
+    if volume_geometry.dimensions > 2:
+        w = sp.Matrix([[k["ux"]], [k["uy"]], [k["uz"]]])
+    else:
+        w = sp.Matrix([[k["ux"]], [k["uy"]]])
+
+    return (
+        k["delta"] 
+        * u
+        * dot(jac_b_inv.T * tabulate(jac_a_inv.T * grad_v), w)
+        * tabulate(jac_a_abs_det)
+        * jac_b_abs_det
+    )
diff --git a/hog_tests/operator_generation/test_indexing.py b/hog_tests/operator_generation/test_indexing.py
index 1387e7506ac20bfa5e2d006894835e44f992e7dc..da31d156a8f71b06512c5a3d93e2da3963c30eb1 100644
--- a/hog_tests/operator_generation/test_indexing.py
+++ b/hog_tests/operator_generation/test_indexing.py
@@ -277,7 +277,6 @@ def test_micro_volume_to_volume_indices():
         array_index = sp.simplify(
             dof_indices[intra_primitive_index].array_index(geometry, indexing_info)
         )
-        print(array_index)
         assert array_index == target_array_index
 
     # 2D, P0: