implement inpainting

src/run.jl

@@ -6,7 +6,7 @@ using LinearAlgebra: norm
 
 function myrun()
     name = "test"
-    mesh = init_grid(200)
+    mesh = init_grid(50)
     d = 2
     m = 1
 
@@ -16,6 +16,7 @@ function myrun()
     Vp2 = FeSpace(mesh, DP1(), (m, d))
 
     g = FeFunction(Vg, name="g")
+    mask = FeFunction(Vg, name="mask")
     u = FeFunction(Vu, name="u")
     p1 = FeFunction(Vp1, name="p1")
     p2 = FeFunction(Vp2, name="p2")
@@ -41,18 +42,25 @@ function myrun()
     gamma2 = 1e-3
 
     interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
+    interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
+
+    #T(data, u) = u
+    #tdata = u # FIXME
+
+    @inline T(data, u) = iszero(data) ? zero(u) : u
+    tdata = mask
+
+    #T(data, u) = dot(data, u)
+    #tdata = nabla(fw)
 
-    # TODO: fixme
-    #T(u) = [u[1] + u[2]]
-    T(u) = u
     S(u) = u
 
-    function dp1_update(x; g, u, p1, du)
-	m1 = max(gamma1, norm(T(u) - g))
-	cond = norm(T(u) - g) > gamma1 ?
-	    dot(T(u) - g, T(du)) / norm(T(u) - g)^2 * p1 :
+    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 :
 	    zeros(size(p1))
-	return -p1 + alpha1 / m1 * (T(u) + T(du) - g) - cond
+	return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
     end
 
 
@@ -80,13 +88,13 @@ function myrun()
 	end
     end
 
-    @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2)
-	m1 = max(gamma1, norm(T(u) - g))
-	cond1 = norm(T(u) - g) > gamma1 ?
-	    dot(T(u) - g, T(du)) / norm(T(u) - g)^2 * p1 :
+    @inline 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 :
 	    zeros(size(p1))
-	a1 = alpha1 / m1 * dot(T(du), T(phi)) -
-	    dot(cond1, T(phi))
+	a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
+	    dot(cond1, T(tdata, phi))
 
 	m2 = max(gamma2, norm(nablau))
 	cond2 = norm(nablau) > gamma2 ?
@@ -95,20 +103,20 @@ function myrun()
 	a2 = lambda / m2 * dot(nabladu, nablaphi) -
 	    dot(cond2, nablaphi)
 
-	aB = alpha2 * dot(T(du), T(phi)) +
+	aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
 	    beta * dot(S(du), S(phi))
 
 	return a1 + a2 + aB
     end
 
-    @inline function du_l(x, phi, nablaphi; g, u, nablau)
-	aB = alpha2 * dot(T(u), T(phi)) +
+    @inline function du_l(x, phi, nablaphi; g, u, nablau, tdata)
+	aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
 	    beta * dot(S(u), S(phi))
-	m1 = max(gamma1, norm(T(u) - g))
-	p1part = alpha1 / m1 * dot(T(u) - g, T(phi))
+	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(phi))
+	gpart = alpha2 * dot(g, T(tdata, phi))
 
 	return -aB - p1part - p2part + gpart
     end
@@ -122,17 +130,17 @@ function myrun()
 
     pvd = paraview_collection("$(name).pvd")
     save_step!(pvd, 0)
-    for i = 1:5
-	print("newton step $i ...")
+    for i = 1:20
+	print("[$(lpad(i, 3, '0'))] ")
 	# solve du
-	A = assemble(Vu, du_a; g, u, nablau, p1, p2)
-	b = assemble_rhs(Vu, du_l; g, u, nablau)
-	print(" assembled ...")
+	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
-	println(" solved ...")
 
 	# solve dp1, dp2
-	interpolate!(dp1, dp1_update; g, u, p1, du)
+	interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
 	interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
 
 	# newton update
@@ -144,7 +152,9 @@ function myrun()
 	p1_project!()
 	p2_project!()
 
+	print("save ... ")
 	save_step!(pvd, i)
+	println()
     end
     vtk_save(pvd)
     return