implement l2 norm and introduce stopping criteria

src/function.jl

 using Statistics: mean
 using StaticArrays: SA, SArray, MVector, MMatrix
 export FeSpace, Mapper, FeFunction, P1, DP0, DP1
-export interpolate!, sample, bind!, evaluate, nabla
+export interpolate!, sample, bind!, evaluate, integrate, nabla
 
 # Finite Elements
 struct P1 end
@@ -182,7 +182,32 @@ function project!(dst::FeFunction, ::DP0, expr; params...)
 	gdofs = space.dofmap[:, 1, cell]
 	dst.data[gdofs] .= myvec(expr(centroid; opvalues...)) .* intel
     end
+end
+
+function integrate(mesh::Mesh, expr; params...)
+    opparams = NamedTuple(params)
+
+    qw, qx = quadrature()
 
+    # only for type inference
+    opvalues = map(f -> evaluate(f, SA[0., 0.]), opparams)
+    val = zero(expr(SA[0., 0.]; opvalues...))
+    for cell in cells(mesh)
+	foreach(f -> bind!(f, cell), opparams)
+
+	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+	intel = abs(det(delmap))
+
+	# quadrature points
+	for k in axes(qx, 2)
+	    xhat = qx[:, k]
+	    x = elmap(mesh, cell)(xhat)
+	    opvalues = map(f -> evaluate(f, x), opparams)
+
+	    val += qw[k] * expr(x; opvalues...) * intel
+	end
+    end
+    return val
 end
 
 function bind!(f::FeFunction, cell)

src/run.jl

@@ -323,6 +323,8 @@ function test_refine(mesh = nothing)
     return ctx
 end
 
+norm_l2(f) = sqrt(integrate(f.space.mesh, (x; f) -> dot(f, f); f))
+
 function denoise(img; name, params...)
     m = 1
     mesh = init_grid(img; type=:vertex)
@@ -334,7 +336,7 @@ function denoise(img; name, params...)
 
     interpolate!(ctx.g, x -> interpolate_bilinear(img, x))
     m = (size(img) .- 1) ./ 2 .+ 1
-    interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
+    interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3)
 
     save_denoise(ctx, i) =
 	save(ctx, "$(ctx.name)_$(lpad(i, 5, '0')).vtu",
@@ -345,20 +347,24 @@ function denoise(img; name, params...)
 
     k = 0
     while true
-	for i in 1:10
+	while true
 	    k += 1
 	    step!(ctx)
 	    estimate!(ctx)
 	    pvd[k] = save_denoise(ctx, k)
 	    println()
+	    println("ndofs: $(ndofs(ctx.u.space)), est: $(norm_l2(ctx.est)))")
+
+	    norm_step = sqrt(norm_l2(ctx.du)^2 + norm_l2(ctx.dp1)^2 + norm_l2(ctx.dp2))
+            norm_step <= 1e-3 && break
 	end
-	marked_cells = mark(ctx; theta = 0.7)
+	marked_cells = mark(ctx; theta = 0.5)
 	#println(marked_cells)
 	println("refining ...")
 	ctx = refine(ctx, marked_cells)
 	test_mesh(ctx.mesh)
-	interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
-	k >= 50 && break
+	interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3)
+	k >= 200 && break
     end
     vtk_save(pvd)
     return ctx
@@ -384,6 +390,8 @@ function inpaint(img, imgmask; name, params...)
     interpolate!(mask, x -> imgmask[round.(Int, x)...])
     #interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
     interpolate!(ctx.g, x -> imgmask[round.(Int, x)...] ? img[round.(Int, x)...] :  0.)
+    m = (size(img) .- 1) ./ 2 .+ 1
+    interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
 
     save_inpaint(i) =
 	save(ctx, "$(ctx.name)_$(lpad(i, 5, '0')).vtu",