Feature: DGMulti support.

src/Trixi2Vtk.jl

@@ -8,7 +8,7 @@ using ProgressMeter: @showprogress, Progress, next!
 using StaticArrays: SVector
 using TimerOutputs
 using Trixi: Trixi, transfinite_mapping, coordinates2mapping, polynomial_interpolation_matrix,
-             gauss_lobatto_nodes_weights, TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh
+             gauss_lobatto_nodes_weights, TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, DGMultiMesh
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, vtk_save, paraview_collection
 
 # Include all top-level submodule files

src/convert.jl

@@ -117,7 +117,12 @@ function trixi2vtk(filename::AbstractString...;
 
     # Read mesh
     verbose && println("| Reading mesh file...")
-    @timeit "read mesh" mesh = Trixi.load_mesh_serial(meshfile; n_cells_max=0, RealT=Float64)
+    @timeit "read mesh" mesh = Trixi.load_mesh_serial(meshfile; n_cells_max=0,
+                                                      RealT=Float64)
+
+    # For DGMulti, we will reconstruct the `basis`, which is of type `RefElemData`. For 
+    # other mesh types, we do not need to do this, and will just set `basis = nothing`.
+    @timeit "reconstruct basis" basis = load_basis(meshfile, mesh)
 
     # Check compatibility of the mesh type and the output format
     if format === :vti && !(mesh isa Trixi.TreeMesh{2})
@@ -127,9 +132,11 @@ function trixi2vtk(filename::AbstractString...;
     # Read data only if it is a data file
     if is_datafile
       verbose && println("| Reading data file...")
-      @timeit "read data" (labels, data, n_elements, n_nodes,
-                           element_variables, node_variables, time) = read_datafile(filename)
 
+      @timeit "read data" (labels, data, n_elements, n_nodes,
+                           element_variables, node_variables, 
+                           time) = read_datafile(filename, mesh)
+
       assert_cells_elements(n_elements, mesh, filename, meshfile)
 
       # Determine resolution for data interpolation
@@ -156,6 +163,7 @@ function trixi2vtk(filename::AbstractString...;
     else
       # If file is a mesh file, do not interpolate data as detailed
       n_visnodes = get_default_nvisnodes_mesh(nvisnodes, mesh)
+
       # Create an "empty" node set that is unused in the mesh conversion
       node_set = Array{Float64}(undef, n_visnodes)
     end
@@ -164,16 +172,17 @@ function trixi2vtk(filename::AbstractString...;
     mkpath(output_directory)
 
     # Build VTK grids
-    vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, node_set, n_visnodes, verbose,
-                                                 output_directory, is_datafile, filename, Val(reinterpolate))
+    vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, node_set, basis,
+                                                 n_visnodes, verbose, output_directory, is_datafile, filename, Val(reinterpolate))
 
     # Interpolate data
     if is_datafile
       verbose && println("| Interpolating data...")
       if reinterpolate
         @timeit "interpolate data" interpolated_data = interpolate_data(Val(format),
                                                                         data, mesh,
-                                                                        n_visnodes, verbose)
+                                                                        basis, n_visnodes,
+                                                                        verbose)
       else # Copy the raw solution data; only works for `vtu` format
         # Extract data shape information
         ndims_ = ndims(data) - 2
@@ -220,7 +229,8 @@ function trixi2vtk(filename::AbstractString...;
               end
               @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format),
                                                                     reshape(variable, size(variable)..., 1),
-                                                                    mesh, n_visnodes, verbose)
+                                                                    mesh, basis, 
+                                                                    n_visnodes, verbose)
             else
               @timeit "interpolate data" interpolated_cell_data = reshape(variable,
                                                                           n_visnodes^ndims_ * n_elements)
@@ -314,7 +324,7 @@ function assert_cells_elements(n_elements, mesh::UnstructuredMesh2D, filename, m
 end
 
 
-function assert_cells_elements(n_elements, mesh::Union{P4estMesh, T8codeMesh}, filename, meshfile)
+function assert_cells_elements(n_elements, mesh::Union{P4estMesh, T8codeMesh, DGMultiMesh}, filename, meshfile)
   # Check if dimensions match
   if Trixi.ncells(mesh) != n_elements
     error("number of elements in '$(filename)' do not match number of cells in " *
@@ -336,7 +346,7 @@ function get_default_nvisnodes_solution(nvisnodes, n_nodes, mesh::TreeMesh)
 end
 
 function get_default_nvisnodes_solution(nvisnodes, n_nodes,
-                                        mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh})
+                                        mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, DGMultiMesh})
   if nvisnodes === nothing || nvisnodes == 0
     return n_nodes
   else
@@ -356,7 +366,7 @@ function get_default_nvisnodes_mesh(nvisnodes, mesh::TreeMesh)
 end
 
 function get_default_nvisnodes_mesh(nvisnodes,
-                                    mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh})
+                                    mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh, DGMultiMesh})
   if nvisnodes === nothing
     # for curved meshes, we need to get at least the vertices
     return 2
@@ -463,6 +473,16 @@ function add_celldata!(vtk_celldata, mesh::T8codeMesh, verbose)
   return vtk_celldata
 end
 
+function add_celldata!(vtk_celldata, mesh::DGMultiMesh, verbose)
+  @timeit "add data to VTK file" begin
+    # Add tree/element data to celldata VTK file
+    verbose && println("| | element_ids...")
+    @timeit "element_ids" vtk_celldata["element_ids"] = collect(1:Trixi.ncells(mesh))
+  end
+
+  return vtk_celldata
+end
+
 function expand_filename_patterns(patterns)
   filenames = String[]
 

src/interpolate.jl

@@ -1,6 +1,7 @@
 
 # Interpolate data from input format to desired output format (vtu version)
-function interpolate_data(::Val{:vtu}, input_data, mesh::TreeMesh, n_visnodes, verbose)
+function interpolate_data(::Val{:vtu}, input_data, mesh::TreeMesh, basis, 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))
@@ -11,7 +12,7 @@ end
 
 # Interpolate data from input format to desired output format (StructuredMesh or UnstructuredMesh2D version)
 function interpolate_data(::Val{:vtu}, input_data,
-                          mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh},
+                          mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh}, basis,
                           n_visnodes, verbose)
   # Calculate equidistant output nodes
   nodes_out = collect(range(-1, 1, length=n_visnodes))
@@ -21,8 +22,9 @@ end
 
 
 # Interpolate data from input format to desired output format (vti version)
-function interpolate_data(::Val{:vti}, input_data, mesh::TreeMesh, n_visnodes, verbose)
-  coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh)
+function interpolate_data(::Val{:vti}, input_data, mesh::TreeMesh, basis, n_visnodes, 
+                          verbose)
+coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh)
 
   # Normalize element coordinates: move center to origin and domain size to [-1, 1]²
   normalized_coordinates = similar(coordinates)
@@ -47,6 +49,33 @@ function interpolate_data(::Val{:vti}, input_data, mesh::TreeMesh, n_visnodes, v
   return structured_data
 end
 
+# Interpolate data from input format to desired output format (DGMulti version)
+function interpolate_data(::Val{:vtu}, input_data,
+                          mesh::DGMultiMesh, basis,
+                          n_visnodes, verbose)
+  if length(basis.N) > 1
+      @assert length(Set(basis.N)) == 1 "`order` must have equal elements."
+      order = first(basis.N)
+  else
+      order = basis.N
+  end
+
+  visualization_nodes = Trixi.StartUpDG.equi_nodes(basis.element_type, order)
+
+  interpolator = Trixi.StartUpDG.vandermonde(basis.element_type, order, 
+                                             visualization_nodes...) / basis.VDM
+
+  dof_per_elem, n_elements, n_variables = size(input_data)
+  result = zeros((interpolator |> size |> first, n_elements, n_variables))
+
+  for v = 1:n_variables
+    for e = 1:n_elements
+      @views result[:,e,v] = interpolator * input_data[:,e,v]
+    end
+  end
+
+  return reshape(result, (interpolator |> size |> first) * n_elements, n_variables)
+end
 
 # Interpolate unstructured DG data to structured data (cell-centered)
 function unstructured2structured(unstructured_data::AbstractArray{Float64},

src/io.jl

@@ -32,7 +32,8 @@ end
 
 
 # Read in data file and return all relevant information
-function read_datafile(filename::String)
+function read_datafile(filename::String, ::Union{TreeMesh, StructuredMesh,
+                                                 UnstructuredMesh2D, P4estMesh, T8codeMesh})
   # Open file for reading
   h5open(filename, "r") do file
     # Extract basic information
@@ -79,3 +80,104 @@ function read_datafile(filename::String)
     return labels, data, n_elements, n_nodes, element_variables, node_variables, time
   end
 end
+
+# Read in data file and return all relevant information for a DGMulti solution
+function read_datafile(filename::String, ::DGMultiMesh)
+  # Open file for reading
+  h5open(filename, "r") do file
+    # Extract basic information
+    ndims_ = read(attributes(file)["ndims"])
+
+    etype_str = read(attributes(file)["element_type"])
+    etype = Trixi.get_element_type_from_string(etype_str)()
+
+    if etype isa Trixi.Wedge && haskey(attributes(file), "polydeg_tri")
+        polydeg = tuple(read(attributes(file)["polydeg_tri"]),
+                        read(attributes(file)["polydeg_line"]))
+    else
+        polydeg = read(attributes(file)["polydeg"])
+    end
+
+    n_elements = read(attributes(file)["n_elements"])
+    dof_per_elem = read(attributes(file)["dof_per_elem"])
+    n_variables = read(attributes(file)["n_vars"])
+    time = read(attributes(file)["time"])
+
+    # Extract labels for legend
+    labels = Array{String}(undef, 1, n_variables)
+    for v = 1:n_variables
+      labels[1, v] = read(attributes(file["variables_$v"])["name"])
+    end
+
+    # Extract data arrays
+    n_nodes = maximum(polydeg) + 1
+
+    data = Array{Float64}(undef, dof_per_elem, n_elements, n_variables)
+    for v = 1:n_variables
+      vardata = read(file["variables_$v"])
+      @views data[:, :, v] .= vardata
+    end
+
+    # Extract element variable arrays
+    element_variables = Dict{String, Union{Vector{Float64}, Vector{Int}}}()
+    index = 1
+    while haskey(file, "element_variables_$index")
+      varname = read(attributes(file["element_variables_$index"])["name"])
+      element_variables[varname] = read(file["element_variables_$index"])
+      index +=1
+    end
+
+    # Extract node variable arrays
+    node_variables = Dict{String, Union{Array{Float64}, Array{Int}}}()
+    index = 1
+    while haskey(file, "node_variables_$index")
+      varname = read(attributes(file["node_variables_$index"])["name"])
+      nodedata = read(file["node_variables_$index"])
+      if ndims_ == 2
+        node_variables[varname] = Array{Float64}(undef, n_nodes, n_nodes, n_elements)
+        @views node_variables[varname][:, :, :] .= nodedata
+      else
+        error("Unsupported number of spatial dimensions: ", ndims_)
+      end
+      index +=1
+    end
+
+    return labels, data, n_elements, n_nodes, element_variables, node_variables, time, etype, polydeg
+  end
+end
+
+# Reconstruct a `RefElemData` from information in `mesh_file`
+function load_basis(mesh_file, ::DGMultiMesh)
+  etype_str = h5open(mesh_file, "r") do file
+      return read(attributes(file)["element_type"])
+  end
+
+  etype = Trixi.get_element_type_from_string(etype_str)()
+  polydeg = h5open(mesh_file, "r") do file
+      if etype isa Trixi.Wedge && haskey(attributes(file), "polydeg_tri")
+          return tuple(read(attributes(file)["polydeg_tri"]),
+                      read(attributes(file)["polydeg_line"]))
+      else
+          return read(attributes(file)["polydeg"])
+      end
+  end
+
+  if etype isa Trixi.Wedge && polydeg isa NTuple{2}
+    factor_a = Trixi.RefElemData(Trixi.StartUpDG.Tri(), Trixi.Polynomial(), polydeg[1])
+    factor_b = Trixi.RefElemData(Trixi.StartUpDG.Line(), Trixi.Polynomial(), polydeg[2])
+
+    tensor = Trixi.TensorProductWedge(factor_a, factor_b)
+    rd = Trixi.RefElemData(etype, tensor)
+  else
+    rd = Trixi.RefElemData(etype, Trixi.Polynomial(), polydeg)
+  end
+
+  return rd
+end
+
+# For other mesh types, return nothing
+function load_basis(mesh_file, ::Union{TreeMesh, StructuredMesh,
+                                       UnstructuredMesh2D, P4estMesh, 
+                                       T8codeMesh})
+  return nothing
+end