Feature t8code: Extending to 3D

.github/workflows/ci.yml

@@ -70,6 +70,7 @@ jobs:
           - p4est_part1
           - p4est_part2
           - t8code_part1
+          - t8code_part2
           - unstructured_dgmulti
           - parabolic
           - paper_self_gravitating_gas_dynamics

examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl

@@ -93,12 +93,10 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 coordinates_min = (-5.0, -5.0)
 coordinates_max = (5.0, 5.0)
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (1, 1)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 1)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl

@@ -41,9 +41,9 @@ isfile(mesh_file) ||
 # we can create a t8code mesh.
 conn = Trixi.read_inp_p4est(mesh_file, Val(2))
 
-mesh = T8codeMesh{2}(conn, polydeg = 3,
-                     mapping = mapping_flag,
-                     initial_refinement_level = 1)
+mesh = T8codeMesh(conn, polydeg = 3,
+                  mapping = mapping_flag,
+                  initial_refinement_level = 1)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)

examples/t8code_2d_dgsem/elixir_advection_basic.jl

@@ -16,12 +16,10 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
 coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (8, 8)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 1)
 
 # A semidiscretization collects data structures and functions for the spatial discretization

examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl

@@ -20,31 +20,28 @@ f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s))
 faces = (f1, f2, f3, f4)
 mapping = Trixi.transfinite_mapping(faces)
 
-# Create P4estMesh with 3 x 2 trees and 6 x 4 elements,
+# Create T8codeMesh with 3 x 2 trees and 6 x 4 elements,
 # approximate the geometry with a smaller polydeg for testing.
 trees_per_dimension = (3, 2)
 mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
                   mapping = mapping,
                   initial_refinement_level = 1)
 
-function adapt_callback(forest,
-                        forest_from,
-                        which_tree,
-                        lelement_id,
-                        ts,
-                        is_family,
-                        num_elements,
-                        elements_ptr)::Cint
-    vertex = Vector{Cdouble}(undef, 3)
-
-    elements = unsafe_wrap(Array, elements_ptr, num_elements)
-
-    Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex))
+# Note: This is actually a `p4est_quadrant_t` which is much bigger than the
+# following struct. But we only need the first three fields for our purpose.
+struct t8_dquad_t
+    x::Int32
+    y::Int32
+    level::Int8
+    # [...] # See `p4est.h` in `p4est` for more info.
+end
 
-    level = Trixi.t8_element_level(ts, elements[1])
+# Refine quadrants of each tree at lower left edge to level 4.
+function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family,
+                        user_data)
+    el = unsafe_load(Ptr{t8_dquad_t}(elements[1]))
 
-    # TODO: Make this condition more general.
-    if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 4
+    if el.x == 0 && el.y == 0 && el.level < 4
         # return true (refine)
         return 1
     else
@@ -53,26 +50,7 @@ function adapt_callback(forest,
     end
 end
 
-Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0);
-
-# Init new forest.
-new_forest_ref = Ref{Trixi.t8_forest_t}()
-Trixi.t8_forest_init(new_forest_ref);
-new_forest = new_forest_ref[]
-
-# Check out `examples/t8_step4_partition_balance_ghost.jl` in
-# https://github.com/DLR-AMR/T8code.jl for detailed explanations.
-let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0
-    Trixi.t8_forest_set_user_data(new_forest, C_NULL)
-    Trixi.t8_forest_set_adapt(new_forest, mesh.forest,
-                              Trixi.@t8_adapt_callback(adapt_callback), recursive)
-    Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition)
-    Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening)
-    Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES)
-    Trixi.t8_forest_commit(new_forest)
-end
-
-mesh.forest = new_forest
+adapt!(mesh, adapt_callback)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,

examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl

@@ -38,9 +38,9 @@ isfile(mesh_file) ||
 # we can create a t8code mesh.
 conn = Trixi.read_inp_p4est(mesh_file, Val(2))
 
-mesh = T8codeMesh{2}(conn, polydeg = 3,
-                     mapping = mapping_flag,
-                     initial_refinement_level = 2)
+mesh = T8codeMesh(conn, polydeg = 3,
+                  mapping = mapping_flag,
+                  initial_refinement_level = 2)
 
 # A semidiscretization collects data structures and functions for the spatial discretization.
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

examples/t8code_2d_dgsem/elixir_euler_free_stream.jl

@@ -40,25 +40,17 @@ isfile(mesh_file) ||
 # we can create a t8code mesh.
 conn = Trixi.read_inp_p4est(mesh_file, Val(2))
 
-mesh = T8codeMesh{2}(conn, polydeg = 3,
-                     mapping = mapping,
-                     initial_refinement_level = 1)
-
-function adapt_callback(forest,
-                        forest_from,
-                        which_tree,
-                        lelement_id,
-                        ts,
-                        is_family,
-                        num_elements,
-                        elements_ptr)::Cint
-    vertex = Vector{Cdouble}(undef, 3)
+mesh = T8codeMesh(conn, polydeg = 3,
+                  mapping = mapping,
+                  initial_refinement_level = 1)
 
-    elements = unsafe_wrap(Array, elements_ptr, num_elements)
+function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family,
+                        user_data)
+    vertex = Vector{Cdouble}(undef, 3)
 
