diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..feca861
Binary files /dev/null and b/.DS_Store differ
diff --git a/src/Trixi2Vtk.jl b/src/Trixi2Vtk.jl
index e636ed8..d73319a 100644
--- a/src/Trixi2Vtk.jl
+++ b/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
diff --git a/src/convert.jl b/src/convert.jl
index 5da6159..6a08a97 100644
--- a/src/convert.jl
+++ b/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,8 +172,8 @@ 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
@@ -173,7 +181,8 @@ function trixi2vtk(filename::AbstractString...;
       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[]
 
diff --git a/src/interpolate.jl b/src/interpolate.jl
index c5af52d..e48d4f5 100644
--- a/src/interpolate.jl
+++ b/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},
diff --git a/src/io.jl b/src/io.jl
index 8748d59..b1a94d7 100644
--- a/src/io.jl
+++ b/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
\ No newline at end of file
diff --git a/src/vtktools.jl b/src/vtktools.jl
index 270d9ce..015b6a6 100644
--- a/src/vtktools.jl
+++ b/src/vtktools.jl
@@ -1,6 +1,6 @@
 # Create and return VTK grids on enriched uniform nodes
 # that are ready to be filled with reinterpolated data (vtu version)
-function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, nodes, n_visnodes, verbose,
+function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh, nodes, basis, n_visnodes, verbose,
                          output_directory, is_datafile, filename, reinterpolate::Val{true})
   coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh)
 
@@ -47,7 +47,7 @@ end
 # Create and return VTK grids with nodes on which the raw was created
 # ready to be filled with raw data (vtu version)
 function build_vtk_grids(::Val{:vtu}, mesh::TreeMesh,
-                         nodes, n_visnodes, verbose,
+                         nodes, basis, n_visnodes, verbose,
                          output_directory, is_datafile, filename, reinterpolate::Val{false})
 
   @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, nodes, n_visnodes)
@@ -88,7 +88,7 @@ end
 
 # Create and return VTK grids on enriched uniform nodes
 # that are ready to be filled with reinterpolated data (vti version)
-function build_vtk_grids(::Val{:vti}, mesh::TreeMesh, nodes, n_visnodes, verbose,
+function build_vtk_grids(::Val{:vti}, mesh::TreeMesh, nodes, basis, n_visnodes, verbose,
                          output_directory, is_datafile, filename, reinterpolate::Val{true})
 
   coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh)
@@ -138,7 +138,7 @@ end
 # Routine is agnostic with respect to reinterpolation.
 function build_vtk_grids(::Val{:vtu},
                          mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh},
-                         nodes, n_visnodes, verbose, output_directory, is_datafile, filename,
+                         nodes, basis, n_visnodes, verbose, output_directory, is_datafile, filename,
                          reinterpolate::Union{Val{true}, Val{false}})
 
   @timeit "prepare coordinate information" node_coordinates = calc_node_coordinates(mesh, nodes, n_visnodes)
@@ -176,6 +176,43 @@ function build_vtk_grids(::Val{:vtu},
   return vtk_nodedata, vtk_celldata
 end
 
+function build_vtk_grids(::Val{:vtu},
+                         mesh::DGMultiMesh,
+                         nodes, basis, n_visnodes, verbose, output_directory, is_datafile, filename,
+                         reinterpolate::Union{Val{true}, Val{false}})
+
+  # 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(mesh.md, basis, reinterpolate)
+    end
+  end
+
+  # Prepare VTK points and cells for celldata file.
+  @timeit "prepare VTK cells (cell data)" begin
+    vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells(mesh.md, basis, reinterpolate)
+  end
+
+  # 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
+    @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, vtk_points...,
+                                                                 vtk_cells)
+  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
 
 function calc_node_coordinates(mesh::TreeMesh, nodes, n_visnodes)
   coordinates, levels, _, _ = extract_mesh_information(mesh)
@@ -256,7 +293,6 @@ function calc_node_coordinates(mesh::T8codeMesh, nodes, n_visnodes)
   return Trixi.calc_node_coordinates!(node_coordinates, mesh, nodes)
 end
 
-
 # Calculation of the node coordinates for `TreeMesh` in 2D
 function calc_node_coordinates!(node_coordinates, nodes, mesh::TreeMesh{2})
   _, levels, _, length_level_0 = extract_mesh_information(mesh)
@@ -620,7 +656,6 @@ function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,4})
   return vtk_points, vtk_cells
 end
 
-
 # Convert coordinates and level information to a list of points and VTK cells for `StructuredMesh` (3D version)
 function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,5})
   n_elements = size(node_coordinates, 5)
@@ -677,3 +712,34 @@ function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,5})
 
   return vtk_points, vtk_cells
 end
+
+# Convert coordinates and level information to a list of points and VTK cells for
+# DGMultiMesh
+function calc_vtk_points_cells(md::Trixi.MeshData, rd::Trixi.RefElemData, 
+                               reinterpolate = Val(true))
+  # Compute the permutation between the StartUpDG order of points and vtk.
+  perm = Trixi.StartUpDG.SUD_to_vtk_order(rd)
+
+  if length(rd.N) > 1
+      @assert length(Set(rd.N)) == 1 "`order` must have equal elements."
+      order = first(rd.N)
+  else
+      order = rd.N
+  end
+
+  # The number of points per element.
+  num_lagrange_points = length(perm)
+  vtk_cell_type = Trixi.StartUpDG.type_to_vtk(rd.element_type)
+
+  # Construction of VTK mesh cell list.
+  vtk_cells = [MeshCell(vtk_cell_type, perm .+ ((i-1) * num_lagrange_points)) for i in 1:md.num_elements]
+
+  if reinterpolate == Val(true)
+      interpolate = Trixi.StartUpDG.vandermonde(rd.element_type, order, Trixi.StartUpDG.equi_nodes(rd.element_type, order)...) / rd.VDM
+      vtk_points = map(x -> vec(interpolate * x), md.xyz)
+  else # don't interpolate
+      vtk_points = vec.(md.xyz) 
+  end   
+  
+  return vtk_points, vtk_cells
+end