Skip to content
Snippets Groups Projects
Commit fbf2276a authored by Stephan Hilb's avatar Stephan Hilb
Browse files

implement mesh area

parent 59c92154
No related branches found
No related tags found
No related merge requests found
...@@ -340,6 +340,14 @@ function save(filename::String, mesh::Mesh, fs...) ...@@ -340,6 +340,14 @@ function save(filename::String, mesh::Mesh, fs...)
vtk_save(vtk) vtk_save(vtk)
end end
function area(mesh, cell)
A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
view(mesh.vertices, :, view(mesh.cells, :, cell)))
return det(A[:, 2:3] .- A[:, 1]) / 2
end
area(mesh) = mapreduce(cell -> area(mesh, cell), +, cells(mesh))
function diam(mesh, cell) function diam(mesh, cell)
A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
view(mesh.vertices, :, view(mesh.cells, :, cell))) view(mesh.vertices, :, view(mesh.cells, :, cell)))
......
...@@ -310,12 +310,14 @@ function denoise(img; name, params...) ...@@ -310,12 +310,14 @@ function denoise(img; name, params...)
estimate!(ctx) estimate!(ctx)
pvd[k] = save_denoise(ctx, k) pvd[k] = save_denoise(ctx, k)
println() println()
norm_step = sqrt((norm_l2(ctx.du)^2 + norm_l2(ctx.dp1)^2 + norm_l2(ctx.dp2)^2) / area(mesh))
println("ndofs: $(ndofs(ctx.u.space)), est: $(norm_l2(ctx.est)))") println("ndofs: $(ndofs(ctx.u.space)), est: $(norm_l2(ctx.est)))")
println("primal energy: $(primal_energy(ctx))") println("primal energy: $(primal_energy(ctx))")
println("norm_step: $(norm_step)")
norm_step = sqrt(norm_l2(ctx.du)^2 + norm_l2(ctx.dp1)^2 + norm_l2(ctx.dp2)) norm_step <= 1e-1 && break
norm_step = sqrt(norm_l2(ctx.du)^2)
norm_step <= 1e-2 && break
end end
marked_cells = mark(ctx; theta = 0.5) marked_cells = mark(ctx; theta = 0.5)
#println(marked_cells) #println(marked_cells)
...@@ -327,7 +329,7 @@ function denoise(img; name, params...) ...@@ -327,7 +329,7 @@ function denoise(img; name, params...)
ctx.g.data .= gnew.data ctx.g.data .= gnew.data
#interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3) #interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3)
k >= 50 && break k >= 100 && break
end end
vtk_save(pvd) vtk_save(pvd)
return ctx return ctx
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment