Support matrix input for L2 projection. (#386)

src/L2_projection.jl

@@ -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 :
@@ -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)).
@@ -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}
 
@@ -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̂.
 
@@ -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

test/test_l2_projection.jl

@@ -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,
@@ -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
@@ -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
@@ -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
@@ -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
@@ -189,4 +205,4 @@ end
     test_projection(2, RefTetrahedron)
     test_projection_mixedgrid()
     test_node_reordering()
-end
+end