improve performance

src/function.jl

@@ -187,6 +187,7 @@ function project!(dst::FeFunction, ::DP0, expr; params...)
     end
 end
 
+# TODO: interface is lacking, which quadrature should be used?
 function integrate(mesh::Mesh, expr; params...)
     opparams = NamedTuple(params)
 
@@ -196,15 +197,16 @@ function integrate(mesh::Mesh, expr; params...)
     opvalues = map(f -> evaluate(f, SA[0., 0.]), opparams)
     val = zero(expr(SA[0., 0.]; opvalues...))
     for cell in cells(mesh)
+        cellmap = elmap(mesh, cell)
 	foreach(f -> bind!(f, cell), opparams)
 
-	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+	delmap = jacobian(cellmap, SA[0., 0.])
 	intel = abs(det(delmap))
 
 	# quadrature points
 	for k in axes(qx, 2)
 	    xhat = qx[:, k]
-	    x = elmap(mesh, cell)(xhat)
+	    x = cellmap(xhat)
 	    opvalues = map(f -> evaluate(f, xhat), opparams)
 
 	    val += qw[k] * expr(x; opvalues...) * intel

src/operator.jl

@@ -66,9 +66,12 @@ function assemble(space::FeSpace, a, l; params...)
     qphi = SArray{Tuple{nrdims, nrdims, nldofs, nqpts}}(qphi_)
     dqphi = SArray{Tuple{nrdims, d, nrdims, nldofs, nqpts}}(dqphi_)
 
-    I = Float64[]
-    J = Float64[]
-    V = Float64[]
+    n = ncells(mesh) * size(qx, 2) *
+        nrdims * nldofs * nrdims * nldofs
+    I = zeros(n)
+    J = zeros(n)
+    V = zeros(n)
+    spidx = 0
     b = zeros(ndofs(space))
     gdof = LinearIndices((nrdims, ndofs(space)))
 
@@ -104,13 +107,15 @@ function assemble(space::FeSpace, a, l; params...)
 		    dphii = SArray{Tuple{space.size..., d}}(dqphi[:, :, idim, ldofi, k] * delmapinv)
 
 		    av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel
-		    push!(I, gdofi)
-		    push!(J, gdofj)
-		    push!(V, av)
+                    spidx += 1
+                    I[spidx] = gdofi
+                    J[spidx] = gdofj
+                    V[spidx] = av
 		end
 	    end
 	end
     end
+    @assert spidx == n
 
     ngdofs = ndofs(space)
     A = sparse(I, J, V, ngdofs, ngdofs)