diff --git a/src/SemiSmoothNewton.jl b/src/SemiSmoothNewton.jl
index 47b41372fbecc66049c15a2184f0bfc212fa6425..eae1c6550f9cb6f6be892a1b22f6253291577616 100644
--- a/src/SemiSmoothNewton.jl
+++ b/src/SemiSmoothNewton.jl
@@ -5,6 +5,8 @@ include("mesh.jl")
include("function.jl")
include("operator.jl")
+include("run.jl")
+
function clip_line(r0, r1, l0, l1)
d_l = -(l1[1] - l0[1])
d_r = +(l1[1] - l0[1])
diff --git a/src/function.jl b/src/function.jl
index fdfa40de53f4480be323e504c60ad11ff7d09178..8a9295cdd1cca866f613ab947f3163a631f862b9 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -1,6 +1,6 @@
using Statistics: mean
export FeSpace, Mapper, FeFunction, P1, DP0, DP1
-export interpolate!, sample, bind!, evaluate
+export interpolate!, sample, bind!, evaluate, nabla
# Finite Elements
struct P1 end
@@ -10,7 +10,7 @@ struct DP1 end
# FIXME: should be static vectors
# evaluate all (1-dim) local basis functions against x
evaluate_basis(::P1, x) = [1 - x[1] - x[2], x[1], x[2]]
-evaluate_basis(::DP0, x) = [x]
+evaluate_basis(::DP0, x) = [1]
evaluate_basis(::DP1, x) = evaluate_basis(P1(), x)
# non-mixed
@@ -36,7 +36,7 @@ function FeSpace(mesh, el::P1, size_=(1,))
end
function FeSpace(mesh, el::DP0, size_=(1,))
- # special case for P1: dofs correspond to vertices
+ # special case for DP1: dofs correspond to cells
rdims = prod(size_)
ncells = size(mesh.cells, 2)
dofmap = Array{Int, 3}(undef, rdims, 1, ncells)
@@ -46,6 +46,17 @@ function FeSpace(mesh, el::DP0, size_=(1,))
return FeSpace(mesh, el, dofmap, rdims * ncells, size_)
end
+function FeSpace(mesh, el::DP1, size_=(1,))
+ # special case for DP1: dofs correspond to vertices
+ rdims = prod(size_)
+ ncells = size(mesh.cells, 2)
+ dofmap = Array{Int, 3}(undef, rdims, 3, ncells)
+ for i in LinearIndices((rdims, 3, ncells))
+ dofmap[i] = i
+ end
+ return FeSpace(mesh, el, dofmap, rdims * 3 * ncells, size_)
+end
+
Base.show(io::IO, ::MIME"text/plain", x::FeSpace) =
print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(x.ndofs) dofs")
@@ -78,38 +89,54 @@ end
Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
print("$(nameof(typeof(f))), size $(f.space.size) with $(length(f.data)) dofs")
-interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.space.element, src)
+interpolate!(dst::FeFunction, expr::Function; params...) = interpolate!(dst, dst.space.element, expr; params...)
myvec(x) = vec(x)
myvec(x::Number) = x
-function interpolate!(dst::FeFunction, ::P1, src::Function)
+function interpolate!(dst::FeFunction, ::P1, expr::Function; params...)
+ params = NamedTuple(params)
space = dst.space
mesh = space.mesh
for cell in axes(mesh.cells, 2)
+ for f in params
+ bind!(f, cell)
+ end
for eldof in axes(mesh.cells, 1)
xid = mesh.cells[eldof, cell]
x = mesh.vertices[:, xid]
+ xloc = [0. 1. 0.; 0. 0. 1.][:, eldof]
+
+ opvalues = map(f -> evaluate(f, xloc), params)
gdofs = space.dofmap[:, eldof, cell]
- dst.data[gdofs] .= myvec(src(x))
+ dst.data[gdofs] .= myvec(expr(x; opvalues...))
end
end
end
-function interpolate!(dst::FeFunction, ::DP0, src::Function)
+function interpolate!(dst::FeFunction, ::DP0, expr::Function; params...)
+ params = NamedTuple(params)
space = dst.space
mesh = space.mesh
for cell in axes(mesh.cells, 2)
+ for f in params
+ bind!(f, cell)
+ end
vertices = mesh.vertices[:, mesh.cells[:, cell]]
centroid = mean(vertices, dims = 2)
+ lcentroid = [1/3, 1/3]
+
+ opvalues = map(f -> evaluate(f, lcentroid), params)
gdofs = space.dofmap[:, 1, cell]
- dst.data[gdofs] .= myvec(src(centroid))
+ dst.data[gdofs] .= myvec(expr(centroid; opvalues...))
end
end
+interpolate!(dst::FeFunction, ::DP1, expr::Function; params...) =
+ interpolate!(dst, P1(), expr; params...)
function bind!(f::FeFunction, cell)
f.ldata .= vec(f.data[f.space.dofmap[:, :, cell]])
@@ -119,11 +146,22 @@ end
# evaluate at local point (needs bind! call before)
evaluate(f::FeFunction, x) = evaluate(f.space, f.ldata, x)
-struct CellFunction
- f::FeFunction
- dofs::Vector{Float64}
+struct Derivative{F}
+ f::F
+end
+
+Base.show(io::IO, ::MIME"text/plain", x::Derivative) =
+ print("$(nameof(typeof(x))) of $(typeof(x.f))")
+
+nabla(f) = Derivative(f)
+
+bind!(df::Derivative, cell) = bind!(df.f, cell)
+function evaluate(df::Derivative, x)
+ jac = jacobian(x -> evaluate(df.f.space, df.f.ldata, x), x)
+ return reshape(jac, df.f.space.size..., :)
end
+
function sample(f::FeFunction)
mesh = f.mapper.mesh
for cell in axes(mesh.cells, 2)
@@ -143,6 +181,7 @@ function sample(f::FeFunction)
end
+
append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element)
function append_data!(vtkfile, f::FeFunction, ::P1)
diff --git a/src/run.jl b/src/run.jl
new file mode 100644
index 0000000000000000000000000000000000000000..54f6e0b5272cdd5dee746e85b034a10dbf86e57d
--- /dev/null
+++ b/src/run.jl
@@ -0,0 +1,66 @@
+export myrun
+
+using LinearAlgebra: norm
+
+
+
+function myrun()
+ mesh = init_grid(50)
+ d = 2
+ m = 1
+
+ Vg = FeSpace(mesh, DP0(), (1,))
+ Vu = FeSpace(mesh, P1(), (m,))
+ Vp1 = FeSpace(mesh, DP0(), (1,))
+ Vp2 = FeSpace(mesh, DP1(), (m, d))
+
+ g = FeFunction(Vu)
+ u = FeFunction(Vu)
+ p1 = FeFunction(Vp1)
+ p2 = FeFunction(Vp2)
+ du = FeFunction(Vu)
+ dp1 = FeFunction(Vp1)
+ dp2 = FeFunction(Vp2)
+ nablau = nabla(u)
+ nabladu = nabla(du)
+
+ g.data .= 1
+ u.data .= 0
+ p1.data .= -1.
+ p2.data .= 0
+ du.data .= 0
+ dp1.data .= 0
+ dp2.data .= 0
+
+ alpha1 = 1.
+ alpha2 = 1.
+ beta = 1e-3
+ lambda = 1.
+
+ gamma1 = 1e-3
+ gamma2 = 1e-3
+
+ T(u) = u
+
+ function dp1_update(xloc; 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 :
+ zeros(size(p1))
+ return -p1 + alpha1 / m1 * (T(u + du) - g) - cond
+ end
+
+ interpolate!(p1, dp1_update; g, u, p1, du)
+
+ function dp2_update(xloc; u, nablau, p2, du, nabladu)
+ m2 = max(gamma2, norm(nablau))
+ cond = norm(nablau) >= gamma2 ?
+ dot(nablau, nabladu) / norm(nablau)^2 * p2 :
+ zeros(size(p2))
+ return -p2 + lambda / m2 * (nablau + nabladu) - cond
+ end
+
+ interpolate!(p2, dp2_update; u, nablau, p2, du, nabladu)
+
+ return p1, p2
+end