improve performance

src/function.jl

@@ -115,8 +115,7 @@ myvec(x::Number) = x
 
 function interpolate!(dst::FeFunction, ::P1, expr::Function; params...)
     params = NamedTuple(params)
-    space = dst.space
-    mesh = space.mesh
+    mesh = dst.space.mesh
     for cell in cells(mesh)
 	for f in params
 	    bind!(f, cell)
@@ -130,8 +129,9 @@ function interpolate!(dst::FeFunction, ::P1, expr::Function; params...)
 
 	    opvalues = map(f -> evaluate(f, xloc), params)
 
-	    gdofs = space.dofmap[:, eldof, cell]
-	    dst.data[gdofs] .= myvec(expr(x; opvalues...))
+            res = myvec(expr(x; opvalues...))
+	    gdofs = dst.space.dofmap[:, eldof, cell]
+	    @inbounds dst.data[gdofs] .= res
 	end
     end
 end

src/image.jl

@@ -166,9 +166,13 @@ function project_l2_lagrange!(u::FeFunction, img::AbstractArray)
     a(xloc, u, du, v, dv; f) = dot(u, v)
     l(xloc, v, dv; f) = dot(f, v)
 
-    I = Float64[]
-    J = Float64[]
+    I = Int[]
+    J = Int[]
     V = Float64[]
+    n = length(img) * nrdims * nldofs * nrdims * nldofs
+    sizehint!(I, n)
+    sizehint!(J, n)
+    sizehint!(V, n)
     b = zeros(ndofs(space))
 
     # mesh cells

src/operator.jl

@@ -68,9 +68,9 @@ function assemble(space::FeSpace, a, l; params...)
 
     n = ncells(mesh) * size(qx, 2) *
         nrdims * nldofs * nrdims * nldofs
-    I = zeros(n)
-    J = zeros(n)
-    V = zeros(n)
+    I = Vector{Int}(undef, n)
+    J = Vector{Int}(undef, n)
+    V = Vector{Float64}(undef, n)
     spidx = 0
     b = zeros(ndofs(space))
     gdof = LinearIndices((nrdims, ndofs(space)))