diff --git a/src/Trixi2Vtk.jl b/src/Trixi2Vtk.jl
index 1df205d..52963c8 100644
--- a/src/Trixi2Vtk.jl
+++ b/src/Trixi2Vtk.jl
@@ -7,7 +7,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..ec8e980 100644
--- a/src/convert.jl
+++ b/src/convert.jl
@@ -127,8 +127,15 @@ 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)
+
+      # Check compatibility of the mesh type and the output format
+      if mesh isa Trixi.DGMultiMesh
+        @timeit "read data" (labels, data, n_elements, n_nodes,
+                             element_variables, node_variables, time) = read_datafile_dgmulti(filename)
+      else
+        @timeit "read data" (labels, data, n_elements, n_nodes,
+                             element_variables, node_variables, time) = read_datafile(filename)
+      end
 
       assert_cells_elements(n_elements, mesh, filename, meshfile)
 
@@ -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
@@ -314,7 +322,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 +344,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 +364,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 +471,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 36e8572..ff8b77c 100644
--- a/src/interpolate.jl
+++ b/src/interpolate.jl
@@ -47,6 +47,32 @@ 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 (StructuredMesh or UnstructuredMesh2D version)
+function interpolate_data(::Val{:vtu}, input_data,
+                          mesh::DGMultiMesh,
+                          n_visnodes, verbose)
+  rd = mesh.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
+
+  interpolator = Trixi.StartUpDG.vandermonde(rd.element_type, order, Trixi.StartUpDG.equi_nodes(rd.element_type, order)...) / rd.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 997cdf5..5c90fff 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -94,3 +94,68 @@ 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
+function read_datafile_dgmulti(filename::String)
+  # 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
+        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
+  end
+end
diff --git a/src/vtktools.jl b/src/vtktools.jl
index ade8666..9075d31 100644
--- a/src/vtktools.jl
+++ b/src/vtktools.jl
@@ -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, 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, mesh.rd, 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, mesh.rd, 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)
@@ -241,6 +278,12 @@ function calc_node_coordinates(mesh::Union{P4estMesh,T8codeMesh}, nodes, n_visno
   return Trixi.calc_node_coordinates!(node_coordinates, mesh, nodes)
 end
 
+# function calc_node_coordinates(mesh::DGMultiMesh, varags...)
+#   # Extract number of spatial dimensions
+#   ndims_ = ndims(mesh)
+# 
+#   return mesh.md.x, mesh.md.y
+# end
 
 # Calculation of the node coordinates for `TreeMesh` in 2D
 function calc_node_coordinates!(node_coordinates, nodes, mesh::TreeMesh{2})
@@ -603,7 +646,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)
@@ -660,3 +702,32 @@ 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 `StructuredMesh` (2D version)
+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