diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1041cd6..5108a29 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,7 @@ jobs: fail-fast: false matrix: version: - - '1.5' + - '1.6' # - 'nightly' os: - ubuntu-latest diff --git a/.gitignore b/.gitignore index 7cc8281..f1677ce 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,8 @@ *.mp4 *.mem *.vtu +*.vts +*.vti *.pvd *.avi *.ogv diff --git a/Project.toml b/Project.toml index e89ffaa..b17a1c4 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] @@ -17,5 +18,6 @@ HDF5 = "0.14, 0.15" ProgressMeter = "1.3" StaticArrays = "0.12, 1.0" TimerOutputs = "0.5" +Trixi = "0.3.31" WriteVTK = "1.7" -julia = "1.5" +julia = "1.6" diff --git a/src/Trixi2Vtk.jl b/src/Trixi2Vtk.jl index 280f8d4..8add08e 100644 --- a/src/Trixi2Vtk.jl +++ b/src/Trixi2Vtk.jl @@ -6,10 +6,11 @@ using HDF5: h5open, attributes, haskey using ProgressMeter: @showprogress, Progress, next! using StaticArrays: SVector using TimerOutputs +using Trixi: Trixi, transfinite_mapping, coordinates2mapping, polynomial_interpolation_matrix, + gauss_lobatto_nodes_weights using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, vtk_save, paraview_collection # Include all top-level submodule files -include("interpolation.jl") include("interpolate.jl") include("io.jl") include("pointlocators.jl") diff --git a/src/convert.jl b/src/convert.jl index 5dafbd7..67faac0 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -15,8 +15,13 @@ Convert Trixi-generated output files to VTK files (VTU or VTI). - `pvd`: Use this filename to store PVD file (instead of auto-detecting name). Note that only the name will be used (directory and file extension are ignored). - `output_directory`: Output directory where generated files are stored. -- `nvisnodes`: Number of visualization nodes per element (default: twice the number of DG nodes). +- `nvisnodes`: Number of visualization nodes per element. + (default: number of DG nodes for CurvedMesh or UnstructuredQuadMesh, + twice the number of DG nodes for TreeMesh). A value of `0` (zero) uses the number of nodes in the DG elements. +- `save_celldata`: Boolean value to determine if cell-based data should be saved. + (the default `nothing` is converted to `false` + for `CurvedMesh`/`UnstructuredQuadMesh` and `true` for `TreeMesh`) # Examples ```julia @@ -26,7 +31,7 @@ julia> trixi2vtk("out/solution_000*.h5") """ function trixi2vtk(filename::AbstractString...; format=:vtu, verbose=false, hide_progress=false, pvd=nothing, - output_directory=".", nvisnodes=nothing) + output_directory=".", nvisnodes=nothing, save_celldata=nothing) # Reset timer reset_timer!() @@ -38,6 +43,9 @@ function trixi2vtk(filename::AbstractString...; for pattern in filename append!(filenames, glob(pattern)) end + if isempty(filenames) + error("no such file(s): ", join(filename, ", ")) + end # Ensure valid format if !(format in (:vtu, :vti)) @@ -97,8 +105,13 @@ function trixi2vtk(filename::AbstractString...; # Read mesh verbose && println("| Reading mesh file...") - @timeit "read mesh" (center_level_0, length_level_0, - leaf_cells, coordinates, levels) = read_meshfile(meshfile) + @timeit "read mesh" mesh = Trixi.load_mesh_serial(meshfile; n_cells_max=0, RealT=Float64) + + if save_celldata === nothing + # If no value for `save_celldata` is specified, + # use true for TreeMesh and false for CurvedMesh or UnstructuredQuadMesh + save_celldata = isa(mesh, Trixi.TreeMesh) + end # Read data only if it is a data file if is_datafile @@ -106,21 +119,10 @@ function trixi2vtk(filename::AbstractString...; @timeit "read data" (labels, data, n_elements, n_nodes, element_variables, time) = read_datafile(filename) - # Check if dimensions match - if length(leaf_cells) != n_elements - error("number of elements in '$(filename)' do not match number of leaf cells in " * - "'$(meshfile)' " * - "(did you forget to clean your 'out/' directory between different runs?)") - end + assert_cells_elements(n_elements, mesh, filename, meshfile) # Determine resolution for data interpolation - if nvisnodes == nothing - n_visnodes = 2 * n_nodes - elseif nvisnodes == 0 - n_visnodes = n_nodes - else - n_visnodes = nvisnodes - end + n_visnodes = get_default_nvisnodes(nvisnodes, n_nodes, mesh) else # If file is a mesh file, do not interpolate data n_visnodes = 1 @@ -130,30 +132,23 @@ function trixi2vtk(filename::AbstractString...; mkpath(output_directory) # Build VTK grids - vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), coordinates, levels, center_level_0, - length_level_0, n_visnodes, verbose, + vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, n_visnodes, verbose, output_directory, is_datafile, filename) # Interpolate data if is_datafile verbose && println("| Interpolating data...") @timeit "interpolate data" interpolated_data = interpolate_data(Val(format), - data, coordinates, levels, - center_level_0, - length_level_0, + data, mesh, n_visnodes, verbose) end # Add data to file verbose && println("| Adding data to VTK file...") @timeit "add data to VTK file" begin - # Add cell/element data to celldata VTK file - verbose && println("| | cell_ids...") - @timeit "cell_ids" vtk_celldata["cell_ids"] = leaf_cells - verbose && println("| | element_ids...") - @timeit "element_ids" vtk_celldata["element_ids"] = collect(1:length(leaf_cells)) - verbose && println("| | levels...") - @timeit "levels" vtk_celldata["levels"] = levels + if save_celldata + add_celldata!(vtk_celldata, mesh, verbose) + end # Only add data if it is a data file if is_datafile @@ -163,10 +158,12 @@ function trixi2vtk(filename::AbstractString...; @timeit label vtk_nodedata[label] = @views interpolated_data[:, variable_id] end - # Add element variables - for (label, variable) in element_variables - verbose && println("| | Element variable: $label...") - @timeit label vtk_celldata[label] = variable + if save_celldata + # Add element variables + for (label, variable) in element_variables + verbose && println("| | Element variable: $label...") + @timeit label vtk_celldata[label] = variable + end end end end @@ -177,8 +174,10 @@ function trixi2vtk(filename::AbstractString...; @timeit "save VTK file" vtk_save(vtk_nodedata) end - verbose && println("| Saving VTK file '$(vtk_celldata.path)'...") - @timeit "save VTK file" vtk_save(vtk_celldata) + if save_celldata + verbose && println("| Saving VTK file '$(vtk_celldata.path)'...") + @timeit "save VTK file" vtk_save(vtk_celldata) + end # Add to PVD file only if it is a datafile if !is_single_file @@ -186,7 +185,9 @@ function trixi2vtk(filename::AbstractString...; verbose && println("| Adding to PVD file...") @timeit "add VTK to PVD file" begin pvd[time] = vtk_nodedata - pvd_celldata[time] = vtk_celldata + if save_celldata + pvd_celldata[time] = vtk_celldata + end end has_data = true else @@ -206,8 +207,11 @@ function trixi2vtk(filename::AbstractString...; verbose && println("| Saving PVD file '$(pvd_filename).pvd'...") @timeit "save PVD files" vtk_save(pvd) end - verbose && println("| Saving PVD file '$(pvd_celldata_filename).pvd'...") - @timeit "save PVD files" vtk_save(pvd_celldata) + + if save_celldata + verbose && println("| Saving PVD file '$(pvd_celldata_filename).pvd'...") + @timeit "save PVD files" vtk_save(pvd_celldata) + end end verbose && println("| done.\n") @@ -215,3 +219,68 @@ function trixi2vtk(filename::AbstractString...; println() end + +function assert_cells_elements(n_elements, mesh::Trixi.TreeMesh, filename, meshfile) + # Check if dimensions match + if length(Trixi.leaf_cells(mesh.tree)) != n_elements + error("number of elements in '$(filename)' do not match number of leaf cells in " * + "'$(meshfile)' " * + "(did you forget to clean your 'out/' directory between different runs?)") + end +end + + +function assert_cells_elements(n_elements, mesh::Trixi.CurvedMesh, filename, meshfile) + # Check if dimensions match + if prod(size(mesh)) != n_elements + error("number of elements in '$(filename)' do not match number of cells in " * + "'$(meshfile)' " * + "(did you forget to clean your 'out/' directory between different runs?)") + end +end + + +function assert_cells_elements(n_elements, mesh::Trixi.UnstructuredQuadMesh, filename, meshfile) + # Check if dimensions match + if length(mesh) != n_elements + error("number of elements in '$(filename)' do not match number of cells in " * + "'$(meshfile)' " * + "(did you forget to clean your 'out/' directory between different runs?)") + end +end + + +function get_default_nvisnodes(nvisnodes, n_nodes, mesh::Trixi.TreeMesh) + if nvisnodes === nothing + return 2 * n_nodes + elseif nvisnodes == 0 + return n_nodes + else + return nvisnodes + end +end + + +function get_default_nvisnodes(nvisnodes, n_nodes, mesh::Union{Trixi.CurvedMesh, Trixi.UnstructuredQuadMesh}) + if nvisnodes === nothing || nvisnodes == 0 + return n_nodes + else + return nvisnodes + end +end + + +function add_celldata!(vtk_celldata, mesh::Trixi.TreeMesh, verbose) + @timeit "add data to VTK file" begin + leaf_cells = Trixi.leaf_cells(mesh.tree) + # Add cell/element data to celldata VTK file + verbose && println("| | cell_ids...") + @timeit "cell_ids" vtk_celldata["cell_ids"] = leaf_cells + verbose && println("| | element_ids...") + @timeit "element_ids" vtk_celldata["element_ids"] = collect(1:length(leaf_cells)) + verbose && println("| | levels...") + @timeit "levels" vtk_celldata["levels"] = mesh.tree.levels + end + + return vtk_celldata +end diff --git a/src/interpolate.jl b/src/interpolate.jl index 68a46d9..5696fa7 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -1,14 +1,29 @@ # Interpolate data from input format to desired output format (vtu version) -function interpolate_data(::Val{:vtu}, input_data, coordinates, levels, - center_level_0, length_level_0, n_visnodes, verbose) - return raw2visnodes(input_data, n_visnodes) +function interpolate_data(::Val{:vtu}, input_data, mesh::Trixi.TreeMesh, n_visnodes, verbose) + # Calculate equidistant output nodes (visualization nodes) + dx = 2 / n_visnodes + nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_visnodes)) + + return raw2interpolated(input_data, nodes_out) +end + + +# Interpolate data from input format to desired output format (CurvedMesh or UnstructuredQuadMesh version) +function interpolate_data(::Val{:vtu}, input_data, + mesh::Union{Trixi.CurvedMesh, Trixi.UnstructuredQuadMesh}, + n_visnodes, verbose) + # Calculate equidistant output nodes + nodes_out = collect(range(-1, 1, length=n_visnodes)) + + return raw2interpolated(input_data, nodes_out) end # Interpolate data from input format to desired output format (vti version) -function interpolate_data(::Val{:vti}, input_data, coordinates, levels, - center_level_0, length_level_0, n_visnodes, verbose) +function interpolate_data(::Val{:vti}, input_data, mesh, n_visnodes, verbose) + coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]² normalized_coordinates = similar(coordinates) for element_id in axes(coordinates, 2) @@ -130,13 +145,14 @@ function coordinate2index(coordinate, resolution::Integer) end -# Interpolate to visualization nodes -function raw2visnodes(data_gl::AbstractArray{Float64}, n_visnodes::Int) +# Interpolate to specified output nodes +function raw2interpolated(data_gl::AbstractArray{Float64}, nodes_out) # Extract number of spatial dimensions ndims_ = ndims(data_gl) - 2 # Extract data shape information n_nodes_in = size(data_gl, 1) + n_nodes_out = length(nodes_out) n_elements = size(data_gl, ndims_ + 1) n_variables = size(data_gl, ndims_ + 2) @@ -144,13 +160,11 @@ function raw2visnodes(data_gl::AbstractArray{Float64}, n_visnodes::Int) nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) # Calculate Vandermonde matrix - dx = 2 / n_visnodes - nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_visnodes)) vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out) if ndims_ == 2 # Create output data structure - data_vis = Array{Float64}(undef, n_visnodes, n_visnodes, n_elements, n_variables) + data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_elements, n_variables) # For each variable, interpolate element data and store to global data structure for v in 1:n_variables @@ -161,12 +175,12 @@ function raw2visnodes(data_gl::AbstractArray{Float64}, n_visnodes::Int) for element_id in 1:n_elements @views data_vis[:, :, element_id, v] .= reshape( interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1), - n_visnodes, n_visnodes) + n_nodes_out, n_nodes_out) end end elseif ndims_ == 3 # Create output data structure - data_vis = Array{Float64}(undef, n_visnodes, n_visnodes, n_visnodes, n_elements, n_variables) + data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_nodes_out, n_elements, n_variables) # For each variable, interpolate element data and store to global data structure for v in 1:n_variables @@ -178,7 +192,7 @@ function raw2visnodes(data_gl::AbstractArray{Float64}, n_visnodes::Int) for element_id in 1:n_elements @views data_vis[:, :, :, element_id, v] .= reshape( interpolate_nodes(reshaped_data[:, :, :, :, element_id], vandermonde, 1), - n_visnodes, n_visnodes, n_visnodes) + n_nodes_out, n_nodes_out, n_nodes_out) end end else @@ -186,6 +200,62 @@ function raw2visnodes(data_gl::AbstractArray{Float64}, n_visnodes::Int) end # Return as one 1D array for each variable - return reshape(data_vis, n_visnodes^ndims_ * n_elements, n_variables) + return reshape(data_vis, n_nodes_out^ndims_ * n_elements, n_variables) end +# Interpolate data using the given Vandermonde matrix and return interpolated values (2D version). +function interpolate_nodes(data_in::AbstractArray{T, 3}, + vandermonde, n_vars) where T + n_nodes_out = size(vandermonde, 1) + data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out) + interpolate_nodes!(data_out, data_in, vandermonde, n_vars) +end + +function interpolate_nodes!(data_out::AbstractArray{T, 3}, data_in::AbstractArray{T, 3}, + vandermonde, n_vars) where T + n_nodes_out = size(vandermonde, 1) + n_nodes_in = size(vandermonde, 2) + + for j in 1:n_nodes_out + for i in 1:n_nodes_out + for v in 1:n_vars + acc = zero(eltype(data_out)) + for jj in 1:n_nodes_in + for ii in 1:n_nodes_in + acc += vandermonde[i, ii] * data_in[v, ii, jj] * vandermonde[j, jj] + end + end + data_out[v, i, j] = acc + end + end + end + + return data_out +end + + +# Interpolate data using the given Vandermonde matrix and return interpolated values (3D version). +function interpolate_nodes(data_in::AbstractArray{T, 4}, + vandermonde, n_vars) where T + n_nodes_out = size(vandermonde, 1) + data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out, n_nodes_out) + interpolate_nodes!(data_out, data_in, vandermonde, n_vars) +end + +function interpolate_nodes!(data_out::AbstractArray{T, 4}, data_in::AbstractArray{T, 4}, + vandermonde, n_vars) where T + n_nodes_out = size(vandermonde, 1) + n_nodes_in = size(vandermonde, 2) + + for k in 1:n_nodes_out, j in 1:n_nodes_out, i in 1:n_nodes_out + for v in 1:n_vars + acc = zero(eltype(data_out)) + for kk in 1:n_nodes_in, jj in 1:n_nodes_in, ii in 1:n_nodes_in + acc += vandermonde[i, ii] * vandermonde[j, jj] * vandermonde[k, kk] * data_in[v, ii, jj, kk] + end + data_out[v, i, j, k] = acc + end + end + + return data_out +end diff --git a/src/interpolation.jl b/src/interpolation.jl deleted file mode 100644 index 9c204ef..0000000 --- a/src/interpolation.jl +++ /dev/null @@ -1,419 +0,0 @@ - -# Interpolate data using the given Vandermonde matrix and return interpolated values (1D version). -function interpolate_nodes(data_in::AbstractArray{T, 2}, - vandermonde::AbstractArray{T, 2}, n_vars::Integer) where T - n_nodes_out = size(vandermonde, 1) - n_nodes_in = size(vandermonde, 2) - data_out = zeros(eltype(data_in), n_vars, n_nodes_out) - - for i = 1:n_nodes_out - for ii = 1:n_nodes_in - for v = 1:n_vars - data_out[v, i] += vandermonde[i, ii] * data_in[v, ii] - end - end - end - - return data_out -end - - -# Interpolate data using the given Vandermonde matrix and return interpolated values (2D version). -function interpolate_nodes(data_in::AbstractArray{T, 3}, - vandermonde, n_vars) where T - n_nodes_out = size(vandermonde, 1) - data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out) - interpolate_nodes!(data_out, data_in, vandermonde, n_vars) -end - -function interpolate_nodes!(data_out::AbstractArray{T, 3}, data_in::AbstractArray{T, 3}, - vandermonde, n_vars) where T - n_nodes_out = size(vandermonde, 1) - n_nodes_in = size(vandermonde, 2) - - for j in 1:n_nodes_out - for i in 1:n_nodes_out - for v in 1:n_vars - acc = zero(eltype(data_out)) - for jj in 1:n_nodes_in - for ii in 1:n_nodes_in - acc += vandermonde[i, ii] * data_in[v, ii, jj] * vandermonde[j, jj] - end - end - data_out[v, i, j] = acc - end - end - end - - return data_out -end - - -# Interpolate data using the given Vandermonde matrix and return interpolated values (3D version). -function interpolate_nodes(data_in::AbstractArray{T, 4}, - vandermonde, n_vars) where T - n_nodes_out = size(vandermonde, 1) - data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out, n_nodes_out) - interpolate_nodes!(data_out, data_in, vandermonde, n_vars) -end - -function interpolate_nodes!(data_out::AbstractArray{T, 4}, data_in::AbstractArray{T, 4}, - vandermonde, n_vars) where T - n_nodes_out = size(vandermonde, 1) - n_nodes_in = size(vandermonde, 2) - - for k in 1:n_nodes_out, j in 1:n_nodes_out, i in 1:n_nodes_out - for v in 1:n_vars - acc = zero(eltype(data_out)) - for kk in 1:n_nodes_in, jj in 1:n_nodes_in, ii in 1:n_nodes_in - acc += vandermonde[i, ii] * vandermonde[j, jj] * vandermonde[k, kk] * data_in[v, ii, jj, kk] - end - data_out[v, i, j, k] = acc - end - end - - return data_out -end - - -# Calculate the Dhat matrix -function calc_dhat(nodes, weights) - n_nodes = length(nodes) - dhat = polynomial_derivative_matrix(nodes) - dhat = transpose(dhat) - - for n = 1:n_nodes, j = 1:n_nodes - dhat[j, n] *= -weights[n] / weights[j] - end - - return dhat -end - - -# Calculate the Dsplit matrix for split-form differentiation: dplit = 2D - M⁻¹B -function calc_dsplit(nodes, weights) - # Start with 2 x the normal D matrix - dsplit = polynomial_derivative_matrix(nodes) - dsplit = 2 .* dsplit - - # Modify to account for - dsplit[1, 1] += 1/weights[1] - dsplit[end, end] -= 1/weights[end] - - return dsplit -end - - -# Calculate the polynomial derivative matrix D -function polynomial_derivative_matrix(nodes) - n_nodes = length(nodes) - d = zeros(n_nodes, n_nodes) - wbary = barycentric_weights(nodes) - - for i = 1:n_nodes, j = 1:n_nodes - if j != i - d[i, j] = wbary[j] / wbary[i] * 1 / (nodes[i] - nodes[j]) - d[i, i] -= d[i, j] - end - end - - return d -end - - -# Calculate and interpolation matrix (Vandermonde matrix) between two given sets of nodes -function polynomial_interpolation_matrix(nodes_in, nodes_out) - n_nodes_in = length(nodes_in) - n_nodes_out = length(nodes_out) - wbary_in = barycentric_weights(nodes_in) - vdm = zeros(n_nodes_out, n_nodes_in) - - for k = 1:n_nodes_out - match = false - for j = 1:n_nodes_in - if isapprox(nodes_out[k], nodes_in[j], rtol=eps()) - match = true - vdm[k, j] = 1 - end - end - - if match == false - s = 0.0 - for j = 1:n_nodes_in - t = wbary_in[j] / (nodes_out[k] - nodes_in[j]) - vdm[k, j] = t - s += t - end - for j = 1:n_nodes_in - vdm[k, j] = vdm[k, j] / s - end - end - end - - return vdm -end - - -# Calculate the barycentric weights for a given node distribution. -function barycentric_weights(nodes) - n_nodes = length(nodes) - weights = ones(n_nodes) - - for j = 2:n_nodes, k = 1:(j-1) - weights[k] *= nodes[k] - nodes[j] - weights[j] *= nodes[j] - nodes[k] - end - - for j = 1:n_nodes - weights[j] = 1 / weights[j] - end - - return weights -end - - -# Calculate Lhat. -function calc_lhat(x::Float64, nodes, weights) - n_nodes = length(nodes) - wbary = barycentric_weights(nodes) - - lhat = lagrange_interpolating_polynomials(x, nodes, wbary) - - for i = 1:n_nodes - lhat[i] /= weights[i] - end - - return lhat -end - - -# Calculate Lagrange polynomials for a given node distribution. -function lagrange_interpolating_polynomials(x::Float64, nodes, wbary) - n_nodes = length(nodes) - polynomials = zeros(n_nodes) - - for i = 1:n_nodes - if isapprox(x, nodes[i], rtol=eps(x)) - polynomials[i] = 1 - return polynomials - end - end - - for i = 1:n_nodes - polynomials[i] = wbary[i] / (x - nodes[i]) - end - total = sum(polynomials) - - for i = 1:n_nodes - polynomials[i] /= total - end - - return polynomials -end - - -# From FLUXO (but really from blue book by Kopriva) -function gauss_lobatto_nodes_weights(n_nodes::Integer) - # From Kopriva's book - n_iterations = 10 - tolerance = 1e-15 - - # Initialize output - nodes = zeros(n_nodes) - weights = zeros(n_nodes) - - # Get polynomial degree for convenience - N = n_nodes - 1 - - # Calculate values at boundary - nodes[1] = -1.0 - nodes[end] = 1.0 - weights[1] = 2 / (N * (N + 1)) - weights[end] = weights[1] - - # Calculate interior values - if N > 1 - cont1 = pi/N - cont2 = 3/(8 * N * pi) - - # Use symmetry -> only left side is computed - for i in 1:(div(N + 1, 2) - 1) - # Calculate node - # Initial guess for Newton method - nodes[i+1] = -cos(cont1*(i+0.25) - cont2/(i+0.25)) - - # Newton iteration to find root of Legendre polynomial (= integration node) - for k in 0:n_iterations - q, qder, _ = calc_q_and_l(N, nodes[i+1]) - dx = -q/qder - nodes[i+1] += dx - if abs(dx) < tolerance * abs(nodes[i+1]) - break - end - end - - # Calculate weight - _, _, L = calc_q_and_l(N, nodes[i+1]) - weights[i+1] = weights[1] / L^2 - - # Set nodes and weights according to symmetry properties - nodes[N+1-i] = -nodes[i+1] - weights[N+1-i] = weights[i+1] - end - end - - # If odd number of nodes, set center node to origin (= 0.0) and calculate weight - if n_nodes % 2 == 1 - _, _, L = calc_q_and_l(N, 0) - nodes[div(N, 2) + 1] = 0.0 - weights[div(N, 2) + 1] = weights[1] / L^2 - end - - return nodes, weights -end - - -# From FLUXO (but really from blue book by Kopriva) -function calc_q_and_l(N::Integer, x::Float64) - L_Nm2 = 1.0 - L_Nm1 = x - Lder_Nm2 = 0.0 - Lder_Nm1 = 1.0 - - local L - for i in 2:N - L = ((2 * i - 1) * x * L_Nm1 - (i - 1) * L_Nm2) / i - Lder = Lder_Nm2 + (2 * i - 1) * L_Nm1 - L_Nm2 = L_Nm1 - L_Nm1 = L - Lder_Nm2 = Lder_Nm1 - Lder_Nm1 = Lder - end - - q = (2 * N + 1)/(N + 1) * (x * L - L_Nm2) - qder = (2 * N + 1) * L - - return q, qder, L -end -calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x)) - - -# From FLUXO (but really from blue book by Kopriva) -function gauss_nodes_weights(n_nodes::Integer) - # From Kopriva's book - n_iterations = 10 - tolerance = 1e-15 - - # Initialize output - nodes = ones(n_nodes) * 1000 - weights = zeros(n_nodes) - - # Get polynomial degree for convenience - N = n_nodes - 1 - if N == 0 - nodes .= 0.0 - weights .= 2.0 - return nodes, weights - elseif N == 1 - nodes[1] = -sqrt(1/3) - nodes[end] = -nodes[1] - weights .= 1.0 - return nodes, weights - else # N > 1 - # Use symmetry property of the roots of the Legendre polynomials - for i in 0:(div(N + 1, 2) - 1) - # Starting guess for Newton method - nodes[i+1] = -cos(pi / (2 * N + 2) * (2 * i + 1)) - - # Newton iteration to find root of Legendre polynomial (= integration node) - for k in 0:n_iterations - poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1]) - dx = -poly / deriv - nodes[i+1] += dx - if abs(dx) < tolerance * abs(nodes[i+1]) - break - end - end - - # Calculate weight - poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1]) - weights[i+1] = (2 * N + 3) / ((1 - nodes[i+1]^2) * deriv^2) - - # Set nodes and weights according to symmetry properties - nodes[N+1-i] = -nodes[i+1] - weights[N+1-i] = weights[i+1] - end - - # If odd number of nodes, set center node to origin (= 0.0) and calculate weight - if n_nodes % 2 == 1 - poly, deriv = legendre_polynomial_and_derivative(N + 1, 0.0) - nodes[div(N, 2) + 1] = 0.0 - weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2 - end - - return nodes, weights - end -end - - -# From FLUXO (but really from blue book by Kopriva) -function legendre_polynomial_and_derivative(N::Int, x::Real) - if N == 0 - poly = 1.0 - deriv = 0.0 - elseif N == 1 - poly = convert(Float64, x) - deriv = 1.0 - else - poly_Nm2 = 1.0 - poly_Nm1 = convert(Float64, x) - deriv_Nm2 = 0.0 - deriv_Nm1 = 1.0 - - poly = 0.0 - deriv = 0.0 - for i in 2:N - poly = ((2*i-1) * x * poly_Nm1 - (i-1) * poly_Nm2) / i - deriv=deriv_Nm2 + (2*i-1)*poly_Nm1 - poly_Nm2=poly_Nm1 - poly_Nm1=poly - deriv_Nm2=deriv_Nm1 - deriv_Nm1=deriv - end - end - - # Normalize - poly = poly * sqrt(N+0.5) - deriv = deriv * sqrt(N+0.5) - - return poly, deriv -end - - -# Calculate Legendre vandermonde matrix and its inverse -function vandermonde_legendre(nodes, N) - n_nodes = length(nodes) - n_modes = N + 1 - vandermonde = zeros(n_nodes, n_modes) - - for i in 1:n_nodes - for m in 1:n_modes - vandermonde[i, m], _ = legendre_polynomial_and_derivative(m-1, nodes[i]) - end - end - # for very high polynomial degree, this is not well conditioned - inverse_vandermonde = inv(vandermonde) - return vandermonde, inverse_vandermonde -end -vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1) - - -# Convert nodal to modal representation -function nodal2modal(data_in, vandermonde) - return interpolate_nodes(data_in, vandermonde, 1) -end - -function nodal2modal!(data_out, data_in, vandermonde) - return interpolate_nodes!(data_out, data_in, vandermonde, 1) -end - diff --git a/src/io.jl b/src/io.jl index 1ce28d6..0c5f254 100644 --- a/src/io.jl +++ b/src/io.jl @@ -20,47 +20,14 @@ function extract_mesh_filename(filename::String) end -# Read in mesh file and return relevant data -function read_meshfile(filename::String) - # Open file for reading - h5open(filename, "r") do file - # Extract basic information - if haskey(attributes(file), "ndims") - ndims_ = read(attributes(file)["ndims"]) - else - ndims_ = read(attributes(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released - end - n_cells = read(attributes(file)["n_cells"]) - n_leaf_cells = read(attributes(file)["n_leaf_cells"]) - center_level_0 = read(attributes(file)["center_level_0"]) - length_level_0 = read(attributes(file)["length_level_0"]) - - # Extract coordinates, levels, child cells - coordinates = Array{Float64}(undef, ndims_, n_cells) - coordinates .= read(file["coordinates"]) - levels = Array{Int}(undef, n_cells) - levels .= read(file["levels"]) - child_ids = Array{Int}(undef, 2^ndims_, n_cells) - child_ids .= read(file["child_ids"]) - - # Extract leaf cells (= cells to be plotted) and contract all other arrays accordingly - leaf_cells = similar(levels) - n_cells = 0 - for cell_id in 1:length(levels) - if sum(child_ids[:, cell_id]) > 0 - continue - end +function extract_mesh_information(mesh::Trixi.TreeMesh) + center_level_0 = mesh.tree.center_level_0 + length_level_0 = mesh.tree.length_level_0 - n_cells += 1 - leaf_cells[n_cells] = cell_id - end - leaf_cells = leaf_cells[1:n_cells] - - coordinates = coordinates[:, leaf_cells] - levels = levels[leaf_cells] + coordinates = mesh.tree.coordinates[:, Trixi.leaf_cells(mesh.tree)] + levels = mesh.tree.levels[Trixi.leaf_cells(mesh.tree)] - return center_level_0, length_level_0, leaf_cells, coordinates, levels - end + return coordinates, levels, center_level_0, length_level_0 end @@ -69,16 +36,8 @@ function read_datafile(filename::String) # Open file for reading h5open(filename, "r") do file # Extract basic information - if haskey(attributes(file), "ndims") - ndims_ = read(attributes(file)["ndims"]) - else - ndims_ = read(attributes(file)["ndim"]) - end - if haskey(attributes(file), "polydeg") - polydeg = read(attributes(file)["polydeg"]) - else - polydeg = read(attributes(file)["N"]) - end + ndims_ = read(attributes(file)["ndims"]) + polydeg = read(attributes(file)["polydeg"]) n_elements = read(attributes(file)["n_elements"]) n_variables = read(attributes(file)["n_vars"]) time = read(attributes(file)["time"]) @@ -120,4 +79,3 @@ function read_datafile(filename::String) return labels, data, n_elements, n_nodes, element_variables, time end end - diff --git a/src/vtktools.jl b/src/vtktools.jl index 2b9eb4e..7bd9778 100644 --- a/src/vtktools.jl +++ b/src/vtktools.jl @@ -1,6 +1,8 @@ # Create and return VTK grids that are ready to be filled with data (vtu version) -function build_vtk_grids(::Val{:vtu}, coordinates, levels, center_level_0, length_level_0, - n_visnodes, verbose, output_directory, is_datafile, filename) +function build_vtk_grids(::Val{:vtu}, mesh::Trixi.TreeMesh, n_visnodes, verbose, + output_directory, is_datafile, filename) + coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + # Extract number of spatial dimensions ndims_ = size(coordinates, 1) @@ -42,45 +44,181 @@ end # Create and return VTK grids that are ready to be filled with data (vti version) -function build_vtk_grids(::Val{:vti}, coordinates, levels, center_level_0, length_level_0, +function build_vtk_grids(::Val{:vti}, mesh, n_visnodes, verbose, + output_directory, is_datafile, filename) + + coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh) + + # Extract number of spatial dimensions + ndims_ = size(coordinates, 1) + + # Prepare VTK points and cells for celldata file + @timeit "prepare VTK cells" vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells( + Val(ndims_), coordinates, levels, center_level_0, length_level_0, 1) + + # Determine output file names + base, _ = splitext(splitdir(filename)[2]) + vtk_filename = joinpath(output_directory, base) + vtk_celldata_filename = vtk_filename * "_celldata" + + # Open VTK files + verbose && println("| Building VTK grid...") + if is_datafile + # Determine level-wise resolution + max_level = maximum(levels) + resolution = n_visnodes * 2^max_level + + Nx = Ny = resolution + 1 + dx = dy = length_level_0/resolution + origin = center_level_0 .- 1/2 * length_level_0 + spacing = [dx, dy] + @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, Nx, Ny, + origin=origin, + spacing=spacing) + else + vtk_nodedata = nothing + end + @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename, + vtk_celldata_points, + vtk_celldata_cells) + + return vtk_nodedata, vtk_celldata +end + + +# Create and return VTK grids that are ready to be filled with data (CurvedMesh/UnstructuredQuadMesh version) +function build_vtk_grids(::Val{:vtu}, mesh::Union{Trixi.CurvedMesh, Trixi.UnstructuredQuadMesh}, n_visnodes, verbose, output_directory, is_datafile, filename) - # Extract number of spatial dimensions - ndims_ = size(coordinates, 1) - - # Prepare VTK points and cells for celldata file - @timeit "prepare VTK cells" vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells( - Val(ndims_), coordinates, levels, center_level_0, length_level_0, 1) - - # Determine output file names - base, _ = splitext(splitdir(filename)[2]) - vtk_filename = joinpath(output_directory, base) - vtk_celldata_filename = vtk_filename * "_celldata" - - # Open VTK files - verbose && println("| Building VTK grid...") - if is_datafile - # Determine level-wise resolution - max_level = maximum(levels) - resolution = n_visnodes * 2^max_level - - Nx = Ny = resolution + 1 - dx = dy = length_level_0/resolution - origin = center_level_0 .- 1/2 * length_level_0 - spacing = [dx, dy] - @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, Nx, Ny, - origin=origin, - spacing=spacing) - else - vtk_nodedata = nothing + + @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, n_visnodes) + + # Calculate VTK points and cells + verbose && println("| Preparing VTK cells...") + if is_datafile + @timeit "prepare VTK cells (node data)" begin + vtk_points, vtk_cells = calc_vtk_points_cells(node_coordinates) end - @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename, - vtk_celldata_points, - vtk_celldata_cells) + end + + # Determine output file names + base, _ = splitext(splitdir(filename)[2]) + vtk_filename = joinpath(output_directory, base) + + # Open VTK files + verbose && println("| Building VTK grid...") + if is_datafile + @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, vtk_points, + vtk_cells) + else + vtk_nodedata = nothing + end + + vtk_celldata = nothing return vtk_nodedata, vtk_celldata end +function calc_node_coordinates(mesh::Trixi.CurvedMesh, n_visnodes) + # Extract number of spatial dimensions + ndims_ = ndims(mesh) + n_elements = prod(size(mesh)) + + nodes = range(-1, 1, length=n_visnodes) + f(args...; kwargs...) = Base.invokelatest(mesh.mapping, args...; kwargs...) + + node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., n_elements) + + return calc_node_coordinates!(node_coordinates, f, nodes, mesh) +end + + +function calc_node_coordinates(mesh::Trixi.UnstructuredQuadMesh, n_visnodes) + # Extract number of spatial dimensions + ndims_ = ndims(mesh) + n_elements = length(mesh) + + # get the uniform nodes + nodes = range(-1, 1, length=n_visnodes) + + # intialize the container for the node coordinates + node_coordinates = Array{Float64, ndims_+2}(undef, ndims_, ntuple(_ -> n_visnodes, ndims_)..., n_elements) + + # work container for the corners of elements + four_corners = zeros(eltype(mesh.corners), 4, 2) + + # loop through all elements and call the correct node coordinate constructor based on curved + for element = 1:n_elements + if mesh.element_is_curved[element] + Trixi.calc_node_coordinates!(node_coordinates, element, nodes, view(mesh.surface_curves, :, element)) + else # straight sided element + for i in 1:4, j in 1:2 + # pull the (x,y) values of these corners out of the global corners array + four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]] + end + Trixi.calc_node_coordinates!(node_coordinates, element, nodes, four_corners) + end + end + + return node_coordinates +end + + +# Calculation of the node coordinates for `CurvedMesh` in 2D +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, f, nodes, mesh) + linear_indices = LinearIndices(size(mesh)) + + # Get cell length in reference mesh + dx = 2 / size(mesh, 1) + dy = 2 / size(mesh, 2) + + for cell_y in 1:size(mesh, 2), cell_x in 1:size(mesh, 1) + element = linear_indices[cell_x, cell_y] + + # Calculate node coordinates of reference mesh + cell_x_offset = -1 + (cell_x-1) * dx + dx/2 + cell_y_offset = -1 + (cell_y-1) * dy + dy/2 + + for j in eachindex(nodes), i in eachindex(nodes) + # node_coordinates are the mapped reference node_coordinates + node_coordinates[:, i, j, element] .= f(cell_x_offset + dx/2 * nodes[i], + cell_y_offset + dy/2 * nodes[j]) + end + end + + return node_coordinates +end + + +# Calculation of the node coordinates for `CurvedMesh` in 3D +function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, f, nodes, mesh) + linear_indices = LinearIndices(size(mesh)) + + # Get cell length in reference mesh + dx = 2 / size(mesh, 1) + dy = 2 / size(mesh, 2) + dz = 2 / size(mesh, 3) + + for cell_z in 1:size(mesh, 3), cell_y in 1:size(mesh, 2), cell_x in 1:size(mesh, 1) + element = linear_indices[cell_x, cell_y, cell_z] + + # Calculate node coordinates of reference mesh + cell_x_offset = -1 + (cell_x-1) * dx + dx/2 + cell_y_offset = -1 + (cell_y-1) * dy + dy/2 + cell_z_offset = -1 + (cell_z-1) * dz + dz/2 + + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + # node_coordinates are the mapped reference node_coordinates + node_coordinates[:, i, j, k, element] .= f(cell_x_offset + dx/2 * nodes[i], + cell_y_offset + dy/2 * nodes[j], + cell_z_offset + dz/2 * nodes[k]) + end + end + + return node_coordinates +end + + # Determine and return filenames for PVD fiels function pvd_filenames(filenames, pvd, output_directory) # Determine pvd filename @@ -110,7 +248,7 @@ end # Determine filename for PVD file based on common name -function get_pvd_filename(filenames::AbstractArray) +function get_pvd_filename(filenames) filenames = getindex.(splitdir.(filenames), 2) bases = getindex.(splitext.(filenames), 1) pvd_filename = longest_common_prefix(bases) @@ -119,7 +257,7 @@ end # Determine longest common prefix -function longest_common_prefix(strings::AbstractArray) +function longest_common_prefix(strings) # Return early if array is empty if isempty(strings) return "" @@ -263,3 +401,94 @@ function calc_vtk_points_cells(::Val{3}, coordinates::AbstractMatrix{Float64}, return vtk_points, vtk_cells end + +# Convert coordinates and level information to a list of points and VTK cells for `CurvedMesh` (2D version) +function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,4}) + n_elements = size(node_coordinates, 4) + size_ = size(node_coordinates) + n_points = prod(size_[2:end]) + # Linear indices to access points by node indices and element id + linear_indices = LinearIndices(size_[2:end]) + + # Use lagrange nodes as VTK points + vtk_points = reshape(node_coordinates, (2, n_points)) + vtk_cells = Vector{MeshCell}(undef, n_elements) + + # Create cell for each element + for element in 1:n_elements + vertices = [linear_indices[1, 1, element], + linear_indices[end, 1, element], + linear_indices[end, end, element], + linear_indices[1, end, element]] + + edges = vcat(linear_indices[2:end-1, 1, element], + linear_indices[end, 2:end-1, element], + linear_indices[2:end-1, end, element], + linear_indices[1, 2:end-1, element]) + + faces = vec(linear_indices[2:end-1, 2:end-1, element]) + + point_ids = vcat(vertices, edges, faces) + vtk_cells[element] = MeshCell(VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL, point_ids) + end + + return vtk_points, vtk_cells +end + + +# Convert coordinates and level information to a list of points and VTK cells for `CurvedMesh` (3D version) +function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,5}) + n_elements = size(node_coordinates, 5) + size_ = size(node_coordinates) + n_points = prod(size_[2:end]) + # Linear indices to access points by node indices and element id + linear_indices = LinearIndices(size_[2:end]) + + # Use lagrange nodes as VTK points + vtk_points = reshape(node_coordinates, (3, n_points)) + vtk_cells = Vector{MeshCell}(undef, n_elements) + + # Create cell for each element + for element in 1:n_elements + vertices = [linear_indices[1, 1, 1, element], + linear_indices[end, 1, 1, element], + linear_indices[end, end, 1, element], + linear_indices[1, end, 1, element], + linear_indices[1, 1, end, element], + linear_indices[end, 1, end, element], + linear_indices[end, end, end, element], + linear_indices[1, end, end, element]] + + # This order doesn't make any sense. This is completely different + # from what is shown in + # https://blog.kitware.com/wp-content/uploads/2018/09/Source_Issue_43.pdf + # but this is the way it works. + edges = vcat(linear_indices[2:end-1, 1, 1, element], + linear_indices[end, 2:end-1, 1, element], + linear_indices[2:end-1, end, 1, element], + linear_indices[1, 2:end-1, 1, element], + linear_indices[2:end-1, 1, end, element], + linear_indices[end, 2:end-1, end, element], + linear_indices[2:end-1, end, end, element], + linear_indices[1, 2:end-1, end, element], + linear_indices[1, 1, 2:end-1, element], + linear_indices[end, 1, 2:end-1, element], + linear_indices[1, end, 2:end-1, element], + linear_indices[end, end, 2:end-1, element]) + + # See above + faces = vcat(vec(linear_indices[1, 2:end-1, 2:end-1, element]), + vec(linear_indices[end, 2:end-1, 2:end-1, element]), + vec(linear_indices[2:end-1, 1, 2:end-1, element]), + vec(linear_indices[2:end-1, end, 2:end-1, element]), + vec(linear_indices[2:end-1, 2:end-1, 1, element]), + vec(linear_indices[2:end-1, 2:end-1, end, element])) + + volume = vec(linear_indices[2:end-1, 2:end-1, 2:end-1, element]) + + point_ids = vcat(vertices, edges, faces, volume) + vtk_cells[element] = MeshCell(VTKCellTypes.VTK_LAGRANGE_HEXAHEDRON, point_ids) + end + + return vtk_points, vtk_cells +end diff --git a/test/Project.toml b/test/Project.toml index 012d81b..d367dfd 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl index f27e121..fb15e00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,2 +1,3 @@ include("test_2d.jl") +include("test_3d.jl") include("test_manual.jl") diff --git a/test/test_2d.jl b/test/test_2d.jl index cbf6c19..55289c1 100644 --- a/test/test_2d.jl +++ b/test/test_2d.jl @@ -9,51 +9,122 @@ include("test_trixi2vtk.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive=true) -# Create empty artifacts directory where all files that should be preserved will be stored +# Create artifacts directory where all files that should be preserved will be stored artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") -isdir(artifacts_dir) && rm(artifacts_dir, recursive=true) -mkdir(artifacts_dir) +if !isdir(artifacts_dir) + mkdir(artifacts_dir) +end @testset "2D" begin - run_trixi(joinpath("2d", "elixir_advection_extended.jl"), maxiters=1) - - @testset "uniform mesh" begin - test_trixi2vtk("solution_000000.h5", outdir, - hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"), - ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")]) - - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu") - testname = "uniform-mesh" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile)) + @testset "TreeMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("2d", "elixir_advection_extended.jl"), maxiters=1) + + @testset "uniform mesh" begin + if Sys.iswindows() + # This test fails on Windows due to globbing not working + test_trixi2vtk("solution_000000.h5", outdir, + hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"), + ("solution_000000_celldata.vtu", "9b20ba10df0d2d0fbd15916e5da0ed72ade9890b")]) + outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu") + + else + test_trixi2vtk("solution_00000*.h5", outdir, + hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"), + ("solution_000000_celldata.vtu", "9b20ba10df0d2d0fbd15916e5da0ed72ade9890b"), + ("solution_00000.pvd", "7ba2f8f1927e90ebd4209aab890c58a20acf63f4"), + ("solution_00000_celldata.pvd", "448a7130a608ed9f7e4630033b9e1338b1403f7b")]) + outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu", + "solution_00000.pvd", "solution_00000_celldata.pvd") + end + + # Store output files as artifacts to facilitate debugging of failing tests + testname = "2d-tree-mesh-uniform-mesh" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end + end + + @testset "uniform mesh with vti output" begin + if Sys.isapple() + # This test fails on MacOS due to differing binary VTI files (even though the contents match) + test_trixi2vtk("restart_000001.h5", outdir, + format=:vti) + else + test_trixi2vtk("restart_000001.h5", outdir, + hashes=[("restart_000001.vti", "664f25ab018a373774b5aad69ad3f2f5a3b21649"), + ("restart_000001_celldata.vtu", "9b20ba10df0d2d0fbd15916e5da0ed72ade9890b")], + format=:vti) + end + + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("restart_000001.vti", "restart_000001_celldata.vtu") + testname = "2d-tree-mesh-uniform-mesh-with-vti-output" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end end end - @testset "uniform mesh with vti output" begin - if Sys.isapple() - # This test fails on MacOS due to differing binary VTI files (even though the contents match) - test_trixi2vtk("restart_000001.h5", outdir, - format=:vti) - else - test_trixi2vtk("restart_000001.h5", outdir, - hashes=[("restart_000001.vti", "664f25ab018a373774b5aad69ad3f2f5a3b21649"), - ("restart_000001_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")], - format=:vti) + @testset "CurvedMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("2d", "elixir_advection_waving_flag.jl"), maxiters=1) + + @testset "waving flag" begin + test_trixi2vtk("solution_000000.h5", outdir, + hashes=[("solution_000000.vtu", "a93dbd393647627a861d890568e65598be0062f9")]) + + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("solution_000000.vtu",) + testname = "2d-curved-mesh-waving-flag" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end + end + + @testset "waving flag (supersampling)" begin + test_trixi2vtk("solution_000000.h5", outdir, nvisnodes=6, + hashes=[("solution_000000.vtu", "c6d74ab831bf4b6de2ba8cf537b6653ad611cfe7")]) + + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("solution_000000.vtu",) + testname = "2d-curved-mesh-waving-flag-supersampling" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end end + end + + @testset "UnstructuredQuadMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("2d", "elixir_euler_unstructured_quad_basic.jl"), maxiters=1) + + @testset "basic" begin + test_trixi2vtk("solution_000000.h5", outdir, + hashes=[("solution_000000.vtu", "0daedeea99d03d53b925ce5691bd7924abe88861")]) - # Store output files as artifacts to facilitate debugging of failing tests - outfiles = ("restart_000001.vti", "restart_000001_celldata.vtu") - testname = "uniform-mesh-with-vti-output" - for outfile in outfiles - println("Copying '", abspath(joinpath(outdir, outfile)), - "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), - "'...") - cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile)) + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("solution_000000.vtu",) + testname = "2d-unstructured-quad-basic" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end end end end diff --git a/test/test_3d.jl b/test/test_3d.jl new file mode 100644 index 0000000..0a79417 --- /dev/null +++ b/test/test_3d.jl @@ -0,0 +1,66 @@ +module Test3D + +using Test +using Trixi2Vtk + +include("test_trixi2vtk.jl") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +# Create artifacts directory where all files that should be preserved will be stored +artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts") +if !isdir(artifacts_dir) + mkdir(artifacts_dir) +end + + +@testset "3D" begin + @testset "TreeMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("3d", "elixir_advection_extended.jl"), maxiters=1) + + @testset "uniform mesh" begin + test_trixi2vtk("solution_000000.h5", outdir, + hashes=[("solution_000000.vtu", "6ab3aa525851187ee0839e1d670a254a66be4ad7"), + ("solution_000000_celldata.vtu", "99c782d732a4d1f6764013c3fba2cdaddf3927ab")]) + + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu") + testname = "3d-tree-mesh-uniform-mesh" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end + end + end + + @testset "CurvedMesh" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("3d", "elixir_advection_basic_curved.jl"), maxiters=1) + + @testset "basic" begin + test_trixi2vtk("solution_000000.h5", outdir, + hashes=[("solution_000000.vtu", "57c58480d8f99b7e9d365cb3fa71790db42de750")]) + + # Store output files as artifacts to facilitate debugging of failing tests + outfiles = ("solution_000000.vtu",) + testname = "3d-curved-mesh-basic" + for outfile in outfiles + println("Copying '", abspath(joinpath(outdir, outfile)), + "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)), + "'...") + cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true) + end + end + end +end + +# Clean up afterwards: delete Trixi output directory +@test_skip rm(outdir, recursive=true) + +end + diff --git a/test/test_manual.jl b/test/test_manual.jl index e5ca50c..2e7a010 100644 --- a/test/test_manual.jl +++ b/test/test_manual.jl @@ -1,10 +1,51 @@ module TestManual using Test -using Documenter using Trixi2Vtk +using Documenter + +include("test_trixi2vtk.jl") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +@testset "Manual" begin + @testset "Documenter" begin + DocMeta.setdocmeta!(Trixi2Vtk, :DocTestSetup, :(using Trixi2Vtk); recursive=true) + doctest(Trixi2Vtk, manual=false) + end + + @testset "trixi2vtk error triggers" begin + isdir(outdir) && rm(outdir, recursive=true) + run_trixi(joinpath("2d", "elixir_advection_extended.jl"), maxiters=1) + + @testset "no input file" begin + @test_throws ErrorException trixi2vtk() + end + + @testset "no such files" begin + @test_throws ErrorException trixi2vtk("this_does_not_exist") + end + + @testset "unsupported file format" begin + @test_throws ErrorException trixi2vtk(joinpath(outdir, "solution_000000.h5"); + output_directory=outdir, + format=:does_not_exist) + end + end + + @testset "trixi2vtk for mesh file" begin + @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"); output_directory=outdir) + end + + @testset "pvd_filenames" begin + @test Trixi2Vtk.pvd_filenames("", "manual", "out") == (joinpath("out", "manual"), joinpath("out", "manual_celldata")) + @test_throws ErrorException Trixi2Vtk.pvd_filenames(("a", "b"), nothing, "out") + end +end -DocMeta.setdocmeta!(Trixi2Vtk, :DocTestSetup, :(using Trixi2Vtk); recursive=true) -doctest(Trixi2Vtk, manual=false) +# Clean up afterwards: delete Trixi output directory +@test_nowarn rm(outdir, recursive=true) end # module