diff --git a/Project.toml b/Project.toml index 68f9060b9d3..812d01733a4 100644 --- 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 index 0923e328487..9138586cccf 100644 --- 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 index a39f3a7e195..48f78dd6da3 100644 --- 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 index ba8f1b59b80..e512f328234 100644 --- 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 index 5e6c4193c50..d9d2c65d988 100644 --- 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 index e496eb76729..48684071d4b 100644 --- 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 index ff2e40ae607..592d5b15a85 100644 --- 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 index e7c0f4b7318..1f9aa3449b0 100644 --- 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 index ee27ee117fe..fe6aa48e7d9 100644 --- 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 index b720bfcd375..e1e4d850a86 100644 --- 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 index b70a6091adf..882e3aebebe 100644 --- 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 index 6ae38d20b5a..777cccf7ad7 100644 --- 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 index 6856be36ea1..a06e7927dd0 100644 --- 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 index abe9d9345b5..6bb98196231 100644 --- 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) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index cb2ac787e14..0af4c6ae023 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -26,7 +26,7 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: nmpiinterfaces :: Int nmpimortars :: Int - function T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, + function T8codeMesh{NDIMS}(forest::Ptr{t8_forest}, tree_node_coordinates, nodes, boundary_names, current_filename) where {NDIMS} is_parallel = mpi_isparallel() ? True() : False() @@ -100,6 +100,129 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::T8codeMesh) end end +""" + T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = 1, mapping = nothing) + +Main mesh constructor for the `T8codeMesh` wrapping around a given t8code +`forest` object. This constructor is typically called by other `T8codeMesh` +constructors. + +# Arguments +- `forest`: Pointer to a t8code forest. +- `boundary_names`: List of boundary names. +- `polydeg::Integer`: Polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. +- `mapping`: A function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +""" +function T8codeMesh{NDIMS, RealT}(forest::Ptr{t8_forest}, boundary_names; polydeg = 1, + mapping = nothing) where {NDIMS, RealT} + # In t8code reference space is [0,1]. + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = 0.5 .* (basis.nodes .+ 1.0) + + cmesh = t8_forest_get_cmesh(forest) + number_of_trees = t8_forest_get_num_global_trees(forest) + + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, + ntuple(_ -> length(nodes), NDIMS)..., + number_of_trees) + + reference_coordinates = Vector{Float64}(undef, 3) + + # Calculate node coordinates of reference mesh. + if NDIMS == 2 + number_of_corners = 4 # quadrilateral + + # Testing for negative element volumes. + vertices = zeros(3, number_of_corners) + for itree in 1:number_of_trees + vertices_pointer = t8_cmesh_get_tree_vertices(cmesh, itree - 1) + + # Note, `vertices = unsafe_wrap(Array, vertices_pointer, (3, 1 << NDIMS))` + # sometimes does not work since `vertices_pointer` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:number_of_corners + vertices[1, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 1) + vertices[2, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 2) + end + + # Check if tree's node ordering is right-handed or print a warning. + let z = zero(eltype(vertices)), o = one(eltype(vertices)) + u = vertices[:, 2] - vertices[:, 1] + v = vertices[:, 3] - vertices[:, 1] + w = [z, z, o] + + # Triple product gives signed volume of spanned parallelepiped. + vol = dot(cross(u, v), w) + + if vol < z + error("Discovered negative volumes in `cmesh`: vol = $vol") + end + end + + # Query geometry data from t8code. + for j in eachindex(nodes), i in eachindex(nodes) + reference_coordinates[1] = nodes[i] + reference_coordinates[2] = nodes[j] + reference_coordinates[3] = 0.0 + t8_geometry_evaluate(cmesh, itree - 1, reference_coordinates, 1, + @view(tree_node_coordinates[:, i, j, itree])) + end + end + + elseif NDIMS == 3 + number_of_corners = 8 # hexahedron + + # Testing for negative element volumes. + vertices = zeros(3, number_of_corners) + for itree in 1:number_of_trees + vertices_pointer = t8_cmesh_get_tree_vertices(cmesh, itree - 1) + + # Note, `vertices = unsafe_wrap(Array, vertices_pointer, (3, 1 << NDIMS))` + # sometimes does not work since `vertices_pointer` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:number_of_corners + vertices[1, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 1) + vertices[2, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 2) + vertices[3, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 3) + end + + # Check if tree's node ordering is right-handed or print a warning. + let z = zero(eltype(vertices)) + u = vertices[:, 2] - vertices[:, 1] + v = vertices[:, 3] - vertices[:, 1] + w = vertices[:, 5] - vertices[:, 1] + + # Triple product gives signed volume of spanned parallelepiped. + vol = dot(cross(u, v), w) + + if vol < z + error("Discovered negative volumes in `cmesh`: vol = $vol") + end + end + + # Query geometry data from t8code. + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + reference_coordinates[1] = nodes[i] + reference_coordinates[2] = nodes[j] + reference_coordinates[3] = nodes[k] + t8_geometry_evaluate(cmesh, itree - 1, reference_coordinates, 1, + @view(tree_node_coordinates[:, i, j, k, itree])) + end + end + else + throw(ArgumentError("$NDIMS dimensions are not supported.")) + end + + # Apply user defined mapping. + map_node_coordinates!(tree_node_coordinates, mapping) + + return T8codeMesh{NDIMS}(forest, tree_node_coordinates, basis.nodes, + boundary_names, "") +end + """ T8codeMesh(trees_per_dimension; polydeg, mapping=identity, RealT=Float64, initial_refinement_level=0, periodicity=true) @@ -187,57 +310,10 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, mpi_comm()) - basis = LobattoLegendreBasis(RealT, polydeg) - nodes = basis.nodes - - num_trees = t8_cmesh_get_num_trees(cmesh) - - tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, - ntuple(_ -> length(nodes), NDIMS)..., - num_trees) - - # Get cell length in reference mesh: Omega_ref = [-1,1]^NDIMS. - dx = [2 / n for n in trees_per_dimension] - # Non-periodic boundaries. boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) - if mapping === nothing - mapping_ = coordinates2mapping(ntuple(_ -> -1.0, NDIMS), ntuple(_ -> 1.0, NDIMS)) - else - mapping_ = mapping - end - - for itree in 1:num_trees - veptr = t8_cmesh_get_tree_vertices(cmesh, itree - 1) - verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) - - # Calculate node coordinates of reference mesh. - if NDIMS == 2 - cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] - cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] - - for j in eachindex(nodes), i in eachindex(nodes) - tree_node_coordinates[:, i, j, itree] .= mapping_(cell_x_offset + - dx[1] * nodes[i] / 2, - cell_y_offset + - dx[2] * nodes[j] / 2) - end - elseif NDIMS == 3 - cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] - cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] - cell_z_offset = (verts[3, 1] - 0.5 * (trees_per_dimension[3] - 1)) * dx[3] - - for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) - tree_node_coordinates[:, i, j, k, itree] .= mapping_(cell_x_offset + - dx[1] * nodes[i] / 2, - cell_y_offset + - dx[2] * nodes[j] / 2, - cell_z_offset + - dx[3] * nodes[k] / 2) - end - end - + for itree in 1:t8_forest_get_num_global_trees(forest) if !periodicity[1] boundary_names[1, itree] = :x_neg boundary_names[2, itree] = :x_pos @@ -256,8 +332,13 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, end end - return T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, - boundary_names, "") + # Note, `p*est_connectivity_new_brick` converts a domain of `[0,nx] x [0,ny] x ....`. + # Hence, transform mesh coordinates to reference space [-1,1]^NDIMS before applying user defined mapping. + mapping_(xyz...) = mapping((x * 2.0 / tpd - 1.0 for (x, tpd) in zip(xyz, + trees_per_dimension))...) + + return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = polydeg, + mapping = mapping_) end """ @@ -295,106 +376,11 @@ function T8codeMesh(cmesh::Ptr{t8_cmesh}; forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, mpi_comm()) - basis = LobattoLegendreBasis(RealT, polydeg) - nodes = basis.nodes - - num_trees = t8_cmesh_get_num_trees(cmesh) - - tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, - ntuple(_ -> length(nodes), NDIMS)..., - num_trees) - - nodes_in = [-1.0, 1.0] - matrix = polynomial_interpolation_matrix(nodes_in, nodes) - - num_local_trees = t8_cmesh_get_num_local_trees(cmesh) - - if NDIMS == 2 - data_in = Array{RealT, 3}(undef, 2, 2, 2) - tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) - verts = zeros(3, 4) - - for itree in 0:(num_local_trees - 1) - veptr = t8_cmesh_get_tree_vertices(cmesh, itree) - - # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` - # sometimes does not work since `veptr` is not necessarily properly - # aligned to 8 bytes. - for icorner in 1:4 - verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) - verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) - end - - # Check if tree's node ordering is right-handed or print a warning. - let z = zero(eltype(verts)), o = one(eltype(verts)) - u = verts[:, 2] - verts[:, 1] - v = verts[:, 3] - verts[:, 1] - w = [z, z, o] - - # Triple product gives signed volume of spanned parallelepiped. - vol = dot(cross(u, v), w) - - if vol < z - @warn "Discovered negative volumes in `cmesh`: vol = $vol" - end - end - - # Tree vertices are stored in z-order. - @views data_in[:, 1, 1] .= verts[1:2, 1] - @views data_in[:, 2, 1] .= verts[1:2, 2] - @views data_in[:, 1, 2] .= verts[1:2, 3] - @views data_in[:, 2, 2] .= verts[1:2, 4] - - # Interpolate corner coordinates to specified nodes. - multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, itree + 1), - matrix, matrix, - data_in, - tmp1) - end - - elseif NDIMS == 3 - data_in = Array{RealT, 4}(undef, 3, 2, 2, 2) - tmp1 = zeros(RealT, 3, length(nodes), length(nodes_in), length(nodes_in)) - verts = zeros(3, 8) - - for itree in 0:(num_trees - 1) - veptr = t8_cmesh_get_tree_vertices(cmesh, itree) - - # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` - # sometimes does not work since `veptr` is not necessarily properly - # aligned to 8 bytes. - for icorner in 1:8 - verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) - verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) - verts[3, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 3) - end - - # Tree vertices are stored in z-order. - @views data_in[:, 1, 1, 1] .= verts[1:3, 1] - @views data_in[:, 2, 1, 1] .= verts[1:3, 2] - @views data_in[:, 1, 2, 1] .= verts[1:3, 3] - @views data_in[:, 2, 2, 1] .= verts[1:3, 4] - - @views data_in[:, 1, 1, 2] .= verts[1:3, 5] - @views data_in[:, 2, 1, 2] .= verts[1:3, 6] - @views data_in[:, 1, 2, 2] .= verts[1:3, 7] - @views data_in[:, 2, 2, 2] .= verts[1:3, 8] - - # Interpolate corner coordinates to specified nodes. - multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, :, itree + 1), - matrix, matrix, matrix, - data_in, - tmp1) - end - end - - map_node_coordinates!(tree_node_coordinates, mapping) - - # There's no simple and generic way to distinguish boundaries. Name all of them :all. - boundary_names = fill(:all, 2 * NDIMS, num_trees) + # There's no simple and generic way to distinguish boundaries, yet. Name all of them :all. + boundary_names = fill(:all, 2 * NDIMS, t8_cmesh_get_num_trees(cmesh)) - return T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, - boundary_names, "") + return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = polydeg, + mapping = mapping) end """ @@ -445,36 +431,201 @@ function T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) return T8codeMesh(cmesh; kwargs...) end +# Convenience types for multiple dispatch. Only used in this file. +struct GmshFile{NDIMS} + path::String +end + +struct AbaqusFile{NDIMS} + path::String +end + """ - T8codeMesh(meshfile::String, ndims; kwargs...) + T8codeMesh(meshfile::String, NDIMS; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming -mesh from a Gmsh mesh file (`.msh`). +mesh from either a Gmsh mesh file (`.msh`) or Abaqus mesh file (`.inp`) which is determined +by the file extension. # Arguments -- `meshfile::String`: path to a Gmsh mesh file. +- `filepath::String`: path to a Gmsh or Abaqus mesh file. - `ndims`: Mesh file dimension: `2` or `3`. -- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + +# Optional Keyword Arguments +- `mapping`: A function of `NDIMS` variables to describe the mapping that transforms the imported mesh to the physical domain. Use `nothing` for the identity map. -- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. +- `polydeg::Integer`: Polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. The default of `1` creates an uncurved geometry. Use a higher value if the mapping will curve the imported uncurved mesh. -- `RealT::Type`: the type that should be used for coordinates. -- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +- `RealT::Type`: The type that should be used for coordinates. +- `initial_refinement_level::Integer`: Refine the mesh uniformly to this level before the simulation starts. """ -function T8codeMesh(meshfile::String, ndims; kwargs...) +function T8codeMesh(filepath::String, ndims; kwargs...) # Prevent `t8code` from crashing Julia if the file doesn't exist. - @assert isfile(meshfile) + @assert isfile(filepath) + + meshfile_prefix, meshfile_suffix = splitext(filepath) + + file_extension = lowercase(meshfile_suffix) + + if file_extension == ".msh" + return T8codeMesh(GmshFile{ndims}(filepath); kwargs...) + end + + if file_extension == ".inp" + return T8codeMesh(AbaqusFile{ndims}(filepath); kwargs...) + end + + throw(ArgumentError("Unknown file extension: " * file_extension)) +end - meshfile_prefix, meshfile_suffix = splitext(meshfile) +""" + T8codeMesh(meshfile::GmshFile{NDIMS}; kwargs...) - cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), ndims, 0, 0) +Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming +mesh from a Gmsh mesh file (`.msh`). + +# Arguments +- `meshfile::GmshFile{NDIMS}`: Gmsh mesh file object of dimension NDIMS and give `path` to the file. + +# Optional Keyword Arguments +- `mapping`: A function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: Polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: The type that should be used for coordinates. +- `initial_refinement_level::Integer`: Refine the mesh uniformly to this level before the simulation starts. +""" +function T8codeMesh(meshfile::GmshFile{NDIMS}; kwargs...) where {NDIMS} + # Prevent `t8code` from crashing Julia if the file doesn't exist. + @assert isfile(meshfile.path) + + meshfile_prefix, meshfile_suffix = splitext(meshfile.path) + + cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), NDIMS, 0, 0) return T8codeMesh(cmesh; kwargs...) end +""" + T8codeMesh(meshfile::AbaqusFile{NDIMS}; + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0, unsaved_changes=true, + boundary_symbols = nothing) + +Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming +mesh from an Abaqus mesh file (`.inp`). + +To create a curved unstructured `T8codeMesh` two strategies are available: + +- `HOHQMesh Abaqus`: High-order, curved boundary information created by + [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl) is + available in the `meshfile`. The mesh polynomial degree `polydeg` + of the boundaries is provided from the `meshfile`. The computation of + the mapped tree coordinates is done with transfinite interpolation + with linear blending similar to [`UnstructuredMesh2D`](@ref). Boundary name + information is also parsed from the `meshfile` such that different boundary + conditions can be set at each named boundary on a given tree. + +- `Standard Abaqus`: By default, with `mapping=nothing` and `polydeg=1`, this creates a + straight-sided mesh from the information parsed from the `meshfile`. If a mapping + function is specified then it computes the mapped tree coordinates via polynomial + interpolants with degree `polydeg`. The mesh created by this function will only + have one boundary `:all` if `boundary_symbols` is not specified. + If `boundary_symbols` is specified the mesh file will be parsed for nodesets defining + the boundary nodes from which boundary edges (2D) and faces (3D) will be assigned. + +Note that the `mapping` and `polydeg` keyword arguments are only used by the `HOHQMesh Abaqus` option. +The `Standard Abaqus` routine obtains the mesh `polydeg` directly from the `meshfile` +and constructs the transfinite mapping internally. + +The particular strategy is selected according to the header present in the `meshfile` where +the constructor checks whether or not the `meshfile` was created with +[HOHQMesh.jl](https://github.com/trixi-framework/HOHQMesh.jl). +If the Abaqus file header is not present in the `meshfile` then the `T8codeMesh` is created +by `Standard Abaqus`. + +The default keyword argument `initial_refinement_level=0` corresponds to a forest +where the number of trees is the same as the number of elements in the original `meshfile`. +Increasing the `initial_refinement_level` allows one to uniformly refine the base mesh given +in the `meshfile` to create a forest with more trees before the simulation begins. +For example, if a two-dimensional base mesh contains 25 elements then setting +`initial_refinement_level=1` creates an initial forest of `2^2 * 25 = 100` trees. + +# Arguments +- `meshfile::AbaqusFile{NDIMS}`: Abaqus mesh file object of dimension NDIMS and given `path` to the file. + +# Optional Keyword Arguments +- `mapping`: A function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: Polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: The type that should be used for coordinates. +- `initial_refinement_level::Integer`: Refine the mesh uniformly to this level before the simulation starts. +- `boundary_symbols::Vector{Symbol}`: A vector of symbols that correspond to the boundary names in the `meshfile`. + If `nothing` is passed then all boundaries are named `:all`. +""" +function T8codeMesh(meshfile::AbaqusFile{NDIMS}; + mapping = nothing, polydeg = 1, RealT = Float64, + initial_refinement_level = 0, + boundary_symbols = nothing) where {NDIMS} + # Prevent `t8code` from crashing Julia if the file doesn't exist. + @assert isfile(meshfile.path) + + # Read in the Header of the meshfile to determine which constructor is appropriate. + header = open(meshfile.path, "r") do io + readline(io) # Header of the Abaqus file; discarded + readline(io) # Read in the actual header information + end + + # Check if the meshfile was generated using HOHQMesh. + if header == " File created by HOHQMesh" + # Mesh curvature and boundary naming is handled with additional information available in meshfile + connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_hohqmesh_abaqus(meshfile.path, + initial_refinement_level, + NDIMS, + RealT) + # Apply user defined mapping. + map_node_coordinates!(tree_node_coordinates, mapping) + else + # Mesh curvature is handled directly by applying the mapping keyword argument. + connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_standard_abaqus(meshfile.path, + mapping, + polydeg, + initial_refinement_level, + NDIMS, + RealT, + boundary_symbols) + end + + cmesh = t8_cmesh_new_from_connectivity(connectivity, mpi_comm()) + p4est_connectivity_destroy(connectivity) + + do_face_ghost = mpi_isparallel() + scheme = t8_scheme_new_default_cxx() + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, + mpi_comm()) + + return T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, + boundary_names, "") +end + +function t8_cmesh_new_from_connectivity(connectivity::Ptr{p4est_connectivity}, comm) + return t8_cmesh_new_from_p4est(connectivity, comm, 0) +end + +function t8_cmesh_new_from_connectivity(connectivity::Ptr{p8est_connectivity}, comm) + return t8_cmesh_new_from_p8est(connectivity, comm, 0) +end + struct adapt_callback_passthrough adapt_callback::Function user_data::Any diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index d536a6dd73a..b63d2a105ac 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -30,8 +30,16 @@ mkdir(outdir) end end +@trixi_testset "test load mesh from path" begin + mktempdir() do path + @test_throws "Unknown file extension: .unknown_ext" begin + mesh = T8codeMesh(touch(joinpath(path, "dummy.unknown_ext")), 2) + end + end +end + @trixi_testset "test check_for_negative_volumes" begin - @test_warn "Discovered negative volumes" begin + @test_throws "Discovered negative volumes" begin # Unstructured mesh with six cells which have left-handed node ordering. mesh_file = Trixi.download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", joinpath(EXAMPLES_DIR, @@ -42,6 +50,27 @@ end end end +@trixi_testset "test t8code mesh from p4est connectivity" begin + @test begin + # Here we use the connectivity constructor from `P4est.jl` since the + # method dispatch works only on `Ptr{p4est_connectivity}` which + # actually is `Ptr{P4est.LibP4est.p4est_connectivity}`. + conn = Trixi.P4est.LibP4est.p4est_connectivity_new_brick(2, 3, 1, 1) + mesh = T8codeMesh(conn) + all(size(mesh.tree_node_coordinates) .== (2, 2, 2, 6)) + end +end + +@trixi_testset "test t8code mesh from ABAQUS HOHQMesh file" begin + @test begin + # Unstructured ABAQUS mesh file created with HOHQMesh.. + file_path = Trixi.download("https://gist.githubusercontent.com/jmark/9e0da4306e266617eeb19bc56b0e7feb/raw/e6856e1deb648a807f6bb6d6dcacff9e55d94e2a/round_2d_tank.inp", + joinpath(EXAMPLES_DIR, "round_2d_tank.inp")) + mesh = T8codeMesh(file_path, 2) + all(size(mesh.tree_node_coordinates) .== (2, 4, 4, 340)) + end +end + @trixi_testset "elixir_advection_basic.jl" begin # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index 4232cf04094..81d2a7cdd85 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -13,6 +13,17 @@ isdir(outdir) && rm(outdir, recursive = true) mkdir(outdir) @testset "T8codeMesh3D" begin + @trixi_testset "test t8code mesh from p8est connectivity" begin + @test begin + # Here we use the connectivity constructor from `P4est.jl` since the + # method dispatch works only on `Ptr{p8est_connectivity}` which + # actually is `Ptr{P4est.LibP4est.p8est_connectivity}`. + conn = Trixi.P4est.LibP4est.p8est_connectivity_new_brick(2, 3, 4, 1, 1, 1) + mesh = T8codeMesh(conn) + all(size(mesh.tree_node_coordinates) .== (3, 2, 2, 2, 24)) + end + end + # This test is identical to the one in `test_p4est_3d.jl`. @trixi_testset "elixir_advection_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), @@ -202,7 +213,7 @@ mkdir(outdir) 3.3228975127030935e-13, 9.592326932761353e-13, ], - tspan=(0.0, 0.1)) + tspan=(0.0, 0.1), atol=5.0e-13,) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let