implement DP0 functions

src/function.jl

-struct P1Mapper
+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
 
-function P1Mapper(mesh)
+struct FeFunction
+    mapper::Mapper
+    data::Vector{Float64}
+end
+
+function FeFunction(mapper)
+    data = Vector{Float64}(undef, mapper.ndofs)
+    return FeFunction(mapper, data)
+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 P1Mapper(mesh, ndofs, map)
+    return Mapper{P1}(mesh, ndofs, map)
 end
 
-struct P1Function
-    mapper::P1Mapper
-    data::Vector{Float64}
-end
+# DP0 Elements (0th order discontinuous Lagrange)
 
-function P1Function(mapper)
-    data = Vector{Float64}(undef, mapper.ndofs)
-    return P1Function(mapper, data)
+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
 
-function interpolate!(dst::P1Function, src::Function)
-    mapper = dst.mapper
+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)
@@ -38,3 +57,30 @@ function interpolate!(dst::P1Function, src::Function)
 	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")
+end
+
+function append_data!(vtkfile, f::FeFunction, ::Mapper{DP0})
+    # FIXME: vector-valued data
+    fdata = reshape(f.data, 1, :)
+    vtk_cell_data(vtkfile, fdata, "f")
+end

src/mesh.jl

+export init_grid, save
+
 using WriteVTK
 
 # 2d, simplex grid
@@ -25,12 +27,6 @@ function init_grid(n)
     return Mesh(vertices, cells)
 end
 
-function append_data!(vtkfile, f::P1Function)
-    # FIXME: vector-valued data
-    fdata = reshape(f.data, 1, :)
-    vtk_point_data(vtkfile, fdata, "f")
-end
-
 function save(filename, mesh::Mesh, fs...)
     cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, mesh.cells[:, i]) for i in axes(mesh.cells, 2)]
     #vertices = [mesh.vertices[i, j] for i in axes(mesh.vertices, 1), j in axes(mesh.vertices, 2)]