diff --git a/src/operator.jl b/src/operator.jl
index 39f5e47673c29c2af39c107ca44e4b869facecb7..16305b16dbf2d22934cae926a93338251bd3e30f 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -6,6 +6,10 @@ export Poisson, L2Projection, init_point!, assemble, assemble_rhs
 
 abstract type Operator end
 
+Base.show(io::IO, ::MIME"text/plain", x::Operator) =
+    print("$(nameof(typeof(x)))")
+
+
 struct L2Projection{Sp, F} <: Operator
     space::Sp
     f::F # target
@@ -58,6 +62,11 @@ function assemble(op::Operator)
 	end
     end
 
+    xhat = zeros(d)
+    phii = zeros(space.size...)
+    dphii = zeros(space.size..., d)
+    phij = zeros(space.size...)
+    dphij = zeros(space.size..., d)
     I = Float64[]
     J = Float64[]
     V = Float64[]
@@ -78,13 +87,13 @@ function assemble(op::Operator)
 
 		# quadrature points
 		for k in axes(qx, 2)
-		    xhat = qx[:, k]
+		    xhat .= qx[:, k]
 		    opvalues = map(f -> evaluate(f, xhat), opparams)
 
-		    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..., :))
+		    phii .= reshape(view(qphi, :, k, idim, ldofi), space.size)
+		    dphii .= reshape(view(dqphi, :, :, k, idim, ldofi) * delmapinv, (space.size..., :))
+		    phij .= reshape(view(qphi, :, k, jdim, ldofj), space.size)
+		    dphij .= reshape(view(dqphi, :, :, k, jdim, ldofj) * delmapinv, (space.size..., :))
 
 		    gdofv = qw[k] * a(op, xhat, phii, dphii, phij, dphij; opvalues...) * intel
 
@@ -92,12 +101,6 @@ function assemble(op::Operator)
 		    push!(J, gdofj)
 		    push!(V, gdofv)
 
-		    #display(phii)
-		    #display(dphii)
-		    #display(phij)
-		    #display(dphij)
-
-		    #println("$cell, $gdofi, $gdofj, $gdofv")
 		end
 	    end
 	end
@@ -131,7 +134,10 @@ function assemble_rhs(op::Operator)
 	end
     end
 
-    b = Vector{Float64}(undef, space.ndofs)
+    xhat = zeros(d)
+    phij = zeros(space.size...)
+    dphij = zeros(space.size..., d)
+    b = zeros(space.ndofs)
     gdof = LinearIndices((nrdims, space.ndofs))
     for cell in axes(mesh.cells, 2)
 	for f in opparams
@@ -147,11 +153,11 @@ function assemble_rhs(op::Operator)
 
 	    # quadrature points
 	    for k in axes(qx, 2)
-		xhat = qx[:, k]
+		xhat .= qx[:, k]
 		opvalues = map(f -> evaluate(f, xhat), opparams)
 
-		phij = reshape(qphi[:, k, jdim, ldofj], space.size)
-		dphij = reshape(dqphi[:, :, k, jdim, ldofj] * delmapinv, (space.size..., :))
+		phij .= reshape(qphi[:, k, jdim, ldofj], space.size)
+		dphij .= reshape(dqphi[:, :, k, jdim, ldofj] * delmapinv, (space.size..., :))
 
 		gdofv = qw[k] * l(op, xhat, phij, dphij; opvalues...) * intel