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

remove old run function

parent 8753ec54
Branches
Tags
No related merge requests found
......@@ -264,173 +264,3 @@ function optflow(imgf0, imgf1; name, params...)
println()
end
end
function myrun()
name = "test"
mesh = init_grid(50)
d = 2
m = 2
Vg = FeSpace(mesh, P1(), (1,))
Vu = FeSpace(mesh, P1(), (m,))
Vp1 = FeSpace(mesh, DP0(), (1,))
Vp2 = FeSpace(mesh, DP1(), (m, d))
# inpainting
mask = FeFunction(Vg, name="mask")
g = FeFunction(Vg, name="g")
u = FeFunction(Vu, name="u")
p1 = FeFunction(Vp1, name="p1")
p2 = FeFunction(Vp2, name="p2")
du = FeFunction(Vu)
dp1 = FeFunction(Vp1)
dp2 = FeFunction(Vp2)
nablau = nabla(u)
nabladu = nabla(du)
g.data .= 0
u.data .= 0
p1.data .= 0
p2.data .= 0
du.data .= 0
dp1.data .= 0
dp2.data .= 0
alpha1 = 0.
alpha2 = 10.
beta = 1e-3
lambda = 0.01
gamma1 = 1e-3
gamma2 = 1e-3
interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
interpolate!(f0, x -> x[1])
interpolate!(f1, x -> x[1] - 0.01)
fw.data .= f1.data
g_optflow(x; u, f0, fw, nablafw) =
nablafw * u - (fw - f0)
interpolate!(g, g_optflow; u, f0, fw, nablafw)
#T(data, u) = u
#tdata = u # FIXME
#T(data, u) = iszero(data) ? zero(u) : u
#tdata = mask
T(data, u) = data * u
tdata = nabla(fw)
#S(u, nablau) = u
S(u, nablau) = nablau
function dp1_update(x; g, u, p1, du, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zero(p1)
return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
end
function dp2_update(x; u, nablau, p2, du, nabladu)
m2 = max(gamma2, norm(nablau))
cond = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zero(p2)
return -p2 + lambda / m2 * (nablau + nabladu) - cond
end
function p1_project!()
p1.space.element::DP0
p1.data .= clamp.(p1.data, -alpha1, alpha1)
end
function p2_project!()
p2.space.element::DP1
p2d = reshape(p2.data, m * d * 3, :) # no copy
for i in axes(p2d, 2)
p2in = norm(p2d[:, i])
if p2in > lambda
p2d[:, i] .*= lambda ./ p2in
end
end
end
function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond1 = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zero(p1)
a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
dot(cond1, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
cond2 = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zero(p2)
a2 = lambda / m2 * dot(nabladu, nablaphi) -
dot(cond2, nablaphi)
aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
beta * dot(S(du, nabladu), S(phi, nablaphi))
return a1 + a2 + aB
end
function du_l(x, phi, nablaphi; g, u, nablau, tdata)
aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
beta * dot(S(u, nablau), S(phi, nablaphi))
m1 = max(gamma1, norm(T(tdata, u) - g))
p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
p2part = lambda / m2 * dot(nablau, nablaphi)
gpart = alpha2 * dot(g, T(tdata, phi))
return -aB - p1part - p2part + gpart
end
function save_step!(pvd, i)
vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh)
vtk_append!(vtk, g, u, p1, p2)
vtk_save(vtk)
pvd[i] = vtk
end
pvd = paraview_collection("$(name).pvd")
save_step!(pvd, 0)
for i = 1:20
print("[$(lpad(i, 3, '0'))] ")
# solve du
print("assemble ... ")
A = assemble(Vu, du_a; g, u, nablau, p1, p2, tdata)
b = assemble_rhs(Vu, du_l; g, u, nablau, tdata)
print("solve ... ")
du.data .= A \ b
# solve dp1, dp2
interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
# newton update
u.data .+= du.data
p1.data .+= dp1.data
p2.data .+= dp2.data
# reproject p1, p2
p1_project!()
p2_project!()
print("save ... ")
save_step!(pvd, i)
println()
end
vtk_save(pvd)
return
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment