diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 6b29bfa0fc17d4b90e2b74b7d6107dfce34b2281..cbae5239d62fb867e8bc936cffa3b9638ca7da5d 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -17,7 +17,7 @@ using StaticArrays: MVector, SArray, SMatrix, SVector, SA
 
 using SemiSmoothNewton
 using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs, diam,
-    elmap, assemble_lumped
+    elmap, assemble_lumped, norm_gradient, integrate_lumped
 using SemiSmoothNewton: project!, project_l2_lagrange!, project_qi_lagrange!,
     project_l2_pixel!
 using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
@@ -475,13 +475,33 @@ function step_pd2!(st::L1L2TVState; sigma, tau, theta = 1.)
 	tau * st.alpha2 * dot(st.T(tdata, w), st.T(tdata, phi)) +
 	tau * st.beta * dot(st.S(w, nablaw), st.S(phi, nablaphi))
 
+    u_au(x, w, nablaw, phi, nablaphi; g, u, p1, p2, tdata) =
+	dot(w, phi) +
+	tau * st.alpha2 * dot(st.T(tdata, w), st.T(tdata, phi)) +
+	tau * st.beta * dot(st.S(w, nablaw), st.S(phi, nablaphi))
+
+    u_ap(x, w, nablaw, phi, nablaphi; g, u, p1, p2, tdata) =
+        0.
+
     u_l(x, phi, nablaphi; g, u, p1, p2, tdata) =
 	dot(u, phi) - tau * (
 	    dot(p1, st.T(tdata, phi)) +
 	    dot(p2, nablaphi) -
 	    st.alpha2 * dot(g, st.T(tdata, phi)))
 
-    A, b = assemble(u_new.space, u_a, u_l; st.g, st.u, st.p1, st.p2, st.tdata)
+    u_lu(x, phi, nablaphi; g, u, p1, p2, tdata) =
+	dot(u, phi) - tau * (-
+	    st.alpha2 * dot(g, st.T(tdata, phi)))
+
+    u_lp(x, phi, nablaphi; g, u, p1, p2, tdata) =
+	- tau * (
+	    dot(p1, st.T(tdata, phi)) +
+	    dot(p2, nablaphi))
+
+    Au, bu = assemble(u_new.space, u_au, u_lu; st.g, st.u, st.p1, st.p2, st.tdata)
+    Ap, bp = assemble_lumped(u_new.space, u_ap, u_lp; st.g, st.u, st.p1, st.p2, st.tdata)
+    A = Au + Ap
+    b = bu + bp
     u_new.data .= A \ b
     st.du.data .= u_new.data .- st.u.data
     st.u.data .= u_new.data
@@ -507,18 +527,24 @@ function step_d!(st::L1L2TVState; tau)
 end
 
 function solve_primal!(u::FeFunction, st::L1L2TVState)
-    u_a(x, u, nablau, phi, nablaphi; g, p1, p2, tdata) =
+    u_au(x, u, nablau, phi, nablaphi; g, p1, p2, tdata) =
 	st.alpha2 * dot(st.T(tdata, u), st.T(tdata, phi)) +
 	    st.beta * dot(st.S(u, nablau), st.S(phi, nablaphi))
+    u_ap(x, u, nablau, phi, nablaphi; g, p1, p2, tdata) =
+        0.
 
-    u_l(x, phi, nablaphi; g, p1, p2, tdata) =
-	-dot(p1, st.T(tdata, phi)) - dot(p2, nablaphi) +
-	    st.alpha2 * dot(g, st.T(tdata, phi))
+    u_lu(x, phi, nablaphi; g, p1, p2, tdata) =
+        st.alpha2 * dot(g, st.T(tdata, phi))
+    u_lp(x, phi, nablaphi; g, p1, p2, tdata) =
+	-dot(p1, st.T(tdata, phi)) - dot(p2, nablaphi)
 
     # u = -B^{-1} * (T^* p_1 + nabla^* p_2 - alpha2 * T^* g)
     # solve:
     # <B u, phi> = -<p_1, T phi> - <p_2, nablaphi> + alpha_2 * <g, T phi>
-    A, b = assemble(u.space, u_a, u_l; st.g, st.p1, st.p2, st.tdata)
+    Au, bu = assemble(u.space, u_au, u_lu; st.g, st.p1, st.p2, st.tdata)
+    Ap, bp = assemble_lumped(u.space, u_ap, u_lp; st.g, st.p1, st.p2, st.tdata)
+    A = Au + Ap
+    b = bu + bp
     u.data .= A \ b
 end
 
@@ -754,7 +780,7 @@ function norm_residual(st::L1L2TVState)
 	    st.lambda * nablau
 	return norm(p1part)^2 + norm(p2part)^2
     end
-    ppart2 = integrate(st.mesh, integrand; st.g, st.u,
+    ppart2 = integrate_lumped(st.mesh, integrand; st.g, st.u,
         nablau = nabla(st.u), st.p1, st.p2, st.tdata)
 
     return sqrt(upart2 + ppart2)
diff --git a/src/function.jl b/src/function.jl
index cd03e9eba90e45b5121664ac72a66bdecbb25fc4..98d56234c9882286eb1d519f88035a75153ffa19 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -193,11 +193,18 @@ function project!(dst::FeFunction, ::DP0, expr; params...)
     end
 end
 
-# TODO: interface is lacking, which quadrature should be used?
-function integrate(mesh::Mesh, expr; params...)
+# TODO: better interface for quadrature choice
+
+integrate(mesh::Mesh, expr; params...) =
+    integrate(mesh::Mesh, quadrature(), expr; params...)
+
+integrate_lumped(mesh::Mesh, expr; params...) =
+    integrate(mesh::Mesh, quadrature_lumped(), expr; params...)
+
+function integrate(mesh::Mesh, quadrature, expr; params...)
     opparams = NamedTuple(params)
 
-    qw, qx = quadrature()
+    qw, qx = quadrature
 
     # only for type inference
     opvalues = map(f -> evaluate(f, SA[0., 0.]), opparams)