diff --git a/src/operator.jl b/src/operator.jl index 76fddf00145f17d3cd21d9c8ce3cc0c13969b088..65a481a09801ccd6b70de1465ab0181620bbf02d 100644 --- a/src/operator.jl +++ b/src/operator.jl @@ -36,12 +36,6 @@ a(op::Poisson, xloc, u, du, v, dv) = dot(du, dv) quadrature() = SA[1/6, 1/6, 1/6], SA[1/6 4/6 1/6; 1/6 1/6 4/6] -function elmap(mesh, cell, x) - A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( - view(mesh.vertices, :, view(mesh.cells, :, cell))) - return A * SA[1 - x[1] - x[2], x[1], x[2]] -end - assemble(op::Operator) = assemble(op.space, (x...; y...) -> a(op, x...; y...), (x...; y...) -> l(op, x...; y...); @@ -81,13 +75,14 @@ function assemble(space::FeSpace, a, l; params...) for cell in cells(mesh) foreach(f -> bind!(f, cell), opparams) # cell map is assumed to be constant per cell - delmap = jacobian(x -> elmap(mesh, cell, x), SA[0., 0.]) + delmap = jacobian(elmap(mesh, cell), SA[0., 0.]) delmapinv = inv(delmap) intel = abs(det(delmap)) # quadrature points for k in axes(qx, 2) xhat = qx[:, k]::SArray + x = elmap(mesh, cell)(xhat) opvalues = map(f -> evaluate(f, xhat), opparams) # local test-function dofs @@ -97,7 +92,7 @@ function assemble(space::FeSpace, a, l; params...) phij = SArray{Tuple{space.size...}}(qphi[:, jdim, ldofj, k]) dphij = SArray{Tuple{space.size..., d}}(dqphi[:, :, jdim, ldofj, k] * delmapinv) - lv = qw[k] * l(xhat, phij, dphij; opvalues...) * intel + lv = qw[k] * l(x, phij, dphij; opvalues...) * intel b[gdofj] += lv # local trial-function dofs @@ -107,7 +102,7 @@ function assemble(space::FeSpace, a, l; params...) phii = SArray{Tuple{space.size...}}(qphi[:, idim, ldofi, k]) dphii = SArray{Tuple{space.size..., d}}(dqphi[:, :, idim, ldofi, k] * delmapinv) - av = qw[k] * a(xhat, phii, dphii, phij, dphij; opvalues...) * intel + av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel push!(I, gdofi) push!(J, gdofj) push!(V, av)