-    Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex))
+    Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, vertex)
 
-    level = Trixi.t8_element_level(ts, elements[1])
+    level = Trixi.t8_element_level(eclass_scheme, elements[1])
 
     # TODO: Make this condition more general.
     if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 3
@@ -70,26 +62,7 @@ function adapt_callback(forest,
     end
 end
 
-Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0);
-
-# Init new forest.
-new_forest_ref = Ref{Trixi.t8_forest_t}()
-Trixi.t8_forest_init(new_forest_ref);
-new_forest = new_forest_ref[]
-
-# Check out `examples/t8_step4_partition_balance_ghost.jl` in
-# https://github.com/DLR-AMR/T8code.jl for detailed explanations.
-let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0
-    Trixi.t8_forest_set_user_data(new_forest, C_NULL)
-    Trixi.t8_forest_set_adapt(new_forest, mesh.forest,
-                              Trixi.@t8_adapt_callback(adapt_callback), recursive)
-    Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition)
-    Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening)
-    Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES)
-    Trixi.t8_forest_commit(new_forest)
-end
-
-mesh.forest = new_forest
+adapt!(mesh, adapt_callback)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)))

examples/t8code_2d_dgsem/elixir_euler_sedov.jl

@@ -58,12 +58,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (4, 4)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 4,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 2, periodicity = true)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl

@@ -29,12 +29,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (4, 4)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 4,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 2, periodicity = true)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl

@@ -40,25 +40,17 @@ isfile(mesh_file) ||
 # we can create a t8code mesh.
 conn = Trixi.read_inp_p4est(mesh_file, Val(2))
 
-mesh = T8codeMesh{2}(conn, polydeg = 3,
-                     mapping = mapping_flag,
-                     initial_refinement_level = 1)
-
-function adapt_callback(forest,
-                        forest_from,
-                        which_tree,
-                        lelement_id,
-                        ts,
-                        is_family,
-                        num_elements,
-                        elements_ptr)::Cint
-    vertex = Vector{Cdouble}(undef, 3)
+mesh = T8codeMesh(conn, polydeg = 3,
+                  mapping = mapping_flag,
+                  initial_refinement_level = 1)
 
-    elements = unsafe_wrap(Array, elements_ptr, num_elements)
+function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family,
+                        user_data)
+    vertex = Vector{Cdouble}(undef, 3)
 
-    Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex))
+    Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, pointer(vertex))
 
-    level = Trixi.t8_element_level(ts, elements[1])
+    level = Trixi.t8_element_level(eclass_scheme, elements[1])
 
     # TODO: Make this condition more general.
     if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 2
@@ -70,26 +62,7 @@ function adapt_callback(forest,
     end
 end
 
-@assert(Trixi.t8_forest_is_committed(mesh.forest)!=0);
-
-# Init new forest.
-new_forest_ref = Ref{Trixi.t8_forest_t}()
-Trixi.t8_forest_init(new_forest_ref);
-new_forest = new_forest_ref[]
-
-# Check out `examples/t8_step4_partition_balance_ghost.jl` in
-# https://github.com/DLR-AMR/T8code.jl for detailed explanations.
-let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0
-    Trixi.t8_forest_set_user_data(new_forest, C_NULL)
-    Trixi.t8_forest_set_adapt(new_forest, mesh.forest,
-                              Trixi.@t8_adapt_callback(adapt_callback), recursive)
-    Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition)
-    Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening)
-    Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES)
-    Trixi.t8_forest_commit(new_forest)
-end
-
-mesh.forest = new_forest
+adapt!(mesh, adapt_callback)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms = source_terms,

examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl

@@ -16,10 +16,8 @@ coordinates_max = (2.0, 2.0)
 
 trees_per_dimension = (1, 1)
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 mesh = T8codeMesh(trees_per_dimension, polydeg = 1,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 2)
 
 semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition,

examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl

@@ -11,18 +11,17 @@ initial_condition = initial_condition_convergence_test
 
 # Get the DG approximation space
 volume_flux = (flux_central, flux_nonconservative_powell)
+
 solver = DGSEM(polydeg = 4, surface_flux = (flux_hlle, flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (sqrt(2.0), sqrt(2.0))
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (8, 8)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 0, periodicity = true)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_dgsem/elixir_mhd_rotor.jl

@@ -78,9 +78,9 @@ isfile(mesh_file) ||
 # we can create a t8code mesh.
 conn = Trixi.read_inp_p4est(mesh_file, Val(2))
 
-mesh = T8codeMesh{2}(conn, polydeg = 4,
-                     mapping = mapping_twist,
-                     initial_refinement_level = 1)
+mesh = T8codeMesh(conn, polydeg = 4,
+                  mapping = mapping_twist,
+                  initial_refinement_level = 1)
 
 boundary_condition = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = Dict(:all => boundary_condition)

examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl

@@ -22,12 +22,10 @@ solver = DGSEM(polydeg = 3,
 coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
 coordinates_max = (sqrt(2.0), sqrt(2.0))  # maximum coordinates (max(x), max(y))
 
-mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
-
 trees_per_dimension = (8, 8)
 
 mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
-                  mapping = mapping,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                   initial_refinement_level = 1)
 
 # A semidiscretization collects data structures and functions for the spatial discretization