diff --git a/generate_all_operators.py b/generate_all_operators.py
index a9258e518a969bae128cbfc23c35609bca198ec2..34696da3c4a7220857b73ec6c8bf6d5c3b27b908 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -47,7 +47,6 @@ from hog.element_geometry import (
 from hog.exception import HOGException
 from hog.forms import (
     diffusion,
-    rigid_body_mode_angular,
     divergence,
     gradient,
     div_k_grad,
@@ -61,6 +60,7 @@ from hog.forms import (
     advection,
     supg_advection,
     grad_rho_by_rho_dot_u,
+    rigid_body_mode_angular,
 )
 from hog.forms_boundary import (
     mass_boundary,
@@ -666,7 +666,12 @@ def all_operators(
     ops.append(OperatorInfo("P2PlusBubble", "Diffusion", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=diffusion,
                             type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
 
-    ops.append(OperatorInfo("P2Vector", "RigidModeAngular", TrialSpace(P2Vector), TestSpace(P2Vector), form=rigid_body_mode_angular,
+    ops.append(OperatorInfo("P2Vector", "RigidModeAngular", TrialSpace(P2Vector), TestSpace(P2Vector), 
+                            form=rigid_body_mode_angular,
+                            type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+    
+    ops.append(OperatorInfo("P2VectorToP1", "DivergenceCompressible", TrialSpace(P2Vector), TestSpace(P1), 
+                            form=partial(divergence, compressible=True, density_function_space=P1),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
     # fmt: on
 
diff --git a/hog/forms.py b/hog/forms.py
index 9688729be116c37e1fa46fcbc7588a630b74b9a6..e16c47292c50bba4e905047d61f4570bcc7bda04 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -638,6 +638,8 @@ def divergence(
     blending: GeometryMap = IdentityMap(),
     component_index: int = 0,
     rotation_wrapper: bool = False,
+    compressible: bool = False,
+    density_function_space: Optional[FunctionSpace] = None,
 ) -> Form:
     docstring = f"""
 Divergence.
@@ -651,25 +653,49 @@ Weak formulation
     v: test function  (scalar space:    {test})
 
     ∫ - ( ∇ · u ) v
+    
+if compressible is True
+    ∫ - ( ∇ · ( 𝜌u ) ) v
+    
 """
 
     from hog.recipes.integrands.volume.divergence import integrand
+    from hog.recipes.integrands.volume.divergence import integrand_compressible
     from hog.recipes.integrands.volume.rotation import RotationType
 
-    return process_integrand(
-        integrand,
-        trial,
-        test,
-        geometry,
-        symbolizer,
-        blending=blending,
-        component_index=component_index,
-        is_symmetric=False,
-        docstring=docstring,
-        rot_type=RotationType.POST_MULTIPLY
-        if rotation_wrapper
-        else RotationType.NO_ROTATION,
-    )
+    if compressible:
+        if not trial.is_vectorial:
+            raise HOGException("Compressible version of divergence form only supports vectorial functions")
+        
+        return process_integrand(
+            integrand_compressible,
+            trial,
+            test,
+            geometry,
+            symbolizer,
+            blending=blending,
+            is_symmetric=False,
+            docstring=docstring,
+            fe_coefficients={"rho": density_function_space},
+            rot_type=RotationType.POST_MULTIPLY
+            if rotation_wrapper
+            else RotationType.NO_ROTATION,
+        )
+    else:
+        return process_integrand(
+            integrand,
+            trial,
+            test,
+            geometry,
+            symbolizer,
+            blending=blending,
+            component_index=component_index,
+            is_symmetric=False,
+            docstring=docstring,
+            rot_type=RotationType.POST_MULTIPLY
+            if rotation_wrapper
+            else RotationType.NO_ROTATION,
+        )
 
 def gradient(
     trial: TrialSpace,
diff --git a/hog/recipes/integrands/volume/divergence.py b/hog/recipes/integrands/volume/divergence.py
index f2ea6530cb567f2b2874f35582f080692884e18e..8ac453cd8000b2124d7057d48b1266799d24aa36 100644
--- a/hog/recipes/integrands/volume/divergence.py
+++ b/hog/recipes/integrands/volume/divergence.py
@@ -44,3 +44,29 @@ def integrand(
             * tabulate(v * jac_a_abs_det)
             * jac_b_abs_det
         )
+
+def integrand_compressible(
+    *,
+    jac_a_inv,
+    jac_a_abs_det,
+    jac_b_inv,
+    jac_b_abs_det,
+    u,
+    grad_u,
+    v,
+    k,
+    grad_k,
+    tabulate,
+    **_,
+):
+    rho = k["rho"]
+    grad_rho = jac_b_inv.T * tabulate(jac_a_inv.T * grad_k["rho"])
+
+    div_u = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)).trace()
+    div_rho_u = rho * div_u + dot(grad_rho, u)[0]
+
+    return (
+        (-div_rho_u)
+        * tabulate(v * jac_a_abs_det)
+        * jac_b_abs_det
+    )