partly refactor mesh interface

The mesh interface is a mess. I probably won't have time to clean it up
before implementing refinement though, sadly.

src/function.jl

@@ -31,41 +31,38 @@ end
 function FeSpace(mesh, el::P1, size_=(1,))
     # special case for P1: dofs correspond to vertices
     rdims = prod(size_)
-    ncells = size(mesh.cells, 2)
-    dofmap = Array{Int, 3}(undef, rdims, 3, ncells)
-    for i in CartesianIndices((rdims, 3, ncells))
+    dofmap = Array{Int, 3}(undef, rdims, 3, ncells(mesh))
+    for i in CartesianIndices((rdims, 3, ncells(mesh)))
 	dofmap[i] = rdims * (mesh.cells[i[2], i[3]] - 1) + i[1]
     end
     return FeSpace{typeof(mesh), typeof(el), size_}(
-	mesh, el, dofmap, rdims * size(mesh.vertices, 2))
+	mesh, el, dofmap, rdims * nvertices(mesh))
 end
 
 function FeSpace(mesh, el::DP0, size_=(1,))
     # special case for DP1: dofs correspond to cells
     rdims = prod(size_)
-    ncells = size(mesh.cells, 2)
-    dofmap = Array{Int, 3}(undef, rdims, 1, ncells)
-    for i in CartesianIndices((rdims, 1, ncells))
+    dofmap = Array{Int, 3}(undef, rdims, 1, ncells(mesh))
+    for i in CartesianIndices((rdims, 1, ncells(mesh)))
 	dofmap[i] = rdims * (i[3] - 1) + i[1]
     end
     return FeSpace{typeof(mesh), typeof(el), size_}(
-	mesh, el, dofmap, rdims * ncells)
+	mesh, el, dofmap, rdims * ncells(mesh))
 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 = Array{Int, 3}(undef, rdims, 3, ncells(mesh))
+    for i in LinearIndices((rdims, 3, ncells(mesh)))
 	dofmap[i] = i
     end
     return FeSpace{typeof(mesh), typeof(el), size_}(
-	mesh, el, dofmap, rdims * 3 * ncells)
+	mesh, el, dofmap, rdims * 3 * ncells(mesh))
 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")
+    print("$(nameof(typeof(x))), $(nameof(typeof(x.element))) elements, size $(x.size), $(ndofs(x)) dofs")
 
 function Base.getproperty(obj::FeSpace{<:Any, <:Any, S}, sym::Symbol) where S
     if sym === :size
@@ -75,6 +72,8 @@ function Base.getproperty(obj::FeSpace{<:Any, <:Any, S}, sym::Symbol) where S
     end
 end
 
+ndofs(space::FeSpace) = space.ndofs
+
 # evaluate at local point
 @inline function evaluate(space::FeSpace, ldofs, xloc)
     bv = evaluate_basis(space.element, xloc)
@@ -97,7 +96,7 @@ struct FeFunction{Sp,Ld}
 end
 
 function FeFunction(space; name=string(gensym("f")))
-    data = Vector{Float64}(undef, space.ndofs)
+    data = Vector{Float64}(undef, ndofs(space))
     ldata = zero(MVector{prod(space.size) * ndofs(space.element)})
     return FeFunction(space, data, name, ldata)
 end
@@ -116,13 +115,14 @@ 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 cell in cells(mesh)
 	for f in params
 	    bind!(f, cell)
 	end
-	for eldof in axes(mesh.cells, 1)
-	    xid = mesh.cells[eldof, cell]
-	    x = SArray{Tuple{ndims_domain(mesh)}}(mesh.vertices[:, xid])
+
+	vertices = mesh.vertices[:, mesh.cells[:, cell]]
+	for eldof in 1:nvertices_cell(mesh)
+	    x = SArray{Tuple{ndims_domain(mesh)}}(vertices[:, eldof])
 	    xloc = SA[0. 1. 0.; 0. 0. 1.][:, eldof]
 
 	    opvalues = map(f -> evaluate(f, xloc), params)
@@ -137,10 +137,11 @@ 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 cell in cells(mesh)
 	for f in params
 	    bind!(f, cell)
 	end
+
 	vertices = mesh.vertices[:, mesh.cells[:, cell]]
 	centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(vertices, dims = 2))
 	lcentroid = SA[1/3, 1/3]
@@ -186,7 +187,7 @@ end
 
 function sample(f::FeFunction)
     mesh = f.mapper.mesh
-    for cell in axes(mesh.cells, 2)
+    for cell in cells(mesh)
 	vertices = mesh.vertices[:, mesh.cells[:, cell]]
 	I0 = floor.(Int, vec(minimum(vertices, dims = 2)))
 	I1 = ceil.(Int, vec(maximum(vertices, dims = 2)))

src/mesh.jl

 export init_grid, save
 
 using WriteVTK
+using StaticArrays: SVector
 
 # 2d, simplex grid
 struct Mesh
@@ -14,6 +15,18 @@ Base.show(io::IO, ::MIME"text/plain", f::Mesh) =
 ndims_domain(::Mesh) = 2
 ndims_space(::Mesh) = 2
 nvertices_cell(::Mesh) = 3
+nvertices(x::Mesh) = size(x.vertices, 2)
+ncells(x::Mesh) = size(x.cells, 2)
+
+cells(x::Mesh) = axes(x.cells, 2)
+
+#function Base.getproperty(obj::Mesh, sym::Symbol)
+#    if sym === :vertices
+#	return reinterpret(SVector{ndims_space(obj), Float64}, vec(obj.vertices))
+#    else
+#	return getfield(obj, sym)
+#    end
+#end
 
 function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
     r1 = LinRange(v0[1], v1[1], m + 1)

src/operator.jl

@@ -74,11 +74,11 @@ function assemble(space::FeSpace, a, l; params...)
     I = Float64[]
     J = Float64[]
     V = Float64[]
-    b = zeros(space.ndofs)
-    gdof = LinearIndices((nrdims, space.ndofs))
+    b = zeros(ndofs(space))
+    gdof = LinearIndices((nrdims, ndofs(space)))
 
     # mesh cells
-    for cell in axes(mesh.cells, 2)
+    for cell in cells(mesh)
 	foreach(f -> bind!(f, cell), opparams)
 	# cell map is assumed to be constant per cell
 	delmap = jacobian(x -> elmap(mesh, cell, x), SA[0., 0.])
@@ -116,7 +116,7 @@ function assemble(space::FeSpace, a, l; params...)
 	end
     end
 
-    ngdofs = space.ndofs
+    ngdofs = ndofs(space)
     A = sparse(I, J, V, ngdofs, ngdofs)
     return A, b
 end