finished l2projection

src/function.jl

@@ -80,6 +80,10 @@ Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
 
 interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.space.element, src)
 
+
+myvec(x) = vec(x)
+myvec(x::Number) = x
+
 function interpolate!(dst::FeFunction, ::P1, src::Function)
     space = dst.space
     mesh = space.mesh
@@ -89,7 +93,7 @@ function interpolate!(dst::FeFunction, ::P1, src::Function)
 	    x = mesh.vertices[:, xid]
 
 	    gdofs = space.dofmap[:, eldof, cell]
-	    dst.data[gdofs] .= vec(src(x))
+	    dst.data[gdofs] .= myvec(src(x))
 	end
     end
 end
@@ -102,7 +106,7 @@ function interpolate!(dst::FeFunction, ::DP0, src::Function)
 	centroid = mean(vertices, dims = 2)
 
 	gdofs = space.dofmap[:, 1, cell]
-	dst.data[gdofs] .= vec(src(centroid))
+	dst.data[gdofs] .= myvec(src(centroid))
     end
 end
 

src/operator.jl

@@ -2,7 +2,7 @@ using SparseArrays: sparse
 using LinearAlgebra: det, dot
 using ForwardDiff: jacobian
 
-export Poisson, L2Projection, init_point!, assemble
+export Poisson, L2Projection, init_point!, assemble, assemble_rhs
 
 abstract type Operator end
 
@@ -14,7 +14,7 @@ end
 L2Projection(f) = L2Projection(f.space, f)
 
 params(op::L2Projection) = (f = op.f,)
-a(op::L2Projection, xloc, u, du, v, dv; _...) = dot(u, v)
+a(op::L2Projection, xloc, u, du, v, dv; f) = dot(u, v)
 l(op::L2Projection, xloc, v, dv; f) = dot(f, v)
 
 
@@ -38,6 +38,7 @@ elmap(mesh, cell, x) =
 function assemble(op::Operator)
     space = op.space
     mesh = space.mesh
+    opparams = params(op)
 
     # precompute basis at quadrature points
 
@@ -53,7 +54,7 @@ function assemble(op::Operator)
     for r in 1:nrdims
 	for k in axes(qx, 2)
 	    qphi[r, k, r, :] .= evaluate_basis(space.element, qx[:, k])
-	    dqphi[r, :, k, r, :] .= transpose(jacobian(localbasis, qx[:, k])::Array{Float64,2})
+	    dqphi[r, :, k, r, :] .= transpose(jacobian(x -> evaluate_basis(space.element, x), qx[:, k])::Array{Float64,2})
 	end
     end
 
@@ -62,6 +63,9 @@ function assemble(op::Operator)
     V = Float64[]
     gdof = LinearIndices((nrdims, space.ndofs))
     for cell in axes(mesh.cells, 2)
+	for f in opparams
+	    bind!(f, cell)
+	end
 	delmap = jacobian(x -> elmap(mesh, cell, x), [0., 0.])::Array{Float64,2} # constant on element
 	delmapinv = inv(delmap) # constant on element
 	intel = abs(det(delmap))
@@ -75,12 +79,14 @@ function assemble(op::Operator)
 		# quadrature points
 		for k in axes(qx, 2)
 		    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..., :))
 
-		    gdofv = qw[k] * a(op, xhat, phii, dphii, phij, dphij) * intel
+		    gdofv = qw[k] * a(op, xhat, phii, dphii, phij, dphij; opvalues...) * intel
 
 		    push!(I, gdofi)
 		    push!(J, gdofj)
@@ -102,4 +108,61 @@ function assemble(op::Operator)
     return A
 end
 
+function assemble_rhs(op::Operator)
+    space = op.space
+    mesh = space.mesh
+    opparams = params(op)
+
+    # precompute basis at quadrature points
+
+    qw, qx = quadrature()
+
+    d = size(qx, 1) # domain dimension
+    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, :] .= evaluate_basis(space.element, qx[:, k])
+	    dqphi[r, :, k, r, :] .= transpose(jacobian(x -> evaluate_basis(space.element, x), qx[:, k])::Array{Float64,2})
+	end
+    end
+
+    b = Vector{Float64}(undef, space.ndofs)
+    gdof = LinearIndices((nrdims, space.ndofs))
+    for cell in axes(mesh.cells, 2)
+	for f in opparams
+	    bind!(f, cell)
+	end
+	delmap = jacobian(x -> elmap(mesh, cell, x), [0., 0.])::Array{Float64,2} # constant on element
+	delmapinv = inv(delmap) # constant on element
+	intel = abs(det(delmap))
+
+	# local dofs
+	for jdim in 1:nrdims, ldofj in 1:nldofs
+	    gdofj = space.dofmap[jdim, ldofj, cell]
+
+	    # quadrature points
+	    for k in axes(qx, 2)
+		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..., :))
+
+		gdofv = qw[k] * l(op, xhat, phij, dphij; opvalues...) * intel
+
+		b[gdofj] += gdofv
+	    end
+	end
+    end
+
+    return b
+end
+
+
+