diff --git a/src/feoperator/projection_newapi.jl b/src/feoperator/projection_newapi.jl index d6c667bb0..1d5617523 100644 --- a/src/feoperator/projection_newapi.jl +++ b/src/feoperator/projection_newapi.jl @@ -120,6 +120,8 @@ end _may_scalar(a::SVector{1}) = a[1] _may_scalar(a) = a +cell_mean(u::MeshData{CellData}, ::Measure) = u + function build_mass_matrix(u::AbstractFEFunction, dΩ::AbstractMeasure) U = get_fespace(u) V = TestFESpace(U) diff --git a/src/writers/vtk.jl b/src/writers/vtk.jl index a8be43567..638e69a88 100644 --- a/src/writers/vtk.jl +++ b/src/writers/vtk.jl @@ -321,6 +321,19 @@ function _vtk_coords_from_lagrange(shape::Union{Square, Cube}, degree) end """ + write_vtk_lagrange( + basename::String, + vars::Dict{String, F}, + mesh::AbstractMesh, + it::Integer = -1, + time::Real = 0.0; + mesh_degree::Integer = 1, + discontinuous::Bool = true, + functionSpaceType::AbstractFunctionSpaceType = Lagrange(), + collection_append::Bool = false, + vtk_kwargs..., + ) where {F <: AbstractLazy} + write_vtk_lagrange( basename::String, vars::Dict{String, F}, @@ -328,15 +341,25 @@ end U_export::AbstractFESpace, it::Integer = -1, time::Real = 0.0; - collection_append = false, + collection_append::Bool = false, vtk_kwargs..., ) where {F <: AbstractLazy} -Write the provided FEFunction/MeshCellData/CellFunction on the mesh with -the precision of the Lagrange FESpace provided. +Write the provided FEFunction/MeshCellData/CellFunction on a mesh of order `mesh_degree` with the nodes +defined by the `functionSpaceType` (only `Lagrange{:Uniform}` supported for now). The boolean `discontinuous` +indicate if the node values should be discontinuous or not. `vars` is a dictionnary of variable name => Union{FEFunction,MeshCellData,CellFunction} to write. +Cell-centered data should be wrapped as `MeshCellData` in the `vars` dict. To convert an `AbstractLazy` to +a `MeshCellData`, simply use `Bcube.cell_mean` or `MeshCellData ∘ Bcube.var_on_centers`. + +# Remarks: +* If `mesh` is of degree `d`, the solution will be written on a mesh of degree `mesh_degree`, even if this +number is different from `d`. +* The degree of the input FEFunction (P1, P2, P3, ...) is not used to define the nodes where the solution is +written, only `mesh` and `mesh_degree` matter. The FEFunction is simply evaluated on the aforementionned nodes. + # Example ```julia mesh = rectangle_mesh(6, 7; xmin = -1, xmax = 1.0, ymin = -1, ymax = 1.0) @@ -346,13 +369,13 @@ projection_l2!(u, f_u, mesh) vars = Dict("f_u" => f_u, "u" => u, "grad_u" => ∇(u)) -for degree_export in 1:5 - U_export = TrialFESpace(FunctionSpace(:Lagrange, degree_export), mesh) +for mesh_degree in 1:5 Bcube.write_vtk_lagrange( joinpath(@__DIR__, "output"), vars, - mesh, - U_export, + mesh; + mesh_degree, + discontinuous = false ) end ``` @@ -362,6 +385,35 @@ end - `collection_append` is not named `append` to enable passing correct `kwargs` to `vtk_grid` - remove (once fully validated) : `write_vtk_discontinuous` """ +function write_vtk_lagrange( + basename::String, + vars::Dict{String, F}, + mesh::AbstractMesh, + it::Integer = -1, + time::Real = 0.0; + mesh_degree::Integer = 1, + discontinuous::Bool = true, + functionSpaceType::AbstractFunctionSpaceType = Lagrange(), + collection_append::Bool = false, + vtk_kwargs..., +) where {F <: AbstractLazy} + U_export = TrialFESpace( + FunctionSpace(functionSpaceType, mesh_degree), + mesh; + isContinuous = !discontinuous, + ) + write_vtk_lagrange( + basename, + vars, + mesh, + U_export, + it, + time; + collection_append, + vtk_kwargs..., + ) +end + function write_vtk_lagrange( basename::String, vars::Dict{String, F}, @@ -372,12 +424,6 @@ function write_vtk_lagrange( collection_append = false, vtk_kwargs..., ) where {F <: AbstractLazy} - # extract all `MeshCellData` in `vars` as this type of variable - # will not be interpolated to nodes and will be written with - # the `VTKCellData` attribute. - vars_cell = filter(((k, v),) -> v isa MeshData{<:CellData}, vars) - vars_point = filter(((k, v),) -> k ∉ keys(vars_cell), vars) - # FE space stuff fs_export = get_function_space(U_export) @assert get_type(fs_export) <: Lagrange "Only FunctionSpace of type Lagrange are supported for now" @@ -385,6 +431,12 @@ function write_vtk_lagrange( dhl_export = Bcube._get_dhl(U_export) nd = get_ndofs(dhl_export) + # extract all `MeshCellData` in `vars` as this type of variable + # will not be interpolated to nodes and will be written with + # the `VTKCellData` attribute. + vars_cell = filter(((k, v),) -> v isa MeshData{<:CellData}, vars) + vars_point = filter(((k, v),) -> k ∉ keys(vars_cell), vars) + # Get ncomps and type of each `point` variable type_dim = map(var -> get_return_type_and_codim(var, mesh), values(vars_point)) diff --git a/test/writers/test_vtk.jl b/test/writers/test_vtk.jl index bdf56f137..cb49b13d9 100644 --- a/test/writers/test_vtk.jl +++ b/test/writers/test_vtk.jl @@ -32,14 +32,14 @@ # bmxam: for some obscur reason, order 3 and 5 lead to different sha1sum # when running in standard mode or in test mode... - for degree_export in (1, 2, 4) - U_export = TrialFESpace(FunctionSpace(:Lagrange, degree_export), mesh) - basename = "write_vtk_lagrange_deg$(degree_export)" + for mesh_degree in (1, 2, 4) + basename = "write_vtk_lagrange_deg$(mesh_degree)" Bcube.write_vtk_lagrange( joinpath(tempdir, basename), vars, - mesh, - U_export; + mesh; + mesh_degree, + discontinuous = false, vtkversion = v"1.0", ) @@ -52,13 +52,13 @@ quad = Quadrature(4) dΩ = Measure(CellDomain(mesh), quad) vars["umean"] = Bcube.cell_mean(u, dΩ) - U_export = TrialFESpace(FunctionSpace(:Lagrange, 4), mesh) basename = "write_vtk_lagrange_deg4_with_mean" Bcube.write_vtk_lagrange( joinpath(tempdir, basename), vars, - mesh, - U_export; + mesh; + mesh_degree = 4, + discontinuous = false, vtkversion = v"1.0", ) @@ -66,13 +66,13 @@ fname = basename * ".vtu" @test fname2sum[fname] == bytes2hex(open(sha1, joinpath(tempdir, fname))) - U_export = TrialFESpace(FunctionSpace(:Lagrange, 4), mesh, :discontinuous) basename = "write_vtk_lagrange_deg4_dg_with_mean" Bcube.write_vtk_lagrange( joinpath(tempdir, basename), vars, - mesh, - U_export; + mesh; + mesh_degree = 4, + discontinuous = true, vtkversion = v"1.0", )