Add capability to not perform reinterpolation

.gitignore

@@ -12,6 +12,7 @@
 **/Manifest.toml
 out/
 artifacts/
+reference_files/
 docs/build
 public/
 coverage/

src/convert.jl

@@ -1,7 +1,8 @@
 """
     trixi2vtk(filename::AbstractString...;
               format=:vtu, verbose=false, hide_progress=false, pvd=nothing,
-              output_directory=".", nvisnodes=nothing)
+              output_directory=".", nvisnodes=nothing, save_celldata=true,
+              reinterpolate=true, data_is_uniform=false)
 
 Convert Trixi-generated output files to VTK files (VTU or VTI).
 
@@ -21,6 +22,13 @@ Convert Trixi-generated output files to VTK files (VTU or VTI).
                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.
                    (default: `true`)
+- `reinterpolate`: Boolean value to determine if data should be reinterpolated
+                   onto uniform points. When `false` the raw data at the compute nodes
+                   is copied into the appropriate format.
+                   (default: `true`)
+- `data_is_uniform`: Boolean to indicate if the data to be converted is from a finite difference
+                     method on a uniform grid of points.
+                     (default: `false`)
 
 # Examples
 ```julia
@@ -30,7 +38,8 @@ julia> trixi2vtk("out/solution_000*.h5")
 """
 function trixi2vtk(filename::AbstractString...;
                    format=:vtu, verbose=false, hide_progress=false, pvd=nothing,
-                   output_directory=".", nvisnodes=nothing, save_celldata=true)
+                   output_directory=".", nvisnodes=nothing, save_celldata=true,
+                   reinterpolate=true, data_is_uniform=false)
   # Reset timer
   reset_timer!()
 
@@ -103,6 +112,11 @@ function trixi2vtk(filename::AbstractString...;
     verbose && println("| Reading mesh file...")
     @timeit "read mesh" mesh = Trixi.load_mesh_serial(meshfile; n_cells_max=0, RealT=Float64)
 
+    # Check compatibility of the mesh type and the output format
+    if format == :vti && !(mesh isa Trixi.TreeMesh)
+      error("VTI format only available for 2D TreeMesh")
+    end
+
     # Read data only if it is a data file
     if is_datafile
       verbose && println("| Reading data file...")
@@ -113,24 +127,53 @@ function trixi2vtk(filename::AbstractString...;
 
       # Determine resolution for data interpolation
       n_visnodes = get_default_nvisnodes_solution(nvisnodes, n_nodes, mesh)
+
+      # If a user requests that no reinterpolation is done automatically set
+      # `n_visnodes` to be the same as the number of nodes in the raw data.
+      if !reinterpolate
+        n_visnodes = n_nodes
+      end
     else
       # If file is a mesh file, do not interpolate data as detailed
       n_visnodes = get_default_nvisnodes_mesh(nvisnodes, mesh)
     end
 
+    # Check if the raw data is uniform (finite difference) or not (dg)
+    # and create the corresponding node set for reinterpolation / copying.
+    if (reinterpolate & !data_is_uniform) | (!reinterpolate & data_is_uniform)
+      # (1) Default settings; presumably the most common
+      # (2) Finite difference data
+      node_set = range(-1, 1, length=n_visnodes)
+    elseif !reinterpolate & !data_is_uniform
+      # raw data is on a set of LGL nodes
+      node_set, _ = gauss_lobatto_nodes_weights(n_visnodes)
+    else # reinterpolate & data_is_uniform
+      error("uniform data should not be reinterpolated! Set reinterpolate=false and try again.")
+    end
+
     # Create output directory if it does not exist
     mkpath(output_directory)
 
     # Build VTK grids
-    vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, n_visnodes, verbose,
-                                                 output_directory, is_datafile, filename)
+    vtk_nodedata, vtk_celldata = build_vtk_grids(Val(format), mesh, node_set, n_visnodes, verbose,
+                                                 output_directory, is_datafile, filename, Val(reinterpolate))
 
     # Interpolate data
     if is_datafile
       verbose && println("| Interpolating data...")
-      @timeit "interpolate data" interpolated_data = interpolate_data(Val(format),
-                                                                      data, mesh,
-                                                                      n_visnodes, verbose)
+      if reinterpolate
+        @timeit "interpolate data" interpolated_data = interpolate_data(Val(format),
+                                                                        data, mesh,
+                                                                        n_visnodes, verbose)
+      else # Copy the raw solution data; only works for `vtu` format
+        # Extract data shape information
+        ndims_ = ndims(data) - 2
+        n_variables = length(labels)
+        # Save raw data as one 1D array for each variable
+        @timeit "interpolate data" interpolated_data = reshape(data,
+                                                               n_visnodes^ndims_ * n_elements,
+                                                               n_variables)
+      end
     end
 
     # Add data to file
@@ -272,19 +315,10 @@ end
 
 
 # default number of visualization nodes if only the mesh should be visualized
-function get_default_nvisnodes_mesh(nvisnodes, mesh::TreeMesh)
-  if nvisnodes === nothing
-    # for a Cartesian mesh, we do not need to interpolate
-    return 1
-  else
-    return nvisnodes
-  end
-end
-
 function get_default_nvisnodes_mesh(nvisnodes,
-                                    mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh})
+                                    mesh::Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh})
   if nvisnodes === nothing
-    # for curved meshes, we need to get at least the vertices
+    # we need to get at least the vertices
     return 2
   else
     return nvisnodes

src/interpolate.jl

@@ -21,10 +21,10 @@ end
 
 
 # Interpolate data from input format to desired output format (vti version)
-function interpolate_data(::Val{:vti}, input_data, mesh, n_visnodes, verbose)
+function interpolate_data(::Val{:vti}, input_data, mesh::TreeMesh, 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]²
+  # Normalize element coordinates: move center to origin and domain size to [-1, 1]²
   normalized_coordinates = similar(coordinates)
   for element_id in axes(coordinates, 2)
     @views normalized_coordinates[:, element_id] .= (
@@ -77,7 +77,6 @@ function unstructured2structured(unstructured_data::AbstractArray{Float64},
 
   # Create output data structure
   structured = Array{Float64}(undef, resolution, resolution, n_variables)
-
   # For each variable, interpolate element data and store to global data structure
   for v in 1:n_variables
     # Reshape data array for use in interpolate_nodes function
@@ -135,7 +134,7 @@ end
 # Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1])
 function coordinate2index(coordinate, resolution::Integer)
   # Calculate 1D normalized coordinates
-  dx = 2/resolution
+  dx = 2 / resolution
   mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))
 
   # Find index
@@ -211,6 +210,7 @@ function interpolate_nodes(data_in::AbstractArray{T, 3},
   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)