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

flesh out DP1 and add derivative function

parent 18063745
No related branches found
No related tags found
No related merge requests found
...@@ -5,6 +5,8 @@ include("mesh.jl") ...@@ -5,6 +5,8 @@ include("mesh.jl")
include("function.jl") include("function.jl")
include("operator.jl") include("operator.jl")
include("run.jl")
function clip_line(r0, r1, l0, l1) function clip_line(r0, r1, l0, l1)
d_l = -(l1[1] - l0[1]) d_l = -(l1[1] - l0[1])
d_r = +(l1[1] - l0[1]) d_r = +(l1[1] - l0[1])
......
using Statistics: mean using Statistics: mean
export FeSpace, Mapper, FeFunction, P1, DP0, DP1 export FeSpace, Mapper, FeFunction, P1, DP0, DP1
export interpolate!, sample, bind!, evaluate export interpolate!, sample, bind!, evaluate, nabla
# Finite Elements # Finite Elements
struct P1 end struct P1 end
...@@ -10,7 +10,7 @@ struct DP1 end ...@@ -10,7 +10,7 @@ struct DP1 end
# FIXME: should be static vectors # FIXME: should be static vectors
# evaluate all (1-dim) local basis functions against x # evaluate all (1-dim) local basis functions against x
evaluate_basis(::P1, x) = [1 - x[1] - x[2], x[1], x[2]] 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) evaluate_basis(::DP1, x) = evaluate_basis(P1(), x)
# non-mixed # non-mixed
...@@ -36,7 +36,7 @@ function FeSpace(mesh, el::P1, size_=(1,)) ...@@ -36,7 +36,7 @@ function FeSpace(mesh, el::P1, size_=(1,))
end end
function FeSpace(mesh, el::DP0, size_=(1,)) 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_) rdims = prod(size_)
ncells = size(mesh.cells, 2) ncells = size(mesh.cells, 2)
dofmap = Array{Int, 3}(undef, rdims, 1, ncells) dofmap = Array{Int, 3}(undef, rdims, 1, ncells)
...@@ -46,6 +46,17 @@ function FeSpace(mesh, el::DP0, size_=(1,)) ...@@ -46,6 +46,17 @@ function FeSpace(mesh, el::DP0, size_=(1,))
return FeSpace(mesh, el, dofmap, rdims * ncells, size_) return FeSpace(mesh, el, dofmap, rdims * ncells, size_)
end 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) = Base.show(io::IO, ::MIME"text/plain", x::FeSpace) =
print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(x.ndofs) dofs") print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(x.ndofs) dofs")
...@@ -78,38 +89,54 @@ end ...@@ -78,38 +89,54 @@ end
Base.show(io::IO, ::MIME"text/plain", f::FeFunction) = Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
print("$(nameof(typeof(f))), size $(f.space.size) with $(length(f.data)) dofs") 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) = vec(x)
myvec(x::Number) = 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 space = dst.space
mesh = space.mesh mesh = space.mesh
for cell in axes(mesh.cells, 2) for cell in axes(mesh.cells, 2)
for f in params
bind!(f, cell)
end
for eldof in axes(mesh.cells, 1) for eldof in axes(mesh.cells, 1)
xid = mesh.cells[eldof, cell] xid = mesh.cells[eldof, cell]
x = mesh.vertices[:, xid] x = mesh.vertices[:, xid]
xloc = [0. 1. 0.; 0. 0. 1.][:, eldof]
opvalues = map(f -> evaluate(f, xloc), params)
gdofs = space.dofmap[:, eldof, cell] gdofs = space.dofmap[:, eldof, cell]
dst.data[gdofs] .= myvec(src(x)) dst.data[gdofs] .= myvec(expr(x; opvalues...))
end end
end end
end end
function interpolate!(dst::FeFunction, ::DP0, src::Function) function interpolate!(dst::FeFunction, ::DP0, expr::Function; params...)
params = NamedTuple(params)
space = dst.space space = dst.space
mesh = space.mesh mesh = space.mesh
for cell in axes(mesh.cells, 2) for cell in axes(mesh.cells, 2)
for f in params
bind!(f, cell)
end
vertices = mesh.vertices[:, mesh.cells[:, cell]] vertices = mesh.vertices[:, mesh.cells[:, cell]]
centroid = mean(vertices, dims = 2) centroid = mean(vertices, dims = 2)
lcentroid = [1/3, 1/3]
opvalues = map(f -> evaluate(f, lcentroid), params)
gdofs = space.dofmap[:, 1, cell] gdofs = space.dofmap[:, 1, cell]
dst.data[gdofs] .= myvec(src(centroid)) dst.data[gdofs] .= myvec(expr(centroid; opvalues...))
end end
end end
interpolate!(dst::FeFunction, ::DP1, expr::Function; params...) =
interpolate!(dst, P1(), expr; params...)
function bind!(f::FeFunction, cell) function bind!(f::FeFunction, cell)
f.ldata .= vec(f.data[f.space.dofmap[:, :, cell]]) f.ldata .= vec(f.data[f.space.dofmap[:, :, cell]])
...@@ -119,11 +146,22 @@ end ...@@ -119,11 +146,22 @@ end
# evaluate at local point (needs bind! call before) # evaluate at local point (needs bind! call before)
evaluate(f::FeFunction, x) = evaluate(f.space, f.ldata, x) evaluate(f::FeFunction, x) = evaluate(f.space, f.ldata, x)
struct CellFunction struct Derivative{F}
f::FeFunction f::F
dofs::Vector{Float64} 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 end
function sample(f::FeFunction) function sample(f::FeFunction)
mesh = f.mapper.mesh mesh = f.mapper.mesh
for cell in axes(mesh.cells, 2) for cell in axes(mesh.cells, 2)
...@@ -143,6 +181,7 @@ function sample(f::FeFunction) ...@@ -143,6 +181,7 @@ function sample(f::FeFunction)
end end
append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element) append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element)
function append_data!(vtkfile, f::FeFunction, ::P1) function append_data!(vtkfile, f::FeFunction, ::P1)
......
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment