diff --git a/generate_all_operators.py b/generate_all_operators.py
index 005b4d30b2a125d42503e1ad733715da8b819e87..c9b2050b5d1e55a7deaa62c36d7d8ab3e920f1ed 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -43,6 +43,7 @@ from hog.forms import (
     divergence,
     gradient,
     div_k_grad,
+    k_grad,
     shear_heating,
     epsilon,
     full_stokes,
@@ -633,6 +634,10 @@ def all_operators(
 
     # fmt: on
 
+    ops.append(OperatorInfo("P2ToP2Vector", "KGrad", TrialSpace(P2), TestSpace(P2Vector),
+                            form=partial(k_grad, coefficient_function_space=P2),
+                            type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
     p2vec_epsilon = partial(
         epsilon,
         variable_viscosity=True,
diff --git a/hog/forms.py b/hog/forms.py
index b959b215a52d2f2f5c48f1d5806e1f221a64510f..7d61d99dd56454dfd99439bebffaba88ca6634c8 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -168,6 +168,47 @@ Weak formulation
         docstring=docstring,
     )
 
+def k_grad(
+    trial: TrialSpace,
+    test: TestSpace,
+    geometry: ElementGeometry,
+    symbolizer: Symbolizer,
+    blending: GeometryMap = IdentityMap(),
+    coefficient_function_space: Optional[FunctionSpace] = None,
+) -> Form:
+    docstring = f"""
+Diffusion operator with a scalar coefficient.
+
+Geometry map: {blending}
+
+Weak formulation
+
+    u: trial function (space: {trial})
+    v: test function  (space: {test})
+    k: coefficient    (space: {coefficient_function_space})
+
+    ∫ k v ∇u
+"""
+
+    from hog.recipes.integrands.volume.k_grad import integrand, integrand_affine
+
+    # mypy type checking supposedly cannot figure out ternaries.
+    # https://stackoverflow.com/a/70832391
+    integr: Callable[..., Any] = integrand
+    if blending == IdentityMap():
+        integr = integrand_affine
+
+    return process_integrand(
+        integr,
+        trial,
+        test,
+        geometry,
+        symbolizer,
+        blending=blending,
+        fe_coefficients={"k": coefficient_function_space},
+        is_symmetric=trial == test,
+        docstring=docstring,
+    )
 
 def nonlinear_diffusion(
     trial: TrialSpace,
diff --git a/hog/recipes/integrands/volume/k_grad.py b/hog/recipes/integrands/volume/k_grad.py
new file mode 100644
index 0000000000000000000000000000000000000000..5780288315a804f24aff1a129480ecfaf40a12a3
--- /dev/null
+++ b/hog/recipes/integrands/volume/k_grad.py
@@ -0,0 +1,51 @@
+# 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(
+    *,
+    v,
+    jac_a_inv,
+    jac_a_abs_det,
+    jac_b_inv,
+    jac_b_abs_det,
+    grad_u,
+    k,
+    tabulate,
+    **_,
+):
+    return k["k"] * (
+            dot((jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)), v)
+        * tabulate(jac_a_abs_det)
+        * jac_b_abs_det
+    )
+
+def integrand_affine(
+    *,
+    v,
+    jac_a_inv,
+    jac_a_abs_det,
+    grad_u,
+    k,
+    tabulate,
+    **_,
+):
+    return k["k"] * (
+            dot((tabulate(jac_a_inv.T * grad_u)), v)
+        * tabulate(jac_a_abs_det)
+    )