Feature: Read geometry information from t8code (#1900)

* Wrap from_abaqus routines. * Implement geometry data transfer from t8code to Trixi. * Updated examples. * Fixed typos. * Applied formatter. * cubed sphere test case, copied from p4est * add baroclinic instability (copy of p4est) * add cubed sphere constructor * Fix indentation. * Fix indentation. * Switching off formatter in two files. * Upgrading T8code.jl. * Fixed examples. * stress different meaning of first argument it refers to level of refinement in lat lon direction, not number of tree as in the p4est version * use lat lon levels * add t8code cubed sphere tests * Remove TODO comments * Relaxing T8code.jl version requirement. * Restricted t8code version requirement. * Removed cubed spherical shell related code. * Removed cubed spherical shell tests. * Increasing code coverage. * Increasing code coverage. * Further increasing code coverage. * Applied formatter. * Update src/meshes/t8code_mesh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update test/test_t8code_2d.jl Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/meshes/t8code_mesh.jl Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> * Fixing minor stuff. * Update src/meshes/t8code_mesh.jl Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> * Appyling comments. * Applied formatter. * Reverting changes to original state. --------- Co-authored-by: Johannes Markert <johannes.markert@dlr.de> Co-authored-by: Benedict Geihe <bgeihe@uni-koeln.de> Co-authored-by: Benedict <135045760+benegee@users.noreply.github.com> Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
trixi-framework · May 10, 2024 · d1c20c6 · d1c20c6
1 parent 8a9fc7b
commit d1c20c6
Show file tree

Hide file tree

Showing 17 changed files with 411 additions and 246 deletions.
diff --git a/Project.toml b/Project.toml
@@ -92,7 +92,7 @@ StaticArrays = "1.5"
 StrideArrays = "0.1.26"
 StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.4.3, 0.5"
+T8code = "0.5"
 TimerOutputs = "0.5.7"
 Triangulate = "2.2"
 TriplotBase = "0.1"

diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
@@ -33,13 +33,8 @@ mapping_flag = Trixi.transfinite_mapping(faces)
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                            joinpath(@__DIR__, "square_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(2))
-
-mesh = T8codeMesh(conn, polydeg = 3,
-                  mapping = mapping_flag,
+mesh = T8codeMesh(mesh_file, 2;
+                  mapping = mapping_flag, polydeg = 3,
                   initial_refinement_level = 1)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

diff --git a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl
@@ -18,14 +18,14 @@ f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
 f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s))
 
 faces = (f1, f2, f3, f4)
-mapping = Trixi.transfinite_mapping(faces)
 
 # 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)
+                  faces = faces,
+                  initial_refinement_level = 1,
+                  periodicity = (true, true))
 
 # 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.

diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl
@@ -30,13 +30,8 @@ mapping_flag = Trixi.transfinite_mapping(faces)
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                            joinpath(@__DIR__, "square_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(2))
-
-mesh = T8codeMesh(conn, polydeg = 3,
-                  mapping = mapping_flag,
+mesh = T8codeMesh(mesh_file, 2;
+                  mapping = mapping_flag, polydeg = 3,
                   initial_refinement_level = 2)
 
 # A semidiscretization collects data structures and functions for the spatial discretization.

diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl
@@ -32,12 +32,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp",
                            joinpath(@__DIR__, "square_unstructured_1.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(2))
-
-mesh = T8codeMesh(conn, polydeg = 3,
+mesh = T8codeMesh(mesh_file, 2; polydeg = 3,
                   mapping = mapping,
                   initial_refinement_level = 1)
 

diff --git a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl
@@ -32,12 +32,7 @@ mapping_flag = Trixi.transfinite_mapping(faces)
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                            joinpath(@__DIR__, "square_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(2))
-
-mesh = T8codeMesh(conn, polydeg = 3,
+mesh = T8codeMesh(mesh_file, 2; polydeg = 3,
                   mapping = mapping_flag,
                   initial_refinement_level = 1)
 

diff --git a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl
@@ -70,12 +70,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                            joinpath(@__DIR__, "square_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(2))
-
-mesh = T8codeMesh(conn, polydeg = 4,
+mesh = T8codeMesh(mesh_file, 2; polydeg = 4,
                   mapping = mapping_twist,
                   initial_refinement_level = 1)
 

diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl
@@ -50,12 +50,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
                            joinpath(@__DIR__, "cube_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connectivity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
-mesh = T8codeMesh(conn, polydeg = 2,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 2,
                   mapping = mapping,
                   initial_refinement_level = 1)
 

diff --git a/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl
@@ -47,12 +47,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp",
                            joinpath(@__DIR__, "cube_unstructured_1.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connectivity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
-mesh = T8codeMesh(conn, polydeg = 3,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 3,
                   mapping = mapping,
                   initial_refinement_level = 2)
 

diff --git a/examples/t8code_3d_dgsem/elixir_euler_ec.jl b/examples/t8code_3d_dgsem/elixir_euler_ec.jl
@@ -47,12 +47,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
                            joinpath(@__DIR__, "cube_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connectivity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
-mesh = T8codeMesh(conn, polydeg = 5,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 5,
                   mapping = mapping,
                   initial_refinement_level = 0)
 

diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl
@@ -48,12 +48,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp",
                            joinpath(@__DIR__, "cube_unstructured_1.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connectivity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
-mesh = T8codeMesh(conn, polydeg = 2,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 2,
                   mapping = mapping,
                   initial_refinement_level = 0)
 

diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl
@@ -37,12 +37,7 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
                            joinpath(@__DIR__, "cube_unstructured_2.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
-mesh = T8codeMesh(conn, polydeg = 3,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 3,
                   mapping = mapping,
                   initial_refinement_level = 0)
 

diff --git a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl
@@ -50,13 +50,8 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp",
                            joinpath(@__DIR__, "cube_unstructured_1.inp"))
 
-# INP mesh files are only support by p4est. Hence, we
-# create a p4est connecvity object first from which
-# we can create a t8code mesh.
-conn = Trixi.read_inp_p4est(mesh_file, Val(3))
-
 # Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above).
-mesh = T8codeMesh(conn, polydeg = 2,
+mesh = T8codeMesh(mesh_file, 3; polydeg = 2,
                   mapping = mapping,
                   initial_refinement_level = 0)
 

diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl
@@ -387,12 +387,44 @@ function P4estMesh{NDIMS}(meshfile::String;
                             p4est_partition_allow_for_coarsening)
 end
 
+# Wrapper for `p4est_connectivity_from_hohqmesh_abaqus`. The latter is used
+# by `T8codeMesh`, too.
+function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level,
+                                         n_dimensions, RealT)
+    connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_hohqmesh_abaqus(meshfile,
+                                                                                                         initial_refinement_level,
+                                                                                                         n_dimensions,
+                                                                                                         RealT)
+
+    p4est = new_p4est(connectivity, initial_refinement_level)
+
+    return p4est, tree_node_coordinates, nodes, boundary_names
+end
+
+# Wrapper for `p4est_connectivity_from_standard_abaqus`. The latter is used
+# by `T8codeMesh`, too.
+function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
+                                         initial_refinement_level, n_dimensions, RealT,
+                                         boundary_symbols)
+    connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_standard_abaqus(meshfile,
+                                                                                                         mapping,
+                                                                                                         polydeg,
+                                                                                                         initial_refinement_level,
+                                                                                                         n_dimensions,
+                                                                                                         RealT,
+                                                                                                         boundary_symbols)
+
+    p4est = new_p4est(connectivity, initial_refinement_level)
+
+    return p4est, tree_node_coordinates, nodes, boundary_names
+end
+
 # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1]
 # and a list of boundary names for the `P4estMesh`. High-order boundary curve information as well as
 # the boundary names on each tree are provided by the `meshfile` created by
 # [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl).
-function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level,
-                                         n_dimensions, RealT)
+function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_level,
+                                                 n_dimensions, RealT)
     # Create the mesh connectivity using `p4est`
     connectivity = read_inp_p4est(meshfile, Val(n_dimensions))
     connectivity_pw = PointerWrapper(connectivity)
@@ -440,18 +472,17 @@ function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level,
         file_idx += 1
     end
 
-    p4est = new_p4est(connectivity, initial_refinement_level)
-
-    return p4est, tree_node_coordinates, nodes, boundary_names
+    return connectivity, tree_node_coordinates, nodes, boundary_names
 end
 
 # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1]
 # and a list of boundary names for the `P4estMesh`. The tree node coordinates are computed according to
 # the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. All boundary
 # names are given the name `:all`.
-function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
-                                         initial_refinement_level, n_dimensions, RealT,
-                                         boundary_symbols)
+function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg,
+                                                 initial_refinement_level, n_dimensions,
+                                                 RealT,
+                                                 boundary_symbols)
     # Create the mesh connectivity using `p4est`
     connectivity = read_inp_p4est(meshfile, Val(n_dimensions))
     connectivity_pw = PointerWrapper(connectivity)
@@ -474,8 +505,6 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
     calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices,
                                 tree_to_vertex)
 
-    p4est = new_p4est(connectivity, initial_refinement_level)
-
     if boundary_symbols === nothing
         # There's no simple and generic way to distinguish boundaries without any information given.
         # Name all of them :all.
@@ -495,7 +524,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
                                            Val(n_dimensions))
     end
 
-    return p4est, tree_node_coordinates, nodes, boundary_names
+    return connectivity, tree_node_coordinates, nodes, boundary_names
 end
 
 function parse_elements(meshfile, n_trees, n_dims)