Skip to content
Snippets Groups Projects
Commit 18063745 authored by Stephan Hilb's avatar Stephan Hilb
Browse files

prevent some allocations

parent fe512a16
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment