Add 3d functionality to Trixi2Img

Project.toml

@@ -5,6 +5,7 @@ version = "0.2.1-pre"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -13,6 +14,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 ArgParse = "1.1"
+EllipsisNotation = "1.0.0"
 GR = "0.51, 0.52"
 Glob = "1.3"
 HDF5 = "0.13"

src/convert.jl

@@ -27,7 +27,8 @@ julia> trixi2img("out/solution_000*.h5")
 """
 function trixi2img(filename::AbstractString...;
                    format=:png, variables=[], verbose=false, grid_lines=false,
-                   output_directory=".", nvisnodes=nothing, max_supported_level=11)
+                   output_directory=".", nvisnodes=nothing, max_supported_level=11,
+                   slice_axis=:x, slice_axis_intersect=0)
   # Reset timer
   reset_timer!()
 
@@ -65,11 +66,11 @@ function trixi2img(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)
+                         leaf_cells, coordinates, levels, ndims_) = read_meshfile(meshfile)
 
     # Read data
     verbose && println("| Reading data file...")
-    @timeit "read data" labels, unstructured_data, n_nodes, time = read_datafile(filename)
+    @timeit "read data" labels, unstructured_data, n_nodes, time = read_datafile(filename, ndims_)
 
     # Determine resolution for data interpolation
     max_level = maximum(levels)
@@ -92,6 +93,13 @@ function trixi2img(filename::AbstractString...;
     # level-0-cell) that contains the number of visualization nodes for any
     # refinement level to visualize on an equidistant grid
 
+    if ndims_ == 3
+      # convert 3d unstructured data to 2d slice
+      unstructured_data, coordinates, levels, center_level_0 = unstructured_2d_to_3d(
+          unstructured_data, coordinates, levels, length_level_0, center_level_0,
+          slice_axis, slice_axis_intersect)
+    end
+
     # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
     n_elements = length(levels)
     normalized_coordinates = similar(coordinates)
@@ -180,4 +188,3 @@ function trixi2img(filename::AbstractString...;
   print_timer()
   println()
 end
-

src/interpolate.jl

@@ -40,6 +40,113 @@ function cell2node(cell_centered_data::AbstractArray{Float64})
 end
 
 
+# Convert 3d unstructured data to 2d slice.
+# Additional to the new unstructured data updated coordinates and levels are returned.
+function unstructured_2d_to_3d(unstructured_data::AbstractArray{Float64},
+                               coordinates::AbstractArray{Float64},
+                               levels::AbstractArray{Int}, length_level_0::Float64,
+                               center_level_0::AbstractArray{Float64},
+                               slice_axis, slice_axis_intersect)
+
+  dimensions = Dict(
+    :x => 1,
+    :y => 2,
+    :z => 3
+  )
+  if !haskey(dimensions, slice_axis)
+    error("illegal dimension")
+  end
+  slice_axis_dimension = dimensions[slice_axis]
+  other_dimensions = [1, 2, 3][1:end .!= slice_axis_dimension]
+
+  # Extract data shape information
+  n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data)
+
+  # Get node coordinates for DG locations on reference element
+  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)
+
+  # New unstructured data has one dimension less.
+  # The redundant element ids are removed later.
+  new_unstructured_data = similar(unstructured_data[1, ..])
+
+  # Declare new empty arrays to fill in new coordinates and levels
+  new_coordinates = Array{Float64}(undef, 2, 0)
+  new_levels = Array{eltype(levels)}(undef, 0)
+
+  # Counter for new element ids
+  new_id = 0
+
+  # Save vandermonde matrices in a Dict to prevent redundant generation
+  vandermonde_to_2d = Dict()
+
+  # Limits of domain in slice_axis dimension
+  lower_limit = center_level_0[slice_axis_dimension] - length_level_0 / 2
+  upper_limit = center_level_0[slice_axis_dimension] + length_level_0 / 2
+
+  if slice_axis_intersect < lower_limit || slice_axis_intersect > upper_limit
+    error("slice_axis_intersect outside of domain")
+  end
+
+  for element_id in 1:n_elements
+    # Distance from center to border of this element (half the length)
+    element_length = length_level_0 / 2^levels[element_id]
+    first_coordinate = coordinates[:, element_id] .- element_length / 2
+    last_coordinate = coordinates[:, element_id] .+ element_length / 2
+
+    # Check if slice plane and current element intersect
+    # The upper limit check is needed because of the > in the first check
+    if (first_coordinate[slice_axis_dimension] <= slice_axis_intersect &&
+          last_coordinate[slice_axis_dimension] > slice_axis_intersect) ||
+        (slice_axis_intersect == upper_limit &&
+          last_coordinate[slice_axis_dimension] == upper_limit)
+      # This element is of interest
+      new_id += 1
+
+      # Add element to new coordinates and levels
+      new_coordinates = hcat(new_coordinates, coordinates[other_dimensions, element_id])
+      push!(new_levels, levels[element_id])
+
+      # Construct vandermonde matrix (or load from Dict if possible)
+      normalized_intersect =
+          (slice_axis_intersect - first_coordinate[slice_axis_dimension]) /
+          element_length * 2 - 1
+
+      if haskey(vandermonde_to_2d, normalized_intersect)
+        vandermonde = vandermonde_to_2d[normalized_intersect]
+      else
+        # Generate vandermonde matrix to interpolate values at nodes_in to one value
+        vandermonde = polynomial_interpolation_matrix(nodes_in, [normalized_intersect])
+        vandermonde_to_2d[normalized_intersect] = vandermonde
+      end
+
+      # 1D interpolation to specified slice plane
+      for i in 1:n_nodes_in
+        for ii in 1:n_nodes_in
+          for v in 1:n_variables
+            if slice_axis == :x
+              data = unstructured_data[:, i, ii, element_id, v]
+            elseif slice_axis == :y
+              data = unstructured_data[i, :, ii, element_id, v]
+            elseif slice_axis == :z
+              data = unstructured_data[i, ii, :, element_id, v]
+            end
+            value = vandermonde * data
+            new_unstructured_data[i, ii, new_id, v] = value[1]
+          end
+        end
+      end
+    end
+  end
+
+  # Remove redundant element ids
+  unstructured_data = new_unstructured_data[:, :, 1:new_id, :]
+
+  center_level_0 = center_level_0[other_dimensions]
+
+  return unstructured_data, new_coordinates, new_levels, center_level_0
+end
+
+
 # Interpolate unstructured DG data to structured data (cell-centered)
 function unstructured2structured(unstructured_data::AbstractArray{Float64},
                                  normalized_coordinates::AbstractArray{Float64},
@@ -241,4 +348,3 @@ function calc_vertices(coordinates::AbstractArray{Float64, 2},
 
   return x, y
 end
-

src/io.jl

@@ -1,3 +1,5 @@
+using EllipsisNotation
+
 # Use data file to extract mesh filename from attributes
 function extract_mesh_filename(filename::String)
   # Open file for reading
@@ -12,7 +14,6 @@ end
 
 # Read in mesh file and return relevant data
 function read_meshfile(filename::String)
-  n_children_per_cell = 2^ndim
   # Open file for reading
   h5open(filename, "r") do file
     # Extract basic information
@@ -21,13 +22,14 @@ function read_meshfile(filename::String)
     else
       ndims_ = read(attrs(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released
     end
+    n_children_per_cell = 2^ndims_
     n_cells = read(attrs(file)["n_cells"])
     n_leaf_cells = read(attrs(file)["n_leaf_cells"])
     center_level_0 = read(attrs(file)["center_level_0"])
     length_level_0 = read(attrs(file)["length_level_0"])
 
     # Extract coordinates, levels, child cells
-    coordinates = Array{Float64}(undef, ndim, n_cells)
+    coordinates = Array{Float64}(undef, ndims_, n_cells)
     coordinates .= read(file["coordinates"])
     levels = Array{Int}(undef, n_cells)
     levels .= read(file["levels"])
@@ -50,12 +52,12 @@ function read_meshfile(filename::String)
     coordinates = coordinates[:, leaf_cells]
     levels = levels[leaf_cells]
 
-    return center_level_0, length_level_0, leaf_cells, coordinates, levels
+    return center_level_0, length_level_0, leaf_cells, coordinates, levels, ndims_
   end
 end
 
 
-function read_datafile(filename::String)
+function read_datafile(filename::String, ndims::Int64)
   # Open file for reading
   h5open(filename, "r") do file
     # Extract basic information
@@ -76,14 +78,19 @@ function read_datafile(filename::String)
 
     # Extract data arrays
     n_nodes = polydeg + 1
-    data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
+
+    if ndims == 3
+      # Read 3d data
+      data = Array{Float64}(undef, n_nodes, n_nodes, n_nodes, n_elements, n_variables)
+    else
+      # Read 2d data
+      data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
+    end
     for v = 1:n_variables
       vardata = read(file["variables_$v"])
-      #=@show minimum(vardata), maximum(vardata)=#
-      @views data[:, :, :, v][:] .= vardata
+      @views data[.., v][:] .= vardata
     end
 
     return labels, data, n_nodes, time
   end
 end
-

test/runtests.jl

@@ -1 +1,2 @@
 include("test_2d.jl")
+include("test_3d.jl")

test/test_3d.jl

@@ -0,0 +1,72 @@
+module Test3D
+
+using Test
+using Trixi2Img
+
+include("test_trixi2img.jl")
+
+# Start with a clean environment: remove Trixi output directory if it exists
+outdir = "out"
+isdir(outdir) && rm(outdir, recursive=true)
+
+
+@testset "3D" begin
+  run_trixi(joinpath("3d", "parameters.toml"), n_steps_max=1)
+
+  @testset "uniform mesh as PNG" begin
+    @testset "default slice" begin
+      test_trixi2img("solution_000000.h5", outdir)
+        # To be able to compare hashes, we would need to make sure that exactly the same versions
+        # of all required libraries are installed everywhere, which does not seem to be a good option
+        # right now.
+    end
+
+    @testset "x-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:x, slice_axis_intersect=0.6)
+    end
+
+    @testset "y-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intersect=-0.35)
+    end
+
+    @testset "z-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intersect=-0.89)
+    end
+
+    @testset "negative border slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intersect=-1)
+    end
+
+    @testset "positive border slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intersect=1)
+    end
+
+    @testset "outside negative border slice" begin
+      @test_throws ErrorException("slice_axis_intersect outside of domain") trixi2img(
+          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
+          slice_axis=:x, slice_axis_intersect=-1.01)
+    end
+
+    @testset "outside positive border slice" begin
+      @test_throws ErrorException("slice_axis_intersect outside of domain") trixi2img(
+          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
+          slice_axis=:y, slice_axis_intersect=1.005)
+    end
+  end
+
+  # PDF tests do not work on Windows
+  if !Sys.iswindows()
+    @testset "uniform mesh as PDF with grid lines" begin
+      test_trixi2img("solution_000000.h5", outdir,
+          format=:pdf, grid_lines=true)
+          # To be able to compare hashes, we would need to make sure that exactly the same versions
+          # of all required libraries are installed everywhere, which does not seem to be a good option
+          # right now.
+    end
+  end
+end
+
+# Clean up afterwards: delete Trixi output directory
+@test_nowarn rm(outdir, recursive=true)
+
+end