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

provide global coordinates for operator

parent 31606b50
No related branches found
No related tags found
No related merge requests found
...@@ -36,12 +36,6 @@ a(op::Poisson, xloc, u, du, v, dv) = dot(du, dv) ...@@ -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] 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, assemble(op::Operator) = assemble(op.space,
(x...; y...) -> a(op, x...; y...), (x...; y...) -> a(op, x...; y...),
(x...; y...) -> l(op, x...; y...); (x...; y...) -> l(op, x...; y...);
...@@ -81,13 +75,14 @@ function assemble(space::FeSpace, a, l; params...) ...@@ -81,13 +75,14 @@ function assemble(space::FeSpace, a, l; params...)
for cell in cells(mesh) for cell in cells(mesh)
foreach(f -> bind!(f, cell), opparams) foreach(f -> bind!(f, cell), opparams)
# cell map is assumed to be constant per cell # 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) delmapinv = inv(delmap)
intel = abs(det(delmap)) intel = abs(det(delmap))
# quadrature points # quadrature points
for k in axes(qx, 2) for k in axes(qx, 2)
xhat = qx[:, k]::SArray xhat = qx[:, k]::SArray
x = elmap(mesh, cell)(xhat)
opvalues = map(f -> evaluate(f, xhat), opparams) opvalues = map(f -> evaluate(f, xhat), opparams)
# local test-function dofs # local test-function dofs
...@@ -97,7 +92,7 @@ function assemble(space::FeSpace, a, l; params...) ...@@ -97,7 +92,7 @@ function assemble(space::FeSpace, a, l; params...)
phij = SArray{Tuple{space.size...}}(qphi[:, jdim, ldofj, k]) phij = SArray{Tuple{space.size...}}(qphi[:, jdim, ldofj, k])
dphij = SArray{Tuple{space.size..., d}}(dqphi[:, :, jdim, ldofj, k] * delmapinv) 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 b[gdofj] += lv
# local trial-function dofs # local trial-function dofs
...@@ -107,7 +102,7 @@ function assemble(space::FeSpace, a, l; params...) ...@@ -107,7 +102,7 @@ function assemble(space::FeSpace, a, l; params...)
phii = SArray{Tuple{space.size...}}(qphi[:, idim, ldofi, k]) phii = SArray{Tuple{space.size...}}(qphi[:, idim, ldofi, k])
dphii = SArray{Tuple{space.size..., d}}(dqphi[:, :, idim, ldofi, k] * delmapinv) 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!(I, gdofi)
push!(J, gdofj) push!(J, gdofj)
push!(V, av) push!(V, av)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment