introduce FeSpace in place of Mapper

it makes sense to store shape information in the space

src/function.jl

 using Statistics: mean
-export Mapper, FeFunction, P1, DP0, DP1, interpolate!, sample, cell_contains
+export FeSpace, Mapper, FeFunction, P1, DP0, DP1
+export interpolate!, sample, bind!, evaluate
 
 # Finite Elements
 struct P1 end
 struct DP0 end
 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(::DP1, x) = evaluate_basis(P1(), x)
 
-# Fe: finite element type
-struct Mapper{M,Fe}
+# non-mixed
+# scalar elements
+
+struct FeSpace{M, Fe, S}
     mesh::M
-    ndofs::Int
-    map::Array{Int, 2} # (ldof, cid) -> gdof
+    element::Fe
+    dofmap::Array{Int, 3} # (rdim, eldof, cell) -> gdof
+    ndofs::Int # = maximum(dofmap)
+    size::S
 end
 
-# P1 Elements (1st order Lagrange)
-
-function Mapper{P1}(mesh)
+function FeSpace(mesh, el::P1, size_=(1,))
     # special case for P1: dofs correspond to vertices
-    ndofs = size(mesh.vertices, 2)
-    map = copy(mesh.cells)
-    return Mapper{typeof(mesh),P1}(mesh, ndofs, map)
+    rdims = prod(size_)
+    ncells = size(mesh.cells, 2)
+    dofmap = Array{Int, 3}(undef, rdims, 3, ncells)
+    for i in CartesianIndices((rdims, 3, ncells))
+	dofmap[i] = rdims * (mesh.cells[i[2], i[3]] - 1) + i[1]
     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{typeof(mesh),DP0}(mesh, ndofs, map)
+    return FeSpace(mesh, el, dofmap, rdims * size(mesh.vertices, 2), size_)
 end
 
-# DP1 Elements (1st order discontinuous Lagrange)
-
-function Mapper{DP1}(mesh)
-    # special case for DP0: dofs correspond to cells
-    ndofs = size(mesh.cells, 2)
-    map = collect(reshape(axes(mesh.cells, 2), 1, :))
-    return Mapper{typeof(mesh),DP0}(mesh, ndofs, map)
+function FeSpace(mesh, el::DP0, size_=(1,))
+    # special case for P1: dofs correspond to vertices
+    rdims = prod(size_)
+    ncells = size(mesh.cells, 2)
+    dofmap = Array{Int, 3}(undef, rdims, 1, ncells)
+    for i in CartesianIndices((rdims, 1, ncells))
+	dofmap[i] = rdims * (i[3] - 1) + i[1]
+    end
+    return FeSpace(mesh, el, dofmap, rdims * ncells, size_)
 end
 
-Base.show(io::IO, ::MIME"text/plain", m::Mapper) =
-    print("$(nameof(typeof(m))), $(element(m)) elements, $(size(m.map, 1)) local dofs each")
+Base.show(io::IO, ::MIME"text/plain", x::FeSpace) =
+    print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(x.ndofs) dofs")
 
-element(::Mapper{<:Any,Fe}) where Fe = Fe
+# evaluate at local point
+function evaluate(space::FeSpace, ldofs, xloc)
+    bv = evaluate_basis(space.element, xloc)
+    v = reshape(ldofs, size(space.dofmap)[1:2]) * bv
+    return reshape(v, space.size)
+end
 
 
 
@@ -52,54 +62,59 @@ element(::Mapper{<:Any,Fe}) where Fe = Fe
 #   (ldof, fdims...)
 
 # Array-valued function
-struct FeFunction
-    mapper::Mapper
-    size::NTuple
-    data::Array{Float64, 2} # (rdim, mdof) -> data
+struct FeFunction{Sp}
+    space::Sp
+    data::Vector{Float64} # gdof -> data
     name::String
+    ldata::Vector{Float64} # ldof -> data
 end
 
-function FeFunction(mapper, size=(1,); name=string(gensym("f")))
-    data = Array{Float64, 2}(undef, prod(size), mapper.ndofs)
-    return FeFunction(mapper, size, data, name)
+function FeFunction(space; name=string(gensym("f")))
+    data = Vector{Float64}(undef, space.ndofs)
+    ldata = Vector{Float64}(undef, prod(size(space.dofmap)[1:2]))
+    return FeFunction(space, data, name, ldata)
 end
 
 Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
-    print("$(nameof(typeof(f))), size $(f.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.mapper, src)
+interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.space.element, src)
 
-function interpolate!(dst::FeFunction, mapper::Mapper{<:Any, 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]
+function interpolate!(dst::FeFunction, ::P1, src::Function)
+    space = dst.space
+    mesh = space.mesh
+    for cell in axes(mesh.cells, 2)
+	for eldof in axes(mesh.cells, 1)
+	    xid = mesh.cells[eldof, cell]
 	    x = mesh.vertices[:, xid]
 
-	    gdof = mapper.map[ldof, cid]
-	    dst.data[:, gdof] .= vec(src(x))
+	    gdofs = space.dofmap[:, eldof, cell]
+	    dst.data[gdofs] .= vec(src(x))
 	end
     end
 end
 
-function interpolate!(dst::FeFunction, mapper::Mapper{<:Any, DP0}, src::Function)
-    mesh = mapper.mesh
-    for cid in axes(mesh.cells, 2)
-	vertices = mesh.vertices[:, mesh.cells[:, cid]]
+function interpolate!(dst::FeFunction, ::DP0, src::Function)
+    space = dst.space
+    mesh = space.mesh
+    for cell in axes(mesh.cells, 2)
+	vertices = mesh.vertices[:, mesh.cells[:, cell]]
 	centroid = mean(vertices, dims = 2)
 
-	gdof = mapper.map[1, cid]
-	dst.data[:, gdof] .= vec(src(centroid))
+	gdofs = space.dofmap[:, 1, cell]
+	dst.data[gdofs] .= vec(src(centroid))
     end
 end
 
-function cell_contains(cell, v)
-    J = jacobian(x -> cell * [1 - x[1] - x[2], x[1], x[2]], [0., 0.])
-    λ = J \ (v - cell[:, 1])
-    return all(λ .>= 0) && sum(λ) <= 1
+
+function bind!(f::FeFunction, cell)
+    f.ldata .= vec(f.data[f.space.dofmap[:, :, cell]])
+    return f
 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}
@@ -117,7 +132,6 @@ function sample(f::FeFunction)
 
 	for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2])
 	    if all(xloc .>= 0) && sum(xloc) <= 1
-		eval_point(f, 
 		# eval point
 	    end
 	end
@@ -125,13 +139,13 @@ function sample(f::FeFunction)
 end
 
 
-append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.mapper)
+append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element)
 
-function append_data!(vtkfile, f::FeFunction, ::Mapper{<:Any, P1})
+function append_data!(vtkfile, f::FeFunction, ::P1)
     vtk_point_data(vtkfile, f.data, f.name)
 end
 
-function append_data!(vtkfile, f::FeFunction, ::Mapper{<:Any, DP0})
+function append_data!(vtkfile, f::FeFunction, ::DP0)
     vtk_cell_data(vtkfile, f.data, f.name)
 end
 

src/mesh.jl

@@ -36,6 +36,13 @@ init_grid(img::Array{<:Any, 2}, type=:vertex) =
 	init_grid(size(img, 1) - 1, size(img, 2) - 1, (1.0, 1.0), size(img)) :
 	init_grid(size(img, 1), size(img, 2), (0.5, 0.5), size(img) .- (0.5, 0.5))
 
+#function cell_contains(mesh, cell, v)
+#    geo = mesh.vertices[:, mesh.cells[:, cell]]
+#    J = jacobian(x -> geo * [1 - x[1] - x[2], x[1], x[2]], [0., 0.])
+#    λ = J \ (v - geo[:, 1])
+#    return all(λ .>= 0) && sum(λ) <= 1
+#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)]

src/operator.jl

@@ -6,22 +6,21 @@ export Poisson, L2Projection, init_point!, assemble
 
 abstract type Operator end
 
-struct L2Projection <: Operator
-    f::FeFunction # target
-    # space info
-    mapper::Mapper
-    size::NTuple
+struct L2Projection{Sp, F} <: Operator
+    space::Sp
+    f::F # target
 end
 
-L2Projection(f, gspace) = L2Projection(f, g.mapper, g.size)
+L2Projection(f) = L2Projection(f.space, f)
 
-a(op::L2Projection, xloc, u, du, v, dv) = dot(u, v)
-l(op::L2Projection, xloc, v, dv) = dot(op.f, v)
+params(op::L2Projection) = (f = op.f,)
+a(op::L2Projection, xloc, u, du, v, dv; _...) = dot(u, v)
+l(op::L2Projection, xloc, v, dv; f) = dot(f, v)
 
 
-struct Poisson{M,S} <: Operator
-    mapper::M
-    size::S
+
+struct Poisson{Sp} <: Operator
+    space::Sp
 end
 
 a(op::Poisson, xloc, u, du, v, dv) = dot(du, dv)
@@ -31,32 +30,29 @@ a(op::Poisson, xloc, u, du, v, dv) = dot(du, dv)
 # degree 2 quadrature on triangle
 quadrature() = [1/6, 1/6, 1/6], [1/6 4/6 1/6; 1/6 1/6 4/6]
 
-# evaluate all (1-dim) local basis functions against x
-# linear 1st order lagrange
-localbasis(x) = [1 - x[1] - x[2], x[1], x[2]]
 
 elmap(mesh, cell, x) =
     mesh.vertices[:, mesh.cells[:, cell]] * [1 - x[1] - x[2], x[1], x[2]]
 
 
 function assemble(op::Operator)
-    mapper = op.mapper
-    mesh = mapper.mesh
+    space = op.space
+    mesh = space.mesh
 
     # precompute basis at quadrature points
 
     qw, qx = quadrature()
 
     d = size(qx, 1) # domain dimension
-    nrdims = prod(op.size)
-    nldofs = size(mapper.map, 1) # number of local dofs (not counting range dimensions)
+    nrdims = prod(space.size)
+    nldofs = size(space.dofmap, 2) # number of local dofs (not counting range dimensions)
     nqpts = length(qw) # number of quadrature points
 
     qphi = zeros(nrdims, nqpts, nrdims, nldofs)
     dqphi = zeros(nrdims, d, nqpts, nrdims, nldofs)
     for r in 1:nrdims
 	for k in axes(qx, 2)
-	    qphi[r, k, r, :] .= localbasis(qx[:, k])
+	    qphi[r, k, r, :] .= evaluate_basis(space.element, qx[:, k])
 	    dqphi[r, :, k, r, :] .= transpose(jacobian(localbasis, qx[:, k])::Array{Float64,2})
 	end
     end
@@ -64,7 +60,7 @@ function assemble(op::Operator)
     I = Float64[]
     J = Float64[]
     V = Float64[]
-    gdof = LinearIndices((nrdims, mapper.ndofs))
+    gdof = LinearIndices((nrdims, space.ndofs))
     for cell in axes(mesh.cells, 2)
 	delmap = jacobian(x -> elmap(mesh, cell, x), [0., 0.])::Array{Float64,2} # constant on element
 	delmapinv = inv(delmap) # constant on element
@@ -72,19 +68,17 @@ function assemble(op::Operator)
 
 	# local dofs
 	for idim in 1:nrdims, ldofi in 1:nldofs
-	    mdofi = mapper.map[ldofi, cell]
-	    gdofi = gdof[idim, mdofi]
+	    gdofi = space.dofmap[idim, ldofi, cell]
 	    for jdim in 1:nrdims, ldofj in 1:nldofs
-		mdofj = mapper.map[ldofj, cell]
-		gdofj = gdof[jdim, mdofj]
+		gdofj = space.dofmap[jdim, ldofj, cell]
 
 		# quadrature points
 		for k in axes(qx, 2)
 		    xhat = qx[:, k]
-		    phii = reshape(qphi[:, k, idim, ldofi], op.size)
-		    dphii = reshape(dqphi[:, :, k, idim, ldofi] * delmapinv, (op.size..., :))
-		    phij = reshape(qphi[:, k, jdim, ldofj], op.size)
-		    dphij = reshape(dqphi[:, :, k, jdim, ldofj] * delmapinv, (op.size..., :))
+		    phii = reshape(qphi[:, k, idim, ldofi], space.size)
+		    dphii = reshape(dqphi[:, :, k, idim, ldofi] * delmapinv, (space.size..., :))
+		    phij = reshape(qphi[:, k, jdim, ldofj], space.size)
+		    dphij = reshape(dqphi[:, :, k, jdim, ldofj] * delmapinv, (space.size..., :))
 
 		    gdofv = qw[k] * a(op, xhat, phii, dphii, phij, dphij) * intel
 
@@ -103,7 +97,7 @@ function assemble(op::Operator)
 	end
     end
 
-    ngdofs = nrdims * mapper.ndofs
+    ngdofs = space.ndofs
     A = sparse(I, J, V, ngdofs, ngdofs)
     return A
 end