diff --git a/generate_all_operators.py b/generate_all_operators.py
index e92c28af80f9a8b5dac97ef16451d7430cfe3b79..bb20e9a895a33f24d1511102804dc78758b9bcdc 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -53,7 +53,7 @@ from hog.forms import (
     supg_advection,
     grad_rho_by_rho_dot_u,
     nonlinear_diffusion_stiffness,
-    nonlinear_diffusion_newton_galerkin_approximated_mass_lumped
+    nonlinear_diffusion_newton_galerkin_approximated_mass_lumped, nonlinear_diffusion_mass_lumping_template
 )
 from hog.forms_boundary import (
     mass_boundary,
@@ -653,6 +653,12 @@ def all_operators(
                             coefficient_function_space=P1, use_a_derivative=True),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
+    ops.append(OperatorInfo("P1", "NonlinearDiffusionMassLumpingTemplate", TrialSpace(P1),
+                            TestSpace(P1), form=partial(nonlinear_diffusion_mass_lumping_template,
+                            coefficient_function_space=P1),
+                            type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
+
     ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinApproximatedMassLumped", TrialSpace(P1),
                             TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin_approximated_mass_lumped,
                                                         coefficient_function_space=P1),
diff --git a/hog/forms.py b/hog/forms.py
index b3338d8c6b98736c66ba8d97cef869368b8ddd6e..357f0e2987ba5f775a9b2972bf40240cac15aca9 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -383,6 +383,72 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
         fe_coefficients={"k": coefficient_function_space},
     )
 
+
+def nonlinear_diffusion_mass_lumping_template(
+        trial: TrialSpace,
+        test: TestSpace,
+        geometry: ElementGeometry,
+        symbolizer: Symbolizer,
+        coefficient_function_space: FunctionSpace,
+        blending: GeometryMap = IdentityMap(),
+) -> Form:
+    docstring = f"""
+todo
+
+Weak formulation
+
+    u: trial function (space: {trial})
+    v: test function  (space: {test})
+    k: FE coefficient function (space: {coefficient_function_space})
+
+    ∫ a(k) ∇k · ∇v
+
+    or ∫ a'(k) ∇u · ∇v if use_a_derivative was set
+
+Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended for :math:`k = u`.
+"""
+
+    if trial != test:
+        raise HOGException(
+            "Trial space must be equal to test space to assemble diffusion matrix."
+        )
+
+    if blending != IdentityMap():
+        raise HOGException(
+            "The nonlinear_diffusion_newton_galerkin form does currently not support blending."
+        )
+
+    def integrand(
+            *,
+            jac_a_inv,
+            jac_a_abs_det,
+            grad_k,
+            grad_v,
+            k,
+            tabulate,
+            **_,
+    ):
+        a_prime = 2 * k["k"]
+
+
+        term = a_prime * tabulate(
+            dot(jac_a_inv.T * grad_k["k"], jac_a_inv.T * grad_v) * jac_a_abs_det
+        )
+
+        return term
+
+    return process_integrand(
+        integrand,
+        trial,
+        test,
+        geometry,
+        symbolizer,
+        blending=blending,
+        is_symmetric=False,
+        docstring=docstring,
+        fe_coefficients={"k": coefficient_function_space},
+    )
+
 def nonlinear_diffusion_newton_galerkin_approximated_mass_lumped(
         trial: TrialSpace,
         test: TestSpace,