From 754a746b9608ba108bccd2265511971ee5c02519 Mon Sep 17 00:00:00 2001 From: Stephan Hilb <stephan@ecshi.net> Date: Thu, 15 Jul 2021 13:24:02 +0200 Subject: [PATCH] implement l2 norm and introduce stopping criteria --- src/function.jl | 27 ++++++++++++++++++++++++++- src/run.jl | 18 +++++++++++++----- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/function.jl b/src/function.jl index 3c1e325..b4acb7a 100644 --- a/src/function.jl +++ b/src/function.jl @@ -1,7 +1,7 @@ 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) diff --git a/src/run.jl b/src/run.jl index e51e2d6..39639fc 100644 --- a/src/run.jl +++ b/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", -- GitLab