Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support matrix input for L2 projection. #386

Merged
merged 1 commit into from
Oct 12, 2021
+53 −21
Merged

Support matrix input for L2 projection. #386

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 30 additions & 14 deletions src/L2_projection.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ function L2Projector(
_, vertex_dict, _, _ = __close!(dh)

M = _assemble_L2_matrix(fe_values_mass, set, dh) # the "mass" matrix
M_cholesky = cholesky(M) # TODO maybe have a lazy eval instead of precomputing? / JB
M_cholesky = cholesky(M)

# For deprecated API
fe_values = qr_rhs === nothing ? nothing :
Expand Down Expand Up @@ -153,7 +153,7 @@ end


"""
project(proj::L2Projector, vals::Vector{Vector{T}}}, qr_rhs::QuadratureRule; project_to_nodes=true)
project(proj::L2Projector, vals, qr_rhs::QuadratureRule; project_to_nodes=true)
Makes a L2 projection of data `vals` to the nodes of the grid using the projector `proj`
(see [`L2Projector`](@ref)).
Expand All @@ -167,22 +167,32 @@ where ``f`` is the data to project, i.e. `vals`.
The data `vals` should be a vector, with length corresponding to number of elements, of vectors,
with length corresponding to number of quadrature points per element, matching the number of points in `qr_rhs`.
Example scalar input data:
Alternatively, `vals` can be a matrix, with number of columns corresponding to number of elements,
and number of rows corresponding to number of points in `qr_rhs`.
Example (scalar) input data:
```julia
vals = [
[0.44, 0.98, 0.32], # data for quadrature point 1, 2, 3 of element 1
[0.29, 0.48, 0.55], # data for quadrature point 1, 2, 3 of element 2
# ...
]
```
or equivalent in matrix form:
```julia
vals = [
0.44 0.29 # ...
0.98 0.48 # ...
0.32 0.55 # ...
]
```
Supported data types to project are `Number`s and `AbstractTensor`s.
If the parameter `project_to_nodes` is `true`, then the projection returns the values in the order of the mesh nodes
(suitable format for exporting). If `false`, it returns the values corresponding to the degrees of freedom for a scalar
field over the domain, which is useful if one wants to interpolate the projected values.
"""
function project(proj::L2Projector,
vars::Vector{Vector{T}},
vars::AbstractVector{<:AbstractVector{T}},
qr_rhs::Union{QuadratureRule,Nothing}=nothing;
project_to_nodes::Bool=true) where T <: Union{Number, AbstractTensor}

Expand All @@ -192,26 +202,29 @@ function project(proj::L2Projector,
CellScalarValues(qr_rhs, proj.func_ip, proj.geom_ip)

M = T <: AbstractTensor ? length(vars[1][1].data) : 1
make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1]

projected_vals = _project(vars, proj, fe_values, M)
projected_vals = _project(vars, proj, fe_values, M, T)::Vector{T}
if project_to_nodes
# NOTE we may have more projected values than verticies in the mesh => not all values are returned
nnodes = getnnodes(proj.dh.grid)
reordered_vals = fill(NaN, nnodes, size(projected_vals, 2))
reordered_vals = fill(convert(T, NaN * zero(T)), nnodes)
for node = 1:nnodes
if haskey(proj.node2dof_map, node)
reordered_vals[node, :] = projected_vals[proj.node2dof_map[node], :]
if (k = get(proj.node2dof_map, node, nothing); k !== nothing)
@assert length(k) == 1
reordered_vals[node] = projected_vals[k[1]]
end
end
return T[make_T(reordered_vals[i,:]) for i=1:size(reordered_vals, 1)]
return reordered_vals
else
# convert back to the original tensor type
return T[make_T(projected_vals[i,:]) for i=1:size(projected_vals, 1)]
return projected_vals
end
end
function project(p::L2Projector, vars::AbstractMatrix, qr_rhs::QuadratureRule; project_to_nodes=true)
# TODO: Random access into vars is required for now, hence the collect
return project(p, collect(eachcol(vars)), qr_rhs; project_to_nodes=project_to_nodes)
end

function _project(vars, proj::L2Projector, fe_values::Values, M::Integer)
function _project(vars, proj::L2Projector, fe_values::Values, M::Integer, ::Type{T}) where {T}
# Assemble the multi-column rhs, f = ∭( v ⋅ x̂ )dΩ
# The number of columns corresponds to the length of the data-tuple in the tensor x̂.

Expand Down Expand Up @@ -251,6 +264,9 @@ function _project(vars, proj::L2Projector, fe_values::Values, M::Integer)
end

# solve for the projected nodal values
return proj.M_cholesky \ f
projected_vals = proj.M_cholesky \ f

# Recast to original input type
make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1]
return T[make_T(x) for x in eachrow(projected_vals)]
end
30 changes: 23 additions & 7 deletions test/test_l2_projection.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,14 +39,16 @@ function test_projection(order, refshape)
# Now recover the nodal values using a L2 projection.
proj = L2Projector(ip, grid; geom_ip=ip_geom)
point_vars = project(proj, qp_values, qr)
qp_values_matrix = reduce(hcat, qp_values)
point_vars_2 = project(proj, qp_values_matrix, qr)
## Old API with fe values as first arg
proj2 = @test_deprecated L2Projector(cv, ip, grid)
point_vars_2 = @test_deprecated project(qp_values, proj2)
point_vars_3 = @test_deprecated project(qp_values, proj2)
## Old API with qr as first arg
proj3 = @test_deprecated L2Projector(qr, ip, grid)
point_vars_3 = @test_deprecated project(qp_values, proj3)
point_vars_4 = @test_deprecated project(qp_values, proj3)

@test point_vars ≈ point_vars_2 ≈ point_vars_3
@test point_vars ≈ point_vars_2 ≈ point_vars_3 ≈ point_vars_4

if order == 1
# A linear approximation can not recover a quadratic solution,
Expand All @@ -72,7 +74,12 @@ function test_projection(order, refshape)
# Tensor
f_tensor(x) = Tensor{2,2,Float64}((f(x),2*f(x),3*f(x),4*f(x)))
qp_values = analytical(f_tensor)
qp_values_matrix = reduce(hcat, qp_values)::Matrix
point_vars = project(proj, qp_values, qr)
point_vars_2 = project(proj, qp_values_matrix, qr)

@test point_vars ≈ point_vars_2

if order == 1
ae = [Tensor{2,2,Float64}((f_approx(i),2*f_approx(i),3*f_approx(i),4*f_approx(i))) for i in 1:4]
elseif order == 2
Expand All @@ -83,7 +90,12 @@ function test_projection(order, refshape)
# SymmetricTensor
f_stensor(x) = SymmetricTensor{2,2,Float64}((f(x),2*f(x),3*f(x)))
qp_values = analytical(f_stensor)
qp_values_matrix = reduce(hcat, qp_values)
point_vars = project(proj, qp_values, qr)
point_vars_2 = project(proj, qp_values_matrix, qr)

@test point_vars ≈ point_vars_2

if order == 1
ae = [SymmetricTensor{2,2,Float64}((f_approx(i),2*f_approx(i),3*f_approx(i))) for i in 1:4]
elseif order == 2
Expand Down Expand Up @@ -132,23 +144,26 @@ function test_projection_mixedgrid()
ae = compute_vertex_values(mesh, f)
# analytical values
qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]]
qp_values_matrix = reduce(hcat, qp_values)

# Now recover the nodal values using a L2 projection.
# Assume f would only exist on the first cell, we project it to the nodes of the
# 1st cell while ignoring the rest of the domain. NaNs should be stored in all
# nodes that do not belong to the 1st cell
proj = L2Projector(ip, mesh; geom_ip=ip_geom, set=1:1)
point_vars = project(proj, qp_values, qr)
point_vars_2 = project(proj, qp_values_matrix, qr)
## Old API with fe values as first arg
proj = @test_deprecated L2Projector(cv, ip, mesh, 1:1)
point_vars_2 = @test_deprecated project(qp_values, proj)
point_vars_3 = @test_deprecated project(qp_values, proj)
## Old API with qr as first arg
proj = @test_deprecated L2Projector(qr, ip, mesh, 1:1)
point_vars_3 = @test_deprecated project(qp_values, proj)
point_vars_4 = @test_deprecated project(qp_values, proj)

# In the nodes of the 1st cell we should recover the field
for node in mesh.cells[1].nodes
@test ae[node] ≈ point_vars[node] ≈ point_vars_2[node] ≈ point_vars_3[node]
@test ae[node] ≈ point_vars[node] ≈ point_vars_2[node] ≈ point_vars_3[node] ≈
point_vars_4[node]
end

# in all other nodes we should have NaNs
Expand All @@ -157,6 +172,7 @@ function test_projection_mixedgrid()
@test isnan(point_vars[node][d1, d2])
@test isnan(point_vars_2[node][d1, d2])
@test isnan(point_vars_3[node][d1, d2])
@test isnan(point_vars_4[node][d1, d2])
end
end
end
Expand Down Expand Up @@ -189,4 +205,4 @@ end
test_projection(2, RefTetrahedron)
test_projection_mixedgrid()
test_node_reordering()
end
end