Select Git revision
function.jl
function.jl 2.07 KiB
using Statistics: mean
export Mapper, FeFunction, P1, DP0, DP1, interpolate!
struct P1 end
struct DP0 end
struct DP1 end
struct Mapper{T}
mesh
ndofs::Int
map::Array{Int, 2}
end
struct FeFunction
mapper::Mapper
data::Vector{Float64}
name::String
end
function FeFunction(mapper, name=string(gensym("f")))
data = Vector{Float64}(undef, mapper.ndofs)
return FeFunction(mapper, data, name)
end
# P1 Elements (1st order Lagrange)
function Mapper{P1}(mesh)
# special case for P1: dofs correspond to vertices
ndofs = size(mesh.vertices, 2)
map = copy(mesh.cells)
return Mapper{P1}(mesh, ndofs, map)
end
# DP0 Elements (0th order discontinuous Lagrange)
function Mapper{DP0}(mesh)
# special case for DP0: dofs correspond to cells
ndofs = size(mesh.cells, 2)
map = collect(reshape(axes(mesh.cells, 2), 1, :))
return Mapper{DP0}(mesh, ndofs, map)
end
interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.mapper, src)
function interpolate!(dst::FeFunction, mapper::Mapper{P1}, src::Function)
mesh = mapper.mesh
for cid in axes(mesh.cells, 2)
for ldof in axes(mapper.map, 1)
# TODO: how do ldofs generically map to position?
xid = mesh.cells[ldof, cid]
x = mesh.vertices[:, xid]
# currently scalar
fx = src(x)
gdof = mapper.map[ldof, cid]
dst.data[gdof] = fx
end
end
end
function interpolate!(dst::FeFunction, mapper::Mapper{DP0}, src::Function)
mesh = mapper.mesh
for cid in axes(mesh.cells, 2)
vertices = mesh.vertices[:, mesh.cells[:, cid]]
centroid = mean(vertices, dims=2)
fx = src(centroid)
gdof = mapper.map[1, cid]
dst.data[gdof] = fx
end
end
append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.mapper)
function append_data!(vtkfile, f::FeFunction, ::Mapper{P1})
# FIXME: vector-valued data
fdata = reshape(f.data, 1, :)
vtk_point_data(vtkfile, fdata, f.name)
end
function append_data!(vtkfile, f::FeFunction, ::Mapper{DP0})
# FIXME: vector-valued data
fdata = reshape(f.data, 1, :)
vtk_cell_data(vtkfile, fdata, f.name)
end