Adding T8codeMesh support.

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
+             gauss_lobatto_nodes_weights, TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, vtk_save, paraview_collection
 
 # Include all top-level submodule files

src/convert.jl

@@ -314,7 +314,7 @@ function assert_cells_elements(n_elements, mesh::UnstructuredMesh2D, filename, m
 end
 
 
-function assert_cells_elements(n_elements, mesh::P4estMesh, filename, meshfile)
+function assert_cells_elements(n_elements, mesh::Union{P4estMesh, T8codeMesh}, 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 +336,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})
+                                        mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh})
   if nvisnodes === nothing || nvisnodes == 0
     return n_nodes
   else
@@ -356,7 +356,7 @@ function get_default_nvisnodes_mesh(nvisnodes, mesh::TreeMesh)
 end
 
 function get_default_nvisnodes_mesh(nvisnodes,
-                                    mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh})
+                                    mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh})
   if nvisnodes === nothing
     # for curved meshes, we need to get at least the vertices
     return 2
@@ -434,6 +434,34 @@ function add_celldata!(vtk_celldata, mesh::P4estMesh, verbose)
   return vtk_celldata
 end
 
+function add_celldata!(vtk_celldata, mesh::T8codeMesh, verbose)
+  # Create temporary storage for the tree_ids and levels.
+  tree_ids = zeros( Trixi.ncells(mesh) )
+
+  elem_counter = 1
+  num_local_trees = Trixi.t8_forest_get_num_local_trees(mesh.forest)
+  for itree in 1:num_local_trees
+      num_elements_in_tree = Trixi.t8_forest_get_tree_num_elements(mesh.forest, itree-1)
+      for ielement in 1:num_elements_in_tree
+          tree_ids[elem_counter] = itree
+          elem_counter += 1
+      end
+  end
+
+  levels = Trixi.trixi_t8_get_local_element_levels(mesh.forest)
+
+  @timeit "add data to VTK file" begin
+    # Add tree/element data to celldata VTK file
+    verbose && println("| | tree_ids...")
+    @timeit "tree_ids" vtk_celldata["tree_ids"] = tree_ids
+    verbose && println("| | element_ids...")
+    @timeit "element_ids" vtk_celldata["element_ids"] = collect(1:Trixi.ncells(mesh))
+    verbose && println("| | levels...")
+    @timeit "levels" vtk_celldata["levels"] = levels
+  end
+
+  return vtk_celldata
+end
 
 function expand_filename_patterns(patterns)
   filenames = String[]

src/interpolate.jl

@@ -11,7 +11,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},
+                          mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh},
                           n_visnodes, verbose)
   # Calculate equidistant output nodes
   nodes_out = collect(range(-1, 1, length=n_visnodes))

src/vtktools.jl

@@ -137,7 +137,7 @@ end
 # (StructuredMesh/UnstructuredMesh2D/P4estMesh version).
 # Routine is agnostic with respect to reinterpolation.
 function build_vtk_grids(::Val{:vtu},
-                         mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh},
+                         mesh::Union{StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh},
                          nodes, n_visnodes, verbose, output_directory, is_datafile, filename,
                          reinterpolate::Union{Val{true}, Val{false}})
 
@@ -230,8 +230,7 @@ function calc_node_coordinates(mesh::UnstructuredMesh2D, nodes, n_visnodes)
   return node_coordinates
 end
 
-
-function calc_node_coordinates(mesh::P4estMesh, nodes, n_visnodes)
+function calc_node_coordinates(mesh::Union{P4estMesh,T8codeMesh}, nodes, n_visnodes)
   # Extract number of spatial dimensions
   ndims_ = ndims(mesh)
 

test/test_2d.jl

@@ -394,6 +394,92 @@ end
         compare_point_data(out_file, ref_file)
       end
     end
+
+    @testset "T8codeMesh" begin
+      isdir(outdir) && rm(outdir, recursive=true)
+      run_trixi(joinpath(examples_dir(), "t8code_2d_dgsem", "elixir_mhd_rotor.jl"), maxiters=5)
+
+      @timed_testset "mesh data" begin
+        # create the output file to be tested
+        @test_nowarn trixi2vtk(joinpath(outdir, "mesh_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "mesh_" * LEADING_ZEROS * "000005_celldata.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-t8codemesh"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/t8codemesh/dgsem_rotor_amr_mesh_05.vtu"
+        ref_file = get_test_reference_file("dgsem_rotor_amr_mesh_05.vtu", remote_filename)
+        compare_cell_data(out_file, ref_file)
+      end
+
+      @timed_testset "solution celldata" begin
+        # create the output file to be tested
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000005_celldata.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-t8codemesh"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/t8codemesh/dgsem_rotor_amr_celldata_05.vtu"
+        ref_file = get_test_reference_file("dgsem_rotor_amr_celldata_05.vtu", remote_filename)
+        compare_cell_data(out_file, ref_file)
+      end
+
+      @timed_testset "reinterpolate with nonuniform data" begin
+        # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-t8codemesh-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/t8codemesh/dgsem_rotor_amr_reinterp_05.vtu"
+        ref_file = get_test_reference_file("dgsem_rotor_amr_reinterp_05.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+
+      @timed_testset "do not reinterpolate with nonuniform data" begin
+        # Create and test output without reinterpolation on LGL nodes
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-t8codemesh-no-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/t8codemesh/dgsem_rotor_amr_no_reinterp_05.vtu"
+        ref_file = get_test_reference_file("dgsem_rotor_amr_no_reinterp_05.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+
+      @timed_testset "do not reinterpolate with uniform data" begin
+        # Create and test output without reinterpolation on uniform nodes
+        # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-t8codemesh-no-reinterp-uniform"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/t8codemesh/dgsem_rotor_amr_no_reinterp_uniform_05.vtu"
+        ref_file = get_test_reference_file("dgsem_rotor_amr_no_reinterp_uniform_05.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+    end
   end
 
   if VERSION >= v"1.8"

test/test_3d.jl

@@ -276,6 +276,92 @@ end
         compare_point_data(out_file, ref_file)
       end
     end
+
+    @testset "T8codeMesh" begin
+      isdir(outdir) && rm(outdir, recursive=true)
+      run_trixi(joinpath(examples_dir(), "t8code_3d_dgsem", "elixir_advection_cubed_sphere.jl"), maxiters=2)
+
+      @timed_testset "mesh data" begin
+        # create the output file to be tested
+        @test_nowarn trixi2vtk(joinpath(outdir, "mesh.h5"), output_directory=outdir)
+        outfilename = "mesh_celldata.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "3d-t8codemesh"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "3d/t8codemesh/dgsem_adv_sphere_mesh_02.vtu"
+        ref_file = get_test_reference_file("dgsem_adv_sphere_mesh_02.vtu", remote_filename)
+        compare_cell_data(out_file, ref_file)
+      end
+
+      @timed_testset "solution celldata" begin
+        # create the output file to be tested
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000002_celldata.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "3d-t8codemesh"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "3d/t8codemesh/dgsem_adv_sphere_celldata_02.vtu"
+        ref_file = get_test_reference_file("dgsem_adv_sphere_celldata_02.vtu", remote_filename)
+        compare_cell_data(out_file, ref_file)
+      end
+
+      @timed_testset "reinterpolate with nonuniform data" begin
+        # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "3d-t8codemesh-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "3d/t8codemesh/dgsem_adv_sphere_reinterp_02.vtu"
+        ref_file = get_test_reference_file("dgsem_adv_sphere_reinterp_02.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+
+      @timed_testset "do not reinterpolate with nonuniform data" begin
+        # Create and test output without reinterpolation on LGL nodes
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "3d-t8codemesh-no-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "3d/t8codemesh/dgsem_adv_sphere_no_reinterp_02.vtu"
+        ref_file = get_test_reference_file("dgsem_adv_sphere_no_reinterp_02.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+
+      @timed_testset "do not reinterpolate with uniform data" begin
+        # Create and test output without reinterpolation on uniform nodes
+        # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "3d-t8codemesh-no-reinterp-uniform"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "3d/t8codemesh/dgsem_adv_sphere_no_reinterp_uniform_02.vtu"
+        ref_file = get_test_reference_file("dgsem_adv_sphere_no_reinterp_uniform_02.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+    end
   end
 end