From d1dfb76d53c9244defef34ae3a147bb0117a2b60 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 4 Jan 2024 14:46:39 +0100 Subject: [PATCH 01/69] first take on BCs for standard Abaqus --- .../p4est_2d_dgsem/elixir_euler_airfoil.jl | 57 ++++++ src/meshes/p4est_mesh.jl | 172 +++++++++++++++++- 2 files changed, 221 insertions(+), 8 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_airfoil.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl new file mode 100644 index 00000000000..4133433ff29 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl @@ -0,0 +1,57 @@ + +using Trixi +using OrdinaryDiffEq, Plots + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +d = 3 + +solver = DGSEM(polydeg = d, surface_flux = flux_lax_friedrichs) + +base_path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Airfoil_geo/" + +mesh_file = base_path * "NACA.inp" + +boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] + +mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols) +boundary_conditions = Dict(:PhysicalLine10 => BoundaryConditionDirichlet(initial_condition), + :PhysicalLine20 => BoundaryConditionDirichlet(initial_condition), + :PhysicalLine30 => BoundaryConditionDirichlet(initial_condition), + :PhysicalLine40 => BoundaryConditionDirichlet(initial_condition)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary + +plot(sol) + +pd = PlotData2D(sol) +plot(pd["rho"]) +plot!(getmesh(pd)) \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 60db285e04f..1d0ac915ae0 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -345,11 +345,14 @@ For example, if a two-dimensional base mesh contains 25 elements then setting - `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity independent of domain partitioning. Should be `false` for static meshes to permit more fine-grained partitioning. +- `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 P4estMesh{NDIMS}(meshfile::String; mapping = nothing, polydeg = 1, RealT = Float64, initial_refinement_level = 0, unsaved_changes = true, - p4est_partition_allow_for_coarsening = true) where {NDIMS} + p4est_partition_allow_for_coarsening = true, + boundary_symbols = nothing) where {NDIMS} # Prevent `p4est` from crashing Julia if the file doesn't exist @assert isfile(meshfile) @@ -367,13 +370,24 @@ function P4estMesh{NDIMS}(meshfile::String; NDIMS, RealT) else - # Mesh curvature is handled directly by applying the mapping keyword argument - p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, - mapping, - polydeg, - initial_refinement_level, - NDIMS, - RealT) + if boundary_symbols === nothing + # Mesh curvature is handled directly by applying the mapping keyword argument + p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, + mapping, + polydeg, + initial_refinement_level, + NDIMS, + RealT) + else + # Mesh curvature is handled directly by applying the mapping keyword argument + p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, + mapping, + polydeg, + initial_refinement_level, + NDIMS, + RealT, + boundary_symbols) + end end return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, @@ -475,6 +489,148 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, return p4est, tree_node_coordinates, nodes, boundary_names end +function parse_elements(meshfile, n_trees) + element_node_matrix = Matrix{Int64}(undef, n_trees, 4) + el_list_follows = false + + tree_id = 1 + open(meshfile, "r") do file + for line in eachline(file) + # Check if the line contains nodes assembled in special set, i.e., boundary + if startswith(line, "*ELEMENT, type=CPS4") + el_list_follows = true + elseif el_list_follows + content = split(line, ",") + if length(content) == 5 + content_int = parse.(Int64, content) + element_node_matrix[tree_id, :] = content_int[2:end] + tree_id += 1 + else + el_list_follows = false + end + end + end + end + + return element_node_matrix +end + +# TODO: NSET probably sufficient +function parse_sets(set_type , meshfile, boundary_symbols) + @assert set_type == "NSET" || set_type == "ELSET" "Only NSET and ELSET are supported" + nodes_dict = Dict{Symbol, Vector{Int64}}() + current_symbol = nothing + current_nodes = Int64[] + + other_set_type = set_type == "NSET" ? "ELSET" : "NSET" + + open(meshfile, "r") do file + for line in eachline(file) + # Check if the line contains nodes assembled in special set, i.e., boundary + if startswith(line, "*" * set_type * "," * set_type * "=") + # Safe the previous nodeset + if current_symbol !== nothing + nodes_dict[current_symbol] = current_nodes + end + current_symbol = Symbol(split(line, "=")[2]) + current_nodes = Int64[] + elseif startswith(line, "*" * other_set_type * "," * other_set_type * "=") + # Safe the previous nodeset + if current_symbol !== nothing + nodes_dict[current_symbol] = current_nodes + end + # Ignore the other set type + current_symbol = nothing + elseif current_symbol !== nothing # Read only if there was already a nodeset specified + append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1])) + end + end + if current_symbol !== nothing + nodes_dict[current_symbol] = current_nodes + end + end + + for key in keys(nodes_dict) + if key ∉ boundary_symbols + delete!(nodes_dict, key) + end + end + + return nodes_dict +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) + # Create the mesh connectivity using `p4est` + connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) + connectivity_pw = PointerWrapper(connectivity) + + # These need to be of the type Int for unsafe_wrap below to work + n_trees::Int = connectivity_pw.num_trees[] + n_vertices::Int = connectivity_pw.num_vertices[] + + vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) + tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, + (2^n_dimensions, n_trees)) + + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = basis.nodes + + tree_node_coordinates = Array{RealT, n_dimensions + 2}(undef, n_dimensions, + ntuple(_ -> length(nodes), + n_dimensions)..., + n_trees) + calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, + tree_to_vertex) + + p4est = new_p4est(connectivity, initial_refinement_level) + + node_set_dict = parse_sets("NSET", meshfile, boundary_symbols) + #println(node_set_dict) + + element_set_dict = parse_sets("ELSET", meshfile, boundary_symbols) + #println(element_set_dict) + + element_node_matrix = parse_elements(meshfile, n_trees) + + boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) + + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + for boundary in keys(node_set_dict) + # Check bottom edge + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] + # Bottom boundary is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check right edge + if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] + # Right boundary is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check top edge + if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + # Top boundary is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check left edge + if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] + # Left boundary is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end + end + end + + println(boundary_names) + + return p4est, tree_node_coordinates, nodes, boundary_names +end + """ P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; polydeg, RealT=Float64, From 23b2e84b363b185cf00c56fc30f4cd36767bd7e4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 4 Jan 2024 14:57:17 +0100 Subject: [PATCH 02/69] simplify --- src/meshes/p4est_mesh.jl | 29 ++++++++--------------------- 1 file changed, 8 insertions(+), 21 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 1d0ac915ae0..2964a713a2d 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -515,33 +515,24 @@ function parse_elements(meshfile, n_trees) return element_node_matrix end -# TODO: NSET probably sufficient -function parse_sets(set_type , meshfile, boundary_symbols) - @assert set_type == "NSET" || set_type == "ELSET" "Only NSET and ELSET are supported" +function parse_node_sets(meshfile, boundary_symbols) nodes_dict = Dict{Symbol, Vector{Int64}}() current_symbol = nothing current_nodes = Int64[] - other_set_type = set_type == "NSET" ? "ELSET" : "NSET" - open(meshfile, "r") do file for line in eachline(file) # Check if the line contains nodes assembled in special set, i.e., boundary - if startswith(line, "*" * set_type * "," * set_type * "=") + if startswith(line, "*NSET,NSET=") # Safe the previous nodeset if current_symbol !== nothing nodes_dict[current_symbol] = current_nodes end + # New nodeset current_symbol = Symbol(split(line, "=")[2]) current_nodes = Int64[] - elseif startswith(line, "*" * other_set_type * "," * other_set_type * "=") - # Safe the previous nodeset - if current_symbol !== nothing - nodes_dict[current_symbol] = current_nodes - end - # Ignore the other set type - current_symbol = nothing elseif current_symbol !== nothing # Read only if there was already a nodeset specified + # There is always a trailing comma, remove the corresponding empty string append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1])) end end @@ -550,6 +541,7 @@ function parse_sets(set_type , meshfile, boundary_symbols) end end + # Remove all nodesets that are not boundaries for key in keys(nodes_dict) if key ∉ boundary_symbols delete!(nodes_dict, key) @@ -590,12 +582,9 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, p4est = new_p4est(connectivity, initial_refinement_level) - node_set_dict = parse_sets("NSET", meshfile, boundary_symbols) - #println(node_set_dict) - - element_set_dict = parse_sets("ELSET", meshfile, boundary_symbols) - #println(element_set_dict) - + # Read in nodes belonging to boundaries + node_set_dict = parse_node_sets(meshfile, boundary_symbols) + # Read in all elements with associated nodes to specify the boundaries element_node_matrix = parse_elements(meshfile, n_trees) boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) @@ -626,8 +615,6 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, end end - println(boundary_names) - return p4est, tree_node_coordinates, nodes, boundary_names end From d86ea1e1648e2322f8e00941e9b1cc2448bcee9f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 4 Jan 2024 15:57:59 +0100 Subject: [PATCH 03/69] Couple comments, example --- .../p4est_2d_dgsem/elixir_euler_airfoil.jl | 57 --------- .../elixir_euler_airfoil_mach2.jl | 120 ++++++++++++++++++ src/meshes/p4est_mesh.jl | 19 ++- 3 files changed, 134 insertions(+), 62 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/elixir_euler_airfoil.jl create mode 100644 examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl deleted file mode 100644 index 4133433ff29..00000000000 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil.jl +++ /dev/null @@ -1,57 +0,0 @@ - -using Trixi -using OrdinaryDiffEq, Plots - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations2D(1.4) - -initial_condition = initial_condition_constant - -d = 3 - -solver = DGSEM(polydeg = d, surface_flux = flux_lax_friedrichs) - -base_path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Airfoil_geo/" - -mesh_file = base_path * "NACA.inp" - -boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] - -mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols) -boundary_conditions = Dict(:PhysicalLine10 => BoundaryConditionDirichlet(initial_condition), - :PhysicalLine20 => BoundaryConditionDirichlet(initial_condition), - :PhysicalLine30 => BoundaryConditionDirichlet(initial_condition), - :PhysicalLine40 => BoundaryConditionDirichlet(initial_condition)) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_conditions) - -tspan = (0.0, 0.5) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() - -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval) - -stepsize_callback = StepsizeCallback(cfl = 2.0) - -callbacks = CallbackSet(summary_callback, - analysis_callback, - stepsize_callback) - -############################################################################### -# run the simulation - -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), - dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); -summary_callback() # print the timer summary - -plot(sol) - -pd = PlotData2D(sol) -plot(pd["rho"]) -plot!(getmesh(pd)) \ No newline at end of file diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl new file mode 100644 index 00000000000..f2ba09728d9 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -0,0 +1,120 @@ + +using Trixi +using OrdinaryDiffEq, Plots +using Downloads: download + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +@inline function initial_condition_mach2_flow(x, t, equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 2.0 + v2 = 0.0 + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach2_flow + +# Supersonic inflow boundary condition. +# Calculate the boundary flux entirely from the external solution state, i.e., set +# external solution state values for everything entering the domain. +@inline function boundary_condition_supersonic_inflow(u_inner, + normal_direction::AbstractVector, + x, t, surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach2_flow(x, t, equations) + flux = Trixi.flux(u_boundary, normal_direction, equations) + + return flux +end + +# Supersonic outflow boundary condition. +# Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow +# except all the solution state values are set from the internal solution as everything leaves the domain +@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + flux = Trixi.flux(u_inner, normal_direction, equations) + + return flux +end + +d = 3 + +surface_flux = flux_lax_friedrichs # LOCAL Lax-Friedrichs = Rusanov flux +volume_flux = flux_ranocha # For advanced DG technique (flux splitting) + +basis = LobattoLegendreBasis(d) +shock_indicator = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +# DG Solver +solver = DGSEM(polydeg = d, surface_flux = surface_flux, + volume_integral = volume_integral) + +mesh_file = joinpath(@__DIR__, "NACA6412.inp") +isfile(mesh_file) || + download("https://gist.github.com/DanielDoehring/b434005800a1b0c0b4b50a8772283019/raw/1f6649b3370163d75d268176b51775f8685dd81d/NACA6412.inp", + mesh_file) + +boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] + +mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols) + +boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, + :PhysicalLine20 => boundary_condition_outflow, # Right boundary + :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil + :PhysicalLine40 => boundary_condition_outflow # Top/Bottom + ) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 4.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# positivity limiter necessary for this example with strong shocks. Very sensitive +# to the order of the limiter variables, pressure must come first. +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6), + variables = (pressure, Trixi.density)) + +# Run the simulation +############################################################################### +sol = solve(ode, #SSPRK104(; thread = OrdinaryDiffEq.True(), stage_limiter!); + SSPRK104(; thread = OrdinaryDiffEq.True()); + dt = 42.0, # overwritten + callback = callbacks); + +summary_callback() # print the timer summary + +plot(sol) + +pd = PlotData2D(sol) +plot(pd["rho"]) +plot!(getmesh(pd)) \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 2964a713a2d..e2de8414921 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -496,16 +496,18 @@ function parse_elements(meshfile, n_trees) tree_id = 1 open(meshfile, "r") do file for line in eachline(file) - # Check if the line contains nodes assembled in special set, i.e., boundary - if startswith(line, "*ELEMENT, type=CPS4") + # Check for quadrilateral elements + if startswith(line, "*ELEMENT, type=CPS4") || startswith(line, "*ELEMENT, type=C2D4") || + startswith(line, "*ELEMENT, type=S4") # TODO: 3D: C3D8 el_list_follows = true elseif el_list_follows content = split(line, ",") - if length(content) == 5 + if length(content) == 5 # Check that we still read in connectivity data content_int = parse.(Int64, content) - element_node_matrix[tree_id, :] = content_int[2:end] + # Add constituent nodes to the element_node_matrix + element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id tree_id += 1 - else + else # Read all elements for this ELSET el_list_follows = false end end @@ -591,6 +593,9 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, for tree in 1:n_trees tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 + # and search for "Node ordering and face numbering on elements" for boundary in keys(node_set_dict) # Check bottom edge if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] @@ -614,6 +619,10 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, end end end + # TODO: 3D + # Admitted 3D element: C3D8 + # See for node numbering: + # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html return p4est, tree_node_coordinates, nodes, boundary_names end From 3b2f545fa01b46aff1a0298e091d353af771dfd6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 4 Jan 2024 15:59:56 +0100 Subject: [PATCH 04/69] comments --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index f2ba09728d9..f2e97826880 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -47,8 +47,8 @@ end d = 3 -surface_flux = flux_lax_friedrichs # LOCAL Lax-Friedrichs = Rusanov flux -volume_flux = flux_ranocha # For advanced DG technique (flux splitting) +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha basis = LobattoLegendreBasis(d) shock_indicator = IndicatorHennemannGassner(equations, basis, @@ -99,7 +99,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -# positivity limiter necessary for this example with strong shocks. Very sensitive +# Positivity limiter might be necessary for examples with strong shocks. Very sensitive # to the order of the limiter variables, pressure must come first. stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6), variables = (pressure, Trixi.density)) From 27e332b566cd6445c90c74a78b166e6b8face0ae Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 08:50:52 +0100 Subject: [PATCH 05/69] shorten code --- src/meshes/p4est_mesh.jl | 150 ++++++++++++++------------------------- 1 file changed, 55 insertions(+), 95 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index e2de8414921..7cb9ce2b6ca 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -370,24 +370,14 @@ function P4estMesh{NDIMS}(meshfile::String; NDIMS, RealT) else - if boundary_symbols === nothing - # Mesh curvature is handled directly by applying the mapping keyword argument - p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, - mapping, - polydeg, - initial_refinement_level, - NDIMS, - RealT) - else - # Mesh curvature is handled directly by applying the mapping keyword argument - p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, - mapping, - polydeg, - initial_refinement_level, - NDIMS, - RealT, - boundary_symbols) - end + # Mesh curvature is handled directly by applying the mapping keyword argument + p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, + mapping, + polydeg, + initial_refinement_level, + NDIMS, + RealT, + boundary_symbols) end return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, @@ -458,7 +448,8 @@ end # 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) + 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) @@ -483,8 +474,51 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, p4est = new_p4est(connectivity, initial_refinement_level) - # There's no simple and generic way to distinguish boundaries. Name all of them :all. - boundary_names = fill(:all, 2 * n_dimensions, n_trees) + if boundary_symbols === nothing + # There's no simple and generic way to distinguish boundaries. Name all of them :all. + boundary_names = fill(:all, 2 * n_dimensions, n_trees) + else # Boundary information supplied + # Read in nodes belonging to boundaries + node_set_dict = parse_node_sets(meshfile, boundary_symbols) + # Read in all elements with associated nodes to specify the boundaries + element_node_matrix = parse_elements(meshfile, n_trees) + + # Initialize boundary information matrix with symbol for no boundary / internal connection + boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) + + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 + # and search for "Node ordering and face numbering on elements" + for boundary in keys(node_set_dict) + # Check bottom edge + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] + # Bottom boundary is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check right edge + if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] + # Right boundary is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check top edge + if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + # Top boundary is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check left edge + if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] + # Left boundary is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end + end + end + # TODO: 3D + # Admitted 3D element: C3D8 + # See for node numbering: + # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html + end return p4est, tree_node_coordinates, nodes, boundary_names end @@ -553,80 +587,6 @@ function parse_node_sets(meshfile, boundary_symbols) return nodes_dict 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) - # Create the mesh connectivity using `p4est` - connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) - connectivity_pw = PointerWrapper(connectivity) - - # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_pw.num_trees[] - n_vertices::Int = connectivity_pw.num_vertices[] - - vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) - tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, - (2^n_dimensions, n_trees)) - - basis = LobattoLegendreBasis(RealT, polydeg) - nodes = basis.nodes - - tree_node_coordinates = Array{RealT, n_dimensions + 2}(undef, n_dimensions, - ntuple(_ -> length(nodes), - n_dimensions)..., - n_trees) - calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, - tree_to_vertex) - - p4est = new_p4est(connectivity, initial_refinement_level) - - # Read in nodes belonging to boundaries - node_set_dict = parse_node_sets(meshfile, boundary_symbols) - # Read in all elements with associated nodes to specify the boundaries - element_node_matrix = parse_elements(meshfile, n_trees) - - boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) - - for tree in 1:n_trees - tree_nodes = element_node_matrix[tree, :] - # For node labeling, see - # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 - # and search for "Node ordering and face numbering on elements" - for boundary in keys(node_set_dict) - # Check bottom edge - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] - # Bottom boundary is position 3 in p4est indexing - boundary_names[3, tree] = boundary - end - # Check right edge - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] - # Right boundary is position 2 in p4est indexing - boundary_names[2, tree] = boundary - end - # Check top edge - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] - # Top boundary is position 4 in p4est indexing - boundary_names[4, tree] = boundary - end - # Check left edge - if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] - # Left boundary is position 1 in p4est indexing - boundary_names[1, tree] = boundary - end - end - end - # TODO: 3D - # Admitted 3D element: C3D8 - # See for node numbering: - # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html - - return p4est, tree_node_coordinates, nodes, boundary_names -end - """ P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; polydeg, RealT=Float64, From db6fb46659ea52b67ba18cfb924e968e717a3d49 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 09:44:21 +0100 Subject: [PATCH 06/69] refactor --- src/meshes/p4est_mesh.jl | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 7cb9ce2b6ca..45df4ba3784 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -481,7 +481,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # Read in nodes belonging to boundaries node_set_dict = parse_node_sets(meshfile, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries - element_node_matrix = parse_elements(meshfile, n_trees) + element_node_matrix = parse_elements(meshfile, n_trees, n_dimensions) # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) @@ -523,20 +523,24 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, return p4est, tree_node_coordinates, nodes, boundary_names end -function parse_elements(meshfile, n_trees) +function parse_elements(meshfile, n_trees, n_dims) element_node_matrix = Matrix{Int64}(undef, n_trees, 4) el_list_follows = false - tree_id = 1 + + @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" + element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : + ["*ELEMENT, type=C3D8"] + # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index + expected_content_length = n_dims == 2 ? 5 : 9 + open(meshfile, "r") do file for line in eachline(file) - # Check for quadrilateral elements - if startswith(line, "*ELEMENT, type=CPS4") || startswith(line, "*ELEMENT, type=C2D4") || - startswith(line, "*ELEMENT, type=S4") # TODO: 3D: C3D8 + if any(startswith(line, el_type) for el_type in element_types) el_list_follows = true elseif el_list_follows content = split(line, ",") - if length(content) == 5 # Check that we still read in connectivity data + if length(content) == expected_content_length # Check that we still read in connectivity data content_int = parse.(Int64, content) # Add constituent nodes to the element_node_matrix element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id @@ -564,6 +568,12 @@ function parse_node_sets(meshfile, boundary_symbols) if current_symbol !== nothing nodes_dict[current_symbol] = current_nodes end + + # Read only boundary node sets + if current_symbol ∉ boundary_symbols + current_symbol = nothing + end + # New nodeset current_symbol = Symbol(split(line, "=")[2]) current_nodes = Int64[] @@ -576,14 +586,14 @@ function parse_node_sets(meshfile, boundary_symbols) nodes_dict[current_symbol] = current_nodes end end - + #= # Remove all nodesets that are not boundaries for key in keys(nodes_dict) if key ∉ boundary_symbols delete!(nodes_dict, key) end end - + =# return nodes_dict end From 244f75fdda1f91a9cf5f24f4a1a2b1b538f1267f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 09:47:13 +0100 Subject: [PATCH 07/69] fmt --- src/meshes/p4est_mesh.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 45df4ba3784..07474c16842 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -372,12 +372,12 @@ function P4estMesh{NDIMS}(meshfile::String; else # Mesh curvature is handled directly by applying the mapping keyword argument p4est, tree_node_coordinates, nodes, boundary_names = p4est_mesh_from_standard_abaqus(meshfile, - mapping, - polydeg, - initial_refinement_level, - NDIMS, - RealT, - boundary_symbols) + mapping, + polydeg, + initial_refinement_level, + NDIMS, + RealT, + boundary_symbols) end return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, From 8f0e184428016aa419e2e1258344c7f0049032db Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 09:55:09 +0100 Subject: [PATCH 08/69] refactor --- src/meshes/p4est_mesh.jl | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 07474c16842..d5aee879e1c 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -569,14 +569,13 @@ function parse_node_sets(meshfile, boundary_symbols) nodes_dict[current_symbol] = current_nodes end - # Read only boundary node sets - if current_symbol ∉ boundary_symbols + if current_symbol in boundary_symbols + # New nodeset + current_symbol = Symbol(split(line, "=")[2]) + current_nodes = Int64[] + else # Read only boundary node sets current_symbol = nothing end - - # New nodeset - current_symbol = Symbol(split(line, "=")[2]) - current_nodes = Int64[] elseif current_symbol !== nothing # Read only if there was already a nodeset specified # There is always a trailing comma, remove the corresponding empty string append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1])) @@ -586,14 +585,7 @@ function parse_node_sets(meshfile, boundary_symbols) nodes_dict[current_symbol] = current_nodes end end - #= - # Remove all nodesets that are not boundaries - for key in keys(nodes_dict) - if key ∉ boundary_symbols - delete!(nodes_dict, key) - end - end - =# + return nodes_dict end From c783247c3ffa833181f8069149aafc2ba13cfd6d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 11:06:43 +0100 Subject: [PATCH 09/69] 3D --- .../elixir_euler_free_stream_boundaries.jl | 60 ++++++++++ src/meshes/p4est_mesh.jl | 109 ++++++++++++------ 2 files changed, 136 insertions(+), 33 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl new file mode 100644 index 00000000000..78e14b33679 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl @@ -0,0 +1,60 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +d = 2 +solver = DGSEM(polydeg = d, surface_flux = flux_lax_friedrichs) + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "cube_boundaries.inp") +isfile(default_mesh_file) || + download("https://gist.github.com/DanielDoehring/eefe73ae5d113bc3944a518b6e88e663/raw/359a58a808790f3c3efc050273270eb1cc8ee353/cube_boundaries.inp", + default_mesh_file) +mesh_file = default_mesh_file + +boundary_symbols = [:PhysicalSurface1] + +# Map the unstructured mesh with the mapping above +mesh = P4estMesh{3}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) + +boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index d5aee879e1c..16bcc38bd27 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -486,54 +486,97 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) - for tree in 1:n_trees - tree_nodes = element_node_matrix[tree, :] - # For node labeling, see - # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 - # and search for "Node ordering and face numbering on elements" - for boundary in keys(node_set_dict) - # Check bottom edge - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] - # Bottom boundary is position 3 in p4est indexing - boundary_names[3, tree] = boundary - end - # Check right edge - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] - # Right boundary is position 2 in p4est indexing - boundary_names[2, tree] = boundary - end - # Check top edge - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] - # Top boundary is position 4 in p4est indexing - boundary_names[4, tree] = boundary + if n_dimensions == 2 + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 + # and search for "Node ordering and face numbering on elements" + for boundary in keys(node_set_dict) + # Check bottom edge + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] + # Bottom boundary is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check right edge + if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] + # Right boundary is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check top edge + if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + # Top boundary is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check left edge + if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] + # Left boundary is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end end - # Check left edge - if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] - # Left boundary is position 1 in p4est indexing - boundary_names[1, tree] = boundary + end + else # n_dimensions == 3 + # Admitted 3D element: C3D8 + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html + for boundary in keys(node_set_dict) + # Check "front face" (y_min) + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] + # Front face is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check "back face" (y_max) + if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + # Front face is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check "left face" (x_min) + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + # Left face is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end + # Check "right face" (x_max) + if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && tree_nodes[7] in node_set_dict[boundary] + # Right face is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check "bottom face" (z_min) + if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + # Bottom face is position 5 in p4est indexing + boundary_names[5, tree] = boundary + end + # Check "top face" (z_max) + if tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + # Top face is position 6 in p4est indexing + boundary_names[6, tree] = boundary + end end end end - # TODO: 3D - # Admitted 3D element: C3D8 - # See for node numbering: - # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html end return p4est, tree_node_coordinates, nodes, boundary_names end function parse_elements(meshfile, n_trees, n_dims) - element_node_matrix = Matrix{Int64}(undef, n_trees, 4) - el_list_follows = false - tree_id = 1 - @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : ["*ELEMENT, type=C3D8"] # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 + element_node_matrix = Matrix{Int64}(undef, n_trees, expected_content_length - 1) + el_list_follows = false + tree_id = 1 + open(meshfile, "r") do file for line in eachline(file) if any(startswith(line, el_type) for el_type in element_types) @@ -569,9 +612,9 @@ function parse_node_sets(meshfile, boundary_symbols) nodes_dict[current_symbol] = current_nodes end + current_symbol = Symbol(split(line, "=")[2]) if current_symbol in boundary_symbols # New nodeset - current_symbol = Symbol(split(line, "=")[2]) current_nodes = Int64[] else # Read only boundary node sets current_symbol = nothing From 06fb60c75ea269e62cd1402bca7d261f9dc39ca2 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 11:19:20 +0100 Subject: [PATCH 10/69] tests and fmt --- .../elixir_euler_airfoil_mach2.jl | 54 ++++++---------- src/meshes/p4est_mesh.jl | 61 +++++++++++++------ test/test_p4est_2d.jl | 21 +++++++ test/test_p4est_3d.jl | 23 +++++++ 4 files changed, 104 insertions(+), 55 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index f2e97826880..55efabcd346 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -1,6 +1,6 @@ using Trixi -using OrdinaryDiffEq, Plots +using OrdinaryDiffEq using Downloads: download ############################################################################### @@ -9,14 +9,14 @@ using Downloads: download equations = CompressibleEulerEquations2D(1.4) @inline function initial_condition_mach2_flow(x, t, equations::CompressibleEulerEquations2D) - # set the freestream flow parameters - rho_freestream = 1.4 - v1 = 2.0 - v2 = 0.0 - p_freestream = 1.0 - - prim = SVector(rho_freestream, v1, v2, p_freestream) - return prim2cons(prim, equations) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 2.0 + v2 = 0.0 + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) end initial_condition = initial_condition_mach2_flow @@ -28,10 +28,10 @@ initial_condition = initial_condition_mach2_flow normal_direction::AbstractVector, x, t, surface_flux_function, equations::CompressibleEulerEquations2D) - u_boundary = initial_condition_mach2_flow(x, t, equations) - flux = Trixi.flux(u_boundary, normal_direction, equations) + u_boundary = initial_condition_mach2_flow(x, t, equations) + flux = Trixi.flux(u_boundary, normal_direction, equations) - return flux + return flux end # Supersonic outflow boundary condition. @@ -40,15 +40,15 @@ end @inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::CompressibleEulerEquations2D) - flux = Trixi.flux(u_inner, normal_direction, equations) + flux = Trixi.flux(u_inner, normal_direction, equations) - return flux + return flux end d = 3 surface_flux = flux_lax_friedrichs -volume_flux = flux_ranocha +volume_flux = flux_ranocha basis = LobattoLegendreBasis(d) shock_indicator = IndicatorHennemannGassner(equations, basis, @@ -67,17 +67,16 @@ solver = DGSEM(polydeg = d, surface_flux = surface_flux, mesh_file = joinpath(@__DIR__, "NACA6412.inp") isfile(mesh_file) || download("https://gist.github.com/DanielDoehring/b434005800a1b0c0b4b50a8772283019/raw/1f6649b3370163d75d268176b51775f8685dd81d/NACA6412.inp", - mesh_file) + mesh_file) boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] -mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols) +mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, :PhysicalLine20 => boundary_condition_outflow, # Right boundary :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil - :PhysicalLine40 => boundary_condition_outflow # Top/Bottom - ) + :PhysicalLine40 => boundary_condition_outflow) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) @@ -96,25 +95,10 @@ callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) -############################################################################### -# run the simulation - -# Positivity limiter might be necessary for examples with strong shocks. Very sensitive -# to the order of the limiter variables, pressure must come first. -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6), - variables = (pressure, Trixi.density)) - # Run the simulation ############################################################################### -sol = solve(ode, #SSPRK104(; thread = OrdinaryDiffEq.True(), stage_limiter!); - SSPRK104(; thread = OrdinaryDiffEq.True()); +sol = solve(ode, SSPRK104(; thread = OrdinaryDiffEq.True()); dt = 42.0, # overwritten callback = callbacks); summary_callback() # print the timer summary - -plot(sol) - -pd = PlotData2D(sol) -plot(pd["rho"]) -plot!(getmesh(pd)) \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 16bcc38bd27..8138ca6c1af 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -482,7 +482,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, node_set_dict = parse_node_sets(meshfile, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries element_node_matrix = parse_elements(meshfile, n_trees, n_dimensions) - + # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) @@ -494,22 +494,26 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # and search for "Node ordering and face numbering on elements" for boundary in keys(node_set_dict) # Check bottom edge - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] # Bottom boundary is position 3 in p4est indexing boundary_names[3, tree] = boundary end # Check right edge - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] # Right boundary is position 2 in p4est indexing boundary_names[2, tree] = boundary end # Check top edge - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] # Top boundary is position 4 in p4est indexing boundary_names[4, tree] = boundary end # Check left edge - if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] + if tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[1] in node_set_dict[boundary] # Left boundary is position 1 in p4est indexing boundary_names[1, tree] = boundary end @@ -523,38 +527,50 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html for boundary in keys(node_set_dict) # Check "front face" (y_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] # Front face is position 3 in p4est indexing boundary_names[3, tree] = boundary end # Check "back face" (y_max) - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Front face is position 4 in p4est indexing boundary_names[4, tree] = boundary end # Check "left face" (x_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Left face is position 1 in p4est indexing boundary_names[1, tree] = boundary end # Check "right face" (x_max) - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[6] in node_set_dict[boundary] && tree_nodes[7] in node_set_dict[boundary] + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] # Right face is position 2 in p4est indexing boundary_names[2, tree] = boundary end # Check "bottom face" (z_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] # Bottom face is position 5 in p4est indexing boundary_names[5, tree] = boundary end # Check "top face" (z_max) - if tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Top face is position 6 in p4est indexing boundary_names[6, tree] = boundary end @@ -568,8 +584,13 @@ end function parse_elements(meshfile, n_trees, n_dims) @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" - element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : - ["*ELEMENT, type=C3D8"] + element_types = n_dims == 2 ? + [ + "*ELEMENT, type=CPS4", + "*ELEMENT, type=C2D4", + "*ELEMENT, type=S4", + ] : + ["*ELEMENT, type=C3D8"] # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 @@ -621,7 +642,7 @@ function parse_node_sets(meshfile, boundary_symbols) end elseif current_symbol !== nothing # Read only if there was already a nodeset specified # There is always a trailing comma, remove the corresponding empty string - append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1])) + append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)])) end end if current_symbol !== nothing diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cebc2917d52..d921090be77 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -391,6 +391,27 @@ end end end +@trixi_testset "elixir_euler_airfoil_mach2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_airfoil_mach2.jl"), + l2=[ + 1.682112049215034e-10, 3.439606229463632e-10, + 1.973815936404324e-10, 7.666245265417042e-10, + ], + linf=[ + 2.636251394960709e-8, 4.4272578936244145e-8, + 2.8049958258193942e-8, 1.0125492977408612e-7, + ], + tspan=(0.0, 0.1)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_eulergravity_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"), l2=[ diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 4a2d2112c99..6596ecd911f 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -234,6 +234,29 @@ end end end +@trixi_testset "elixir_euler_free_stream_boundaries.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_boundaries.jl"), + l2=[ + 3.3824205900551404e-16, 5.761375062854377e-16, + 1.0712466724440564e-15, 1.2802903179864802e-15, + 2.310819174832704e-15, + ], + linf=[ + 2.1094237467877974e-15, 3.3584246494910985e-15, + 9.2148511043888e-15, 1.2101430968414206e-14, + 1.9539925233402755e-14, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), l2=[ From 4a8ff77b4bfe40221f2eca1f3adf2501f9c389c7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 11:29:20 +0100 Subject: [PATCH 11/69] comment --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index 55efabcd346..55005704d1a 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -73,10 +73,10 @@ boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :Physical mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) -boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, +boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary :PhysicalLine20 => boundary_condition_outflow, # Right boundary :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil - :PhysicalLine40 => boundary_condition_outflow) + :PhysicalLine40 => boundary_condition_outflow) # Top and bottom boundary semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) From 82fed9c29a6498417b212b477ec273107b839e88 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 11:33:33 +0100 Subject: [PATCH 12/69] news --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index cf695912ed7..c122e442099 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,8 @@ for human readability. #### Added - AMR for hyperbolic-parabolic equations on 3D `P4estMesh` - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` +- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, + can now be digested by Trixi in 2D and 3D ## Changes when updating to v0.6 from v0.5.x From a23e73ff6cbdf1cb07bb972537d8d07dcbf6f4e7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 11:35:39 +0100 Subject: [PATCH 13/69] comment --- src/meshes/p4est_mesh.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 8138ca6c1af..0afa30d7b80 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -645,6 +645,7 @@ function parse_node_sets(meshfile, boundary_symbols) append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)])) end end + # Safe the previous nodeset if current_symbol !== nothing nodes_dict[current_symbol] = current_nodes end From 46398cb179d8c13296720f4864158a3385b7596f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 14:00:23 +0100 Subject: [PATCH 14/69] stick to polydeg --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index 55005704d1a..e4e68cf0b01 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -45,7 +45,7 @@ end return flux end -d = 3 +polydeg = 3 surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha @@ -61,7 +61,7 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; volume_flux_fv = surface_flux) # DG Solver -solver = DGSEM(polydeg = d, surface_flux = surface_flux, +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) mesh_file = joinpath(@__DIR__, "NACA6412.inp") @@ -71,7 +71,7 @@ isfile(mesh_file) || boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] -mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) +mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary :PhysicalLine20 => boundary_condition_outflow, # Right boundary From 504022ae531726cf8d4e4fb846db0ee09844d08a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 14:04:45 +0100 Subject: [PATCH 15/69] polydeg 3 --- .../elixir_euler_free_stream_boundaries.jl | 8 ++++---- test/test_p4est_3d.jl | 12 ++++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl index 78e14b33679..cc06d61f31a 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl @@ -10,8 +10,8 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition = initial_condition_constant -d = 2 -solver = DGSEM(polydeg = d, surface_flux = flux_lax_friedrichs) +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) ############################################################################### # Get the uncurved mesh from a file (downloads the file if not available locally) @@ -25,7 +25,7 @@ mesh_file = default_mesh_file boundary_symbols = [:PhysicalSurface1] # Map the unstructured mesh with the mapping above -mesh = P4estMesh{3}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) +mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition)) @@ -45,7 +45,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 6596ecd911f..00f7d4e7cc9 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -238,14 +238,14 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_boundaries.jl"), l2=[ - 3.3824205900551404e-16, 5.761375062854377e-16, - 1.0712466724440564e-15, 1.2802903179864802e-15, - 2.310819174832704e-15, + 6.530157034651212e-16, 1.6057829680004379e-15, + 3.31107455378537e-15, 3.908829498281281e-15, + 5.048390610424672e-15, ], linf=[ - 2.1094237467877974e-15, 3.3584246494910985e-15, - 9.2148511043888e-15, 1.2101430968414206e-14, - 1.9539925233402755e-14, + 4.884981308350689e-15, 1.1921019726912618e-14, + 1.5432100042289676e-14, 2.298161660974074e-14, + 6.039613253960852e-14, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From a28314fb35b0ddb54997561dd781b53f39684ad2 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 5 Jan 2024 14:23:15 +0100 Subject: [PATCH 16/69] polydeg --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index e4e68cf0b01..e2611704296 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -50,7 +50,7 @@ polydeg = 3 surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha -basis = LobattoLegendreBasis(d) +basis = LobattoLegendreBasis(polydeg) shock_indicator = IndicatorHennemannGassner(equations, basis, alpha_max = 0.5, alpha_min = 0.001, From 3197aa5bb6804145b50dcc3a3765abe3f39aaedf Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 6 Jan 2024 10:34:28 +0100 Subject: [PATCH 17/69] Update examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl Co-authored-by: Andrew Winters --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index e2611704296..975015bb067 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -37,7 +37,7 @@ end # Supersonic outflow boundary condition. # Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow # except all the solution state values are set from the internal solution as everything leaves the domain -@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, +@inline function boundary_condition_supersonic_outflow(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::CompressibleEulerEquations2D) flux = Trixi.flux(u_inner, normal_direction, equations) From b84cfa9c98f60ed9c474f1c07c9385fec3d80780 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 6 Jan 2024 10:34:40 +0100 Subject: [PATCH 18/69] Update examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl Co-authored-by: Andrew Winters --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index 975015bb067..9d557965881 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -74,9 +74,9 @@ boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :Physical mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary - :PhysicalLine20 => boundary_condition_outflow, # Right boundary + :PhysicalLine20 => boundary_condition_supersonic_outflow, # Right boundary :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil - :PhysicalLine40 => boundary_condition_outflow) # Top and bottom boundary + :PhysicalLine40 => boundary_condition_supersonic_outflow) # Top and bottom boundary semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) From ad5ff149f24280367fd4b281f96bb6731e57f1ed Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 6 Jan 2024 10:36:16 +0100 Subject: [PATCH 19/69] Update src/meshes/p4est_mesh.jl Co-authored-by: Andrew Winters --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 0afa30d7b80..4a899d0e3ea 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -626,7 +626,7 @@ function parse_node_sets(meshfile, boundary_symbols) open(meshfile, "r") do file for line in eachline(file) - # Check if the line contains nodes assembled in special set, i.e., boundary + # Check if the line contains nodes assembled in a special set, i.e., a physical boundary if startswith(line, "*NSET,NSET=") # Safe the previous nodeset if current_symbol !== nothing From b63935bab97cd2f4b2cb74b129211e86a31fdcb2 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 6 Jan 2024 10:38:01 +0100 Subject: [PATCH 20/69] Update NEWS.md Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index c122e442099..3a3a504a911 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,7 +10,7 @@ for human readability. - AMR for hyperbolic-parabolic equations on 3D `P4estMesh` - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, - can now be digested by Trixi in 2D and 3D + can now be digested by Trixi in 2D and 3D. ## Changes when updating to v0.6 from v0.5.x From 02ccb9cb78202ec02d1fe82bdadacd6e619ee84b Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 6 Jan 2024 18:49:45 +0100 Subject: [PATCH 21/69] Update examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl --- examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl index cc06d61f31a..20d374da0f3 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl @@ -24,7 +24,6 @@ mesh_file = default_mesh_file boundary_symbols = [:PhysicalSurface1] -# Map the unstructured mesh with the mapping above mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition)) From c18ad8eb26af06a2fc33e46a30bd0c89ec97b5c7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 10:36:53 +0100 Subject: [PATCH 22/69] improve (?) formatting --- src/meshes/p4est_mesh.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 4a899d0e3ea..aabc1f7f37b 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -289,7 +289,8 @@ end P4estMesh{NDIMS}(meshfile::String; mapping=nothing, polydeg=1, RealT=Float64, initial_refinement_level=0, unsaved_changes=true, - p4est_partition_allow_for_coarsening=true) + p4est_partition_allow_for_coarsening=true, + boundary_symbols = nothing) Main mesh constructor for the `P4estMesh` that imports an unstructured, conforming mesh from an Abaqus mesh file (`.inp`). Each element of the conforming mesh parsed @@ -584,13 +585,10 @@ end function parse_elements(meshfile, n_trees, n_dims) @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" + # Valid element types (that can be processed by p4est) based on dimension element_types = n_dims == 2 ? - [ - "*ELEMENT, type=CPS4", - "*ELEMENT, type=C2D4", - "*ELEMENT, type=S4", - ] : - ["*ELEMENT, type=C3D8"] + ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", + "*ELEMENT, type=S4"] : ["*ELEMENT, type=C3D8"] # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 From 99b96a696b31800a22143fd5e915d7f8c1295695 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 10:37:51 +0100 Subject: [PATCH 23/69] fmt --- examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index 9d557965881..d893482b855 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -37,9 +37,11 @@ end # Supersonic outflow boundary condition. # Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow # except all the solution state values are set from the internal solution as everything leaves the domain -@inline function boundary_condition_supersonic_outflow(u_inner, normal_direction::AbstractVector, x, t, - surface_flux_function, - equations::CompressibleEulerEquations2D) +@inline function boundary_condition_supersonic_outflow(u_inner, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) flux = Trixi.flux(u_inner, normal_direction, equations) return flux From 918a303d0f89966776ad0d0b3d20ed449330188f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 11:02:32 +0100 Subject: [PATCH 24/69] split constructor --- .../elixir_euler_airfoil_mach2.jl | 2 +- src/meshes/p4est_mesh.jl | 202 ++++++++++-------- 2 files changed, 110 insertions(+), 94 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index d893482b855..88799720029 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -100,7 +100,7 @@ callbacks = CallbackSet(summary_callback, # Run the simulation ############################################################################### sol = solve(ode, SSPRK104(; thread = OrdinaryDiffEq.True()); - dt = 42.0, # overwritten + dt = 1.0, # overwritten by the `stepsize_callback` callback = callbacks); summary_callback() # print the timer summary diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index aabc1f7f37b..0d7cecfb569 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -476,7 +476,8 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, p4est = new_p4est(connectivity, initial_refinement_level) if boundary_symbols === nothing - # There's no simple and generic way to distinguish boundaries. Name all of them :all. + # There's no simple and generic way to distinguish boundaries without any information given. + # Name all of them :all. boundary_names = fill(:all, 2 * n_dimensions, n_trees) else # Boundary information supplied # Read in nodes belonging to boundaries @@ -487,97 +488,10 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) - if n_dimensions == 2 - for tree in 1:n_trees - tree_nodes = element_node_matrix[tree, :] - # For node labeling, see - # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 - # and search for "Node ordering and face numbering on elements" - for boundary in keys(node_set_dict) - # Check bottom edge - if tree_nodes[1] in node_set_dict[boundary] && - tree_nodes[2] in node_set_dict[boundary] - # Bottom boundary is position 3 in p4est indexing - boundary_names[3, tree] = boundary - end - # Check right edge - if tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[3] in node_set_dict[boundary] - # Right boundary is position 2 in p4est indexing - boundary_names[2, tree] = boundary - end - # Check top edge - if tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[4] in node_set_dict[boundary] - # Top boundary is position 4 in p4est indexing - boundary_names[4, tree] = boundary - end - # Check left edge - if tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[1] in node_set_dict[boundary] - # Left boundary is position 1 in p4est indexing - boundary_names[1, tree] = boundary - end - end - end - else # n_dimensions == 3 - # Admitted 3D element: C3D8 - for tree in 1:n_trees - tree_nodes = element_node_matrix[tree, :] - # For node labeling, see - # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html - for boundary in keys(node_set_dict) - # Check "front face" (y_min) - if tree_nodes[1] in node_set_dict[boundary] && - tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && - tree_nodes[6] in node_set_dict[boundary] - # Front face is position 3 in p4est indexing - boundary_names[3, tree] = boundary - end - # Check "back face" (y_max) - if tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && - tree_nodes[8] in node_set_dict[boundary] - # Front face is position 4 in p4est indexing - boundary_names[4, tree] = boundary - end - # Check "left face" (x_min) - if tree_nodes[1] in node_set_dict[boundary] && - tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && - tree_nodes[8] in node_set_dict[boundary] - # Left face is position 1 in p4est indexing - boundary_names[1, tree] = boundary - end - # Check "right face" (x_max) - if tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[6] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] - # Right face is position 2 in p4est indexing - boundary_names[2, tree] = boundary - end - # Check "bottom face" (z_min) - if tree_nodes[1] in node_set_dict[boundary] && - tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[4] in node_set_dict[boundary] - # Bottom face is position 5 in p4est indexing - boundary_names[5, tree] = boundary - end - # Check "top face" (z_max) - if tree_nodes[5] in node_set_dict[boundary] && - tree_nodes[6] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && - tree_nodes[8] in node_set_dict[boundary] - # Top face is position 6 in p4est indexing - boundary_names[6, tree] = boundary - end - end - end - end + # Fill `boundary_names` such that it can be processed by p4est + assign_boundaries_standard_abaqus!(boundary_names, n_trees, + element_node_matrix, node_set_dict, + Val(n_dimensions)) end return p4est, tree_node_coordinates, nodes, boundary_names @@ -604,7 +518,8 @@ function parse_elements(meshfile, n_trees, n_dims) content = split(line, ",") if length(content) == expected_content_length # Check that we still read in connectivity data content_int = parse.(Int64, content) - # Add constituent nodes to the element_node_matrix + # Add constituent nodes to the element_node_matrix. + # Important: Do not use idex form the Abaqus file, but the one from p4est element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id tree_id += 1 else # Read all elements for this ELSET @@ -652,6 +567,107 @@ function parse_node_sets(meshfile, boundary_symbols) return nodes_dict end +function assign_boundaries_standard_abaqus!(boundary_names, n_trees, + element_node_matrix, node_set_dict, + ::Val{2}) # 2D version + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 + # and search for "Node ordering and face numbering on elements" + for boundary in keys(node_set_dict) + # Check bottom edge + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] + # Bottom boundary is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check right edge + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] + # Right boundary is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check top edge + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] + # Top boundary is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check left edge + if tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[1] in node_set_dict[boundary] + # Left boundary is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end + end + end + + return boundary_names +end + +function assign_boundaries_standard_abaqus!(boundary_names, n_trees, + element_node_matrix, node_set_dict, + ::Val{3}) # 3D version + for tree in 1:n_trees + tree_nodes = element_node_matrix[tree, :] + # For node labeling, see + # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html + for boundary in keys(node_set_dict) + # Check "front face" (y_min) + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] + # Front face is position 3 in p4est indexing + boundary_names[3, tree] = boundary + end + # Check "back face" (y_max) + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] + # Front face is position 4 in p4est indexing + boundary_names[4, tree] = boundary + end + # Check "left face" (x_min) + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] + # Left face is position 1 in p4est indexing + boundary_names[1, tree] = boundary + end + # Check "right face" (x_max) + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] + # Right face is position 2 in p4est indexing + boundary_names[2, tree] = boundary + end + # Check "bottom face" (z_min) + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] + # Bottom face is position 5 in p4est indexing + boundary_names[5, tree] = boundary + end + # Check "top face" (z_max) + if tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] + # Top face is position 6 in p4est indexing + boundary_names[6, tree] = boundary + end + end + end + + return boundary_names +end + """ P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; polydeg, RealT=Float64, From fb69699da0e79e2af7fae3bd2fc06c23ae18f3be Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 11:03:21 +0100 Subject: [PATCH 25/69] comment --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 0d7cecfb569..56b8a83a97c 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -479,7 +479,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # There's no simple and generic way to distinguish boundaries without any information given. # Name all of them :all. boundary_names = fill(:all, 2 * n_dimensions, n_trees) - else # Boundary information supplied + else # Boundary information given # Read in nodes belonging to boundaries node_set_dict = parse_node_sets(meshfile, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries From 24977275e691c0793484d70948e7f44dbacfa403 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 11:17:57 +0100 Subject: [PATCH 26/69] two boundaries --- .../elixir_euler_free_stream_boundaries.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl index 20d374da0f3..bdc4da26c1f 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl @@ -16,17 +16,18 @@ solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) ############################################################################### # Get the uncurved mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "cube_boundaries.inp") +default_mesh_file = joinpath(@__DIR__, "mesh_cube_with_boundaries.inp") isfile(default_mesh_file) || - download("https://gist.github.com/DanielDoehring/eefe73ae5d113bc3944a518b6e88e663/raw/359a58a808790f3c3efc050273270eb1cc8ee353/cube_boundaries.inp", + download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp", default_mesh_file) mesh_file = default_mesh_file -boundary_symbols = [:PhysicalSurface1] +boundary_symbols = [:PhysicalSurface1, :PhysicalSurface2] mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) -boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition)) +boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition), + :PhysicalSurface2 => BoundaryConditionDirichlet(initial_condition)) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) From a006756c305b0a46352b666fe483d107b342b7da Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 11:21:17 +0100 Subject: [PATCH 27/69] comment --- src/meshes/p4est_mesh.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 56b8a83a97c..6e6aa6d37c8 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -567,6 +567,8 @@ function parse_node_sets(meshfile, boundary_symbols) return nodes_dict end +# This function assigns the edges of elements to boundaries by +# checking if the nodes of the edges are part of nodesets which correspond to boundaries. function assign_boundaries_standard_abaqus!(boundary_names, n_trees, element_node_matrix, node_set_dict, ::Val{2}) # 2D version @@ -606,6 +608,8 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, return boundary_names end +# This function assigns the edges of elements to boundaries by +# checking if the nodes of the edges are part of nodesets which correspond to boundaries. function assign_boundaries_standard_abaqus!(boundary_names, n_trees, element_node_matrix, node_set_dict, ::Val{3}) # 3D version From 43c4f4313bd3b0133d705c91d9e0df658dd2504f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 11:21:46 +0100 Subject: [PATCH 28/69] comment --- src/meshes/p4est_mesh.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 6e6aa6d37c8..9694cf1c161 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -577,7 +577,7 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, # For node labeling, see # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 # and search for "Node ordering and face numbering on elements" - for boundary in keys(node_set_dict) + for boundary in keys(node_set_dict) # Loop over specified boundaries # Check bottom edge if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] @@ -617,7 +617,7 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, tree_nodes = element_node_matrix[tree, :] # For node labeling, see # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html - for boundary in keys(node_set_dict) + for boundary in keys(node_set_dict) # Loop over specified boundaries # Check "front face" (y_min) if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && From aca432d1f4b9a00b11f9a251f301e9a66efc87f1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 15:04:44 +0100 Subject: [PATCH 29/69] rename elixir, change mesh --- ...h2.jl => elixir_euler_NACA6412airfoil_mach2.jl} | 14 +++++++------- test/test_p4est_2d.jl | 12 ++++++------ 2 files changed, 13 insertions(+), 13 deletions(-) rename examples/p4est_2d_dgsem/{elixir_euler_airfoil_mach2.jl => elixir_euler_NACA6412airfoil_mach2.jl} (84%) diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl similarity index 84% rename from examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl rename to examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl index 88799720029..a27836552bd 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl @@ -66,19 +66,19 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) -mesh_file = joinpath(@__DIR__, "NACA6412.inp") +mesh_file = joinpath(@__DIR__, "mesh_NACA6412.inp") isfile(mesh_file) || - download("https://gist.github.com/DanielDoehring/b434005800a1b0c0b4b50a8772283019/raw/1f6649b3370163d75d268176b51775f8685dd81d/NACA6412.inp", + download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp", mesh_file) -boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] +boundary_symbols = [:PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4] mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) -boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary - :PhysicalLine20 => boundary_condition_supersonic_outflow, # Right boundary - :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil - :PhysicalLine40 => boundary_condition_supersonic_outflow) # Top and bottom boundary +boundary_conditions = Dict(:PhysicalLine1 => boundary_condition_supersonic_inflow, # Left boundary + :PhysicalLine2 => boundary_condition_supersonic_outflow, # Right boundary + :PhysicalLine3 => boundary_condition_slip_wall, # Airfoil + :PhysicalLine4 => boundary_condition_supersonic_outflow) # Top and bottom boundary semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index d921090be77..b034091a175 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -391,15 +391,15 @@ end end end -@trixi_testset "elixir_euler_airfoil_mach2.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_airfoil_mach2.jl"), +@trixi_testset "elixir_euler_NACA6412airfoil_mach2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA6412airfoil_mach2.jl"), l2=[ - 1.682112049215034e-10, 3.439606229463632e-10, - 1.973815936404324e-10, 7.666245265417042e-10, + 1.9752162683735258e-9, 3.150450205812513e-9, + 1.8885499402935914e-9, 7.273629602920966e-9, ], linf=[ - 2.636251394960709e-8, 4.4272578936244145e-8, - 2.8049958258193942e-8, 1.0125492977408612e-7, + 6.007577890709825e-7, 1.005273289944597e-6, + 5.948514542597182e-7, 2.3111764217986774e-6, ], tspan=(0.0, 0.1)) # Ensure that we do not have excessive memory allocations From 0ffaba125d1851b86a65ff820121f8e5a2a1b70c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 7 Jan 2024 15:08:09 +0100 Subject: [PATCH 30/69] comments --- examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl index a27836552bd..6673053d88f 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl @@ -66,6 +66,8 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) +# Mesh generated from the following gmsh geometry input file: +# https://gist.githubusercontent.com/DanielDoehring/5ade6d93629f0d8c23a598812dbee2a9/raw/d2bc904fe92146eae1a36156e7f5c535dc1a80f1/NACA6412.geo mesh_file = joinpath(@__DIR__, "mesh_NACA6412.inp") isfile(mesh_file) || download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp", From f306b8681b77edd055d63f19334069b7d1fb9224 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 9 Jan 2024 15:40:15 +0100 Subject: [PATCH 31/69] add doc --- docs/src/meshes/p4est_mesh.md | 254 +++++++++++++++++++++++++++++++++- 1 file changed, 251 insertions(+), 3 deletions(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 1e5d782ebb6..30537a6b35f 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -55,7 +55,7 @@ This heading is used to indicate to the mesh constructor which of the above mapp create a curvilinear mesh. If the Abaqus file header is **not** present then the `P4estMesh` is created with the first strategy above. -#### List of corner nodes +#### List of corner nodes Next, prefaced with `*NODE`, comes a list of the physical `(x,y,z)` coordinates of all the corners. The first integer in the list of the corners provides its id number. @@ -71,7 +71,7 @@ Thus, for the two-dimensional example mesh this block of corner information is 7, 3.0, -1.0, 0.0 ``` -#### List of elements +#### List of elements The element connectivity is given after the list of corners. The header for this information block is ``` @@ -98,7 +98,9 @@ The construction of the element neighbor ids and identifying physical boundary s directly from the [`p4est`](https://github.com/cburstedde/p4est) library. For example, the neighbor connectivity is created in the mesh constructor using the wrapper `read_inp_p4est` function. -#### HOHQMesh boundary information +#### Encoding of boundaries + +##### HOHQMesh boundary information If present, any additional information in the mesh file that was created by `HOHQMesh` is prefaced with `** ` to make it an Abaqus comment. @@ -230,8 +232,36 @@ For completeness, we provide the entire Abaqus mesh file for the example mesh in ** Bottom --- Right --- ``` +##### Standard Abaqus format boundary information + +As an alternative to an Abaqus mesh generated by `HOHQMesh`, `.inp` files with boundary information encoded as nodesets `*NSET,NSET=` can be used to construct a `p4est` mesh. +This is especially useful for usage of existing meshes (consisting of bilinear elements) which could stem from the popular [`gmsh`](https://gmsh.info/) meshing software. + +In addition to the list of [nodes](#nodes) and [elements](#elements) given above, there are nodesets of the form +``` +*NSET,NSET=PhysicalLine1 +1, 4, 52, 53, 54, 55, 56, 57, 58, +``` +required which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. +By looping over every element and its associated edges consisting of two nodes, we query the read in `NSET`s if the current node pair is present. + +To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user has to supply the `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: + +```julia +boundary_symbols = [:PhysicalLine1] + +mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) +``` +By doing so, only nodesets with label in `boundary_symbols` are treated as boundaries and other nodesets which could be used for diagnostics are not treated as boundary. +Note that there is a leading double-colon `:` compared to the label in the `.inp` mesh file. +This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). + +A 2D example for this mesh generation is presented in `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`. + ### Mesh in three spatial dimensions +#### `HOHQMesh`-Extended Abaqus format + The 3D Abaqus file format with high-order boundary information from `HOHQMesh` is very similar to the 2D version discussed above. There are only three changes: @@ -346,4 +376,222 @@ transfinite map of the straight sided hexahedral element to find \mathbf{X}(\boldsymbol{\xi}) = \boldsymbol\Sigma(\boldsymbol{\xi}) - \mathcal{C}_{\texttt{edge}}(\boldsymbol{\xi}) + \mathbf{X}_{linear}(\boldsymbol{\xi}). +``` + +#### Construction from standard Abaqus + +Also for a mesh in standard Abaqus format there are no qualitative changes when going from 2D to 3D. +The most notable difference is that boundaries are formed in 3D by faces defined by four nodes while in 2D boundaries are edges consisting of two elements. +A simple mesh file used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl` is given below: +``` +*Heading + Something DIFFERENT FROM ** ***** HOHQMesh boundary information ***** ** +*NODE +1, -2, 0, 0 +2, -1, 0, 0 +3, -1, 1, 0 +4, -2, 1, 0 +5, -2, 0, 1 +6, -1, 0, 1 +7, -1, 1, 1 +8, -2, 1, 1 +9, -1.75, 1, 0 +10, -1.5, 1, 0 +11, -1.25, 1, 0 +12, -1, 0.75000000000035, 0 +13, -1, 0.50000000000206, 0 +14, -1, 0.25000000000104, 0 +15, -1.25, 0, 0 +16, -1.5, 0, 0 +17, -1.75, 0, 0 +18, -2, 0.24999999999941, 0 +19, -2, 0.49999999999869, 0 +20, -2, 0.74999999999934, 0 +21, -1.75, 0, 1 +22, -1.5, 0, 1 +23, -1.25, 0, 1 +24, -1, 0.24999999999941, 1 +25, -1, 0.49999999999869, 1 +26, -1, 0.74999999999934, 1 +27, -1.25, 1, 1 +28, -1.5, 1, 1 +29, -1.75, 1, 1 +30, -2, 0.75000000000035, 1 +31, -2, 0.50000000000206, 1 +32, -2, 0.25000000000104, 1 +33, -2, 0, 0.24999999999941 +34, -2, 0, 0.49999999999869 +35, -2, 0, 0.74999999999934 +36, -2, 1, 0.24999999999941 +37, -2, 1, 0.49999999999869 +38, -2, 1, 0.74999999999934 +39, -1, 0, 0.24999999999941 +40, -1, 0, 0.49999999999869 +41, -1, 0, 0.74999999999934 +42, -1, 1, 0.24999999999941 +43, -1, 1, 0.49999999999869 +44, -1, 1, 0.74999999999934 +45, -1.25, 0.25000000000063, 0 +46, -1.25, 0.50000000000122, 0 +47, -1.25, 0.7500000000001, 0 +48, -1.5, 0.25000000000023, 0 +49, -1.5, 0.50000000000038, 0 +50, -1.5, 0.74999999999984, 0 +51, -1.75, 0.24999999999982, 0 +52, -1.75, 0.49999999999953, 0 +53, -1.75, 0.74999999999959, 0 +54, -1.75, 0.25000000000063, 1 +55, -1.75, 0.50000000000122, 1 +56, -1.75, 0.7500000000001, 1 +57, -1.5, 0.25000000000023, 1 +58, -1.5, 0.50000000000038, 1 +59, -1.5, 0.74999999999984, 1 +60, -1.25, 0.24999999999982, 1 +61, -1.25, 0.49999999999953, 1 +62, -1.25, 0.74999999999959, 1 +63, -2, 0.24999999999982, 0.24999999999941 +64, -2, 0.49999999999953, 0.24999999999941 +65, -2, 0.74999999999959, 0.24999999999941 +66, -2, 0.25000000000023, 0.49999999999869 +67, -2, 0.50000000000038, 0.49999999999869 +68, -2, 0.74999999999984, 0.49999999999869 +69, -2, 0.25000000000063, 0.74999999999934 +70, -2, 0.50000000000122, 0.74999999999934 +71, -2, 0.7500000000001, 0.74999999999934 +72, -1.25, 1, 0.74999999999934 +73, -1.25, 1, 0.49999999999869 +74, -1.25, 1, 0.24999999999941 +75, -1.5, 1, 0.74999999999934 +76, -1.5, 1, 0.49999999999869 +77, -1.5, 1, 0.24999999999941 +78, -1.75, 1, 0.74999999999934 +79, -1.75, 1, 0.49999999999869 +80, -1.75, 1, 0.24999999999941 +81, -1, 0.25000000000063, 0.24999999999941 +82, -1, 0.50000000000122, 0.24999999999941 +83, -1, 0.7500000000001, 0.24999999999941 +84, -1, 0.25000000000023, 0.49999999999869 +85, -1, 0.50000000000038, 0.49999999999869 +86, -1, 0.74999999999984, 0.49999999999869 +87, -1, 0.24999999999982, 0.74999999999934 +88, -1, 0.49999999999953, 0.74999999999934 +89, -1, 0.74999999999959, 0.74999999999934 +90, -1.75, 0, 0.74999999999934 +91, -1.75, 0, 0.49999999999869 +92, -1.75, 0, 0.24999999999941 +93, -1.5, 0, 0.74999999999934 +94, -1.5, 0, 0.49999999999869 +95, -1.5, 0, 0.24999999999941 +96, -1.25, 0, 0.74999999999934 +97, -1.25, 0, 0.49999999999869 +98, -1.25, 0, 0.24999999999941 +99, -1.75, 0.25000000000043, 0.74999999999934 +100, -1.75, 0.25000000000023, 0.49999999999869 +101, -1.75, 0.25000000000002, 0.24999999999941 +102, -1.75, 0.5000000000008, 0.74999999999934 +103, -1.75, 0.50000000000038, 0.49999999999869 +104, -1.75, 0.49999999999995, 0.24999999999941 +105, -1.75, 0.74999999999997, 0.74999999999934 +106, -1.75, 0.74999999999984, 0.49999999999869 +107, -1.75, 0.74999999999972, 0.24999999999941 +108, -1.5, 0.25000000000023, 0.74999999999934 +109, -1.5, 0.25000000000023, 0.49999999999869 +110, -1.5, 0.25000000000023, 0.24999999999941 +111, -1.5, 0.50000000000038, 0.74999999999934 +112, -1.5, 0.50000000000038, 0.49999999999869 +113, -1.5, 0.50000000000038, 0.24999999999941 +114, -1.5, 0.74999999999984, 0.74999999999934 +115, -1.5, 0.74999999999984, 0.49999999999869 +116, -1.5, 0.74999999999984, 0.24999999999941 +117, -1.25, 0.25000000000002, 0.74999999999934 +118, -1.25, 0.25000000000023, 0.49999999999869 +119, -1.25, 0.25000000000043, 0.24999999999941 +120, -1.25, 0.49999999999995, 0.74999999999934 +121, -1.25, 0.50000000000038, 0.49999999999869 +122, -1.25, 0.5000000000008, 0.24999999999941 +123, -1.25, 0.74999999999972, 0.74999999999934 +124, -1.25, 0.74999999999984, 0.49999999999869 +125, -1.25, 0.74999999999997, 0.24999999999941 +******* E L E M E N T S ************* +*ELEMENT, type=C3D8, ELSET=Volume1 +153, 54, 21, 5, 32, 99, 90, 35, 69 +154, 99, 90, 35, 69, 100, 91, 34, 66 +155, 100, 91, 34, 66, 101, 92, 33, 63 +156, 101, 92, 33, 63, 51, 17, 1, 18 +157, 55, 54, 32, 31, 102, 99, 69, 70 +158, 102, 99, 69, 70, 103, 100, 66, 67 +159, 103, 100, 66, 67, 104, 101, 63, 64 +160, 104, 101, 63, 64, 52, 51, 18, 19 +161, 56, 55, 31, 30, 105, 102, 70, 71 +162, 105, 102, 70, 71, 106, 103, 67, 68 +163, 106, 103, 67, 68, 107, 104, 64, 65 +164, 107, 104, 64, 65, 53, 52, 19, 20 +165, 29, 56, 30, 8, 78, 105, 71, 38 +166, 78, 105, 71, 38, 79, 106, 68, 37 +167, 79, 106, 68, 37, 80, 107, 65, 36 +168, 80, 107, 65, 36, 9, 53, 20, 4 +169, 57, 22, 21, 54, 108, 93, 90, 99 +170, 108, 93, 90, 99, 109, 94, 91, 100 +171, 109, 94, 91, 100, 110, 95, 92, 101 +172, 110, 95, 92, 101, 48, 16, 17, 51 +173, 58, 57, 54, 55, 111, 108, 99, 102 +174, 111, 108, 99, 102, 112, 109, 100, 103 +175, 112, 109, 100, 103, 113, 110, 101, 104 +176, 113, 110, 101, 104, 49, 48, 51, 52 +177, 59, 58, 55, 56, 114, 111, 102, 105 +178, 114, 111, 102, 105, 115, 112, 103, 106 +179, 115, 112, 103, 106, 116, 113, 104, 107 +180, 116, 113, 104, 107, 50, 49, 52, 53 +181, 28, 59, 56, 29, 75, 114, 105, 78 +182, 75, 114, 105, 78, 76, 115, 106, 79 +183, 76, 115, 106, 79, 77, 116, 107, 80 +184, 77, 116, 107, 80, 10, 50, 53, 9 +185, 60, 23, 22, 57, 117, 96, 93, 108 +186, 117, 96, 93, 108, 118, 97, 94, 109 +187, 118, 97, 94, 109, 119, 98, 95, 110 +188, 119, 98, 95, 110, 45, 15, 16, 48 +189, 61, 60, 57, 58, 120, 117, 108, 111 +190, 120, 117, 108, 111, 121, 118, 109, 112 +191, 121, 118, 109, 112, 122, 119, 110, 113 +192, 122, 119, 110, 113, 46, 45, 48, 49 +193, 62, 61, 58, 59, 123, 120, 111, 114 +194, 123, 120, 111, 114, 124, 121, 112, 115 +195, 124, 121, 112, 115, 125, 122, 113, 116 +196, 125, 122, 113, 116, 47, 46, 49, 50 +197, 27, 62, 59, 28, 72, 123, 114, 75 +198, 72, 123, 114, 75, 73, 124, 115, 76 +199, 73, 124, 115, 76, 74, 125, 116, 77 +200, 74, 125, 116, 77, 11, 47, 50, 10 +201, 24, 6, 23, 60, 87, 41, 96, 117 +202, 87, 41, 96, 117, 84, 40, 97, 118 +203, 84, 40, 97, 118, 81, 39, 98, 119 +204, 81, 39, 98, 119, 14, 2, 15, 45 +205, 25, 24, 60, 61, 88, 87, 117, 120 +206, 88, 87, 117, 120, 85, 84, 118, 121 +207, 85, 84, 118, 121, 82, 81, 119, 122 +208, 82, 81, 119, 122, 13, 14, 45, 46 +209, 26, 25, 61, 62, 89, 88, 120, 123 +210, 89, 88, 120, 123, 86, 85, 121, 124 +211, 86, 85, 121, 124, 83, 82, 122, 125 +212, 83, 82, 122, 125, 12, 13, 46, 47 +213, 7, 26, 62, 27, 44, 89, 123, 72 +214, 44, 89, 123, 72, 43, 86, 124, 73 +215, 43, 86, 124, 73, 42, 83, 125, 74 +216, 42, 83, 125, 74, 3, 12, 47, 11 +*NSET,NSET=PhysicalSurface1 +1, 2, 3, 4, 5, 6, 7, 8, 9, 10, +11, 12, 13, 14, 15, 16, 17, 18, 19, 20, +21, 22, 23, 24, 25, 26, 27, 28, 29, 30, +31, 32, 33, 34, 35, 36, 37, 38, 45, 46, +47, 48, 49, 50, 51, 52, 53, 54, 55, 56, +57, 58, 59, 60, 61, 62, 63, 64, 65, 66, +67, 68, 69, 70, 71, +*NSET,NSET=PhysicalSurface2 +1, 2, 3, 4, 5, 6, 7, 8, 9, 10, +11, 12, 13, 14, 15, 16, 17, 21, 22, 23, +24, 25, 26, 27, 28, 29, 33, 34, 35, 36, +37, 38, 39, 40, 41, 42, 43, 44, 72, 73, +74, 75, 76, 77, 78, 79, 80, 81, 82, 83, +84, 85, 86, 87, 88, 89, 90, 91, 92, 93, +94, 95, 96, 97, 98, ``` \ No newline at end of file From 79e110748307622bf593fb0ec943022ab969cec0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 10 Jan 2024 11:03:09 +0100 Subject: [PATCH 32/69] avoid unicode --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 9694cf1c161..25e4dbef566 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -498,7 +498,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, end function parse_elements(meshfile, n_trees, n_dims) - @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" + @assert n_dims in [2, 3] "Only 2D and 3D meshes are supported" # Valid element types (that can be processed by p4est) based on dimension element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", From f00ac143f903db5bce841022800e7de87a1a80ed Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 10 Jan 2024 11:19:06 +0100 Subject: [PATCH 33/69] improve doc --- docs/src/meshes/p4est_mesh.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 30537a6b35f..2a48da2a5b9 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -256,7 +256,8 @@ By doing so, only nodesets with label in `boundary_symbols` are treated as bound Note that there is a leading double-colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). -A 2D example for this mesh generation is presented in `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`. +A 2D example for this mesh generation is presented in +`examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`. ### Mesh in three spatial dimensions From 682adc34406846ebbe1ce4e43af93fb08b5cce52 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 10 Jan 2024 11:35:18 +0100 Subject: [PATCH 34/69] Add tutorial --- docs/literate/src/files/p4est_from_gmsh.jl | 452 +++++++++++++++++++++ docs/make.jl | 1 + 2 files changed, 453 insertions(+) create mode 100644 docs/literate/src/files/p4est_from_gmsh.jl diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl new file mode 100644 index 00000000000..f7925e68026 --- /dev/null +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -0,0 +1,452 @@ +#src # `P4estMesh` from [`gmsh`](https://gmsh.info/) + +# Trixi.jl supports numerical approximations from structured and unstructured quadrilateral meshes +# with the [`P4estMesh`](@ref) mesh type. + +# The purpose of this tutorial is to demonstrate how to use the `P4estMesh` +# functionality of Trixi.jl for existing meshes with straight-sided (bilinear) elements/cells. +# This begins by running and visualizing an available unstructured quadrilateral mesh example. +# Then, the tutorial will cover how to use existing meshes generated by [`gmsh`](https://gmsh.info/) +# or any other meshing software that can export to the Abaqus input `.inp` format. + +# ## Running the simulation of a near-field flow around an airfoil + +# Trixi.jl supports solving hyperbolic problems on several mesh types. +# A somewhat complex example that employs the `P4estMesh` is the near-field simulation of a +# Mach 2 flow around the NACA6412 airfoil. + +using Trixi +redirect_stdio(stdout=devnull, stderr=devnull) do # code that prints annoying stuff we don't want to see here #hide #md +trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_euler_NACA6412airfoil_mach2.jl"), tspan=(0.0, 0.5)) +end #hide #md + +# Conveniently, we can use the Plots package to have a first look at the results: +# ```julia +# using Plots +# pd = PlotData2D(sol) +# plot(pd["rho"]) +# plot!(getmesh(pd)) +# ``` + +# ## Creating a mesh using `gmsh` + +# The creation of an unstructured quadrilateral mesh using `gmsh` is driven by a **geometry file**. +# There are plenty of possibilities for the user, see the [documentation](https://gmsh.info/doc/texinfo/gmsh.html) and [tutorials](https://gitlab.onelab.info/gmsh/gmsh/tree/master/tutorials). + +# To begin, we provide a complete geometry file in this tutorial. After this we give a breakdown +# of the most important parts required for succesfull mesh generation that can later be used by the `p4est` library +# and `Trixi.jl`. + +# The associated `NACA6412.geo` file is given below: +# ```c++ +# // GMSH geometry script for a NACA 6412 airfoil in a box +# // see https://github.com/cfsengineering/GMSH-Airfoil-2D +# +# // outer bounding box +# Point(1) = {-1.25, -0.5, 0, 1.0}; +# Point(2) = {1.25, -0.5, 0, 1.0}; +# Point(3) = {1.25, 0.5, 0, 1.0}; +# Point(4) = {-1.25, 0.5, 0, 1.0}; +# +# // lines of the bounding box +# Line(1) = {1, 2}; +# Line(2) = {2, 3}; +# Line(3) = {3, 4}; +# Line(4) = {4, 1}; +# // outer box +# Line Loop(8) = {1, 2, 3, 4}; +# +# // Settings +# // This value gives the global element size factor (lower -> finer mesh) +# Mesh.CharacteristicLengthFactor = 1.0 * 2^(-3); +# // Insist on quads instead of default triangles +# Mesh.RecombineAll = 1; +# // Violet instead of green base color for better visibility +# Mesh.ColorCarousel = 0; +# +# // points of the airfoil contour +# Point(5) = {-0.4900332889206208, 0.09933466539753061, 0, 0.125}; +# Point(6) = {-0.4900274857651495, 0.1021542752054094, 0, 0.125}; +# Point(7) = {-0.4894921489729144, 0.1049830248247787, 0, 0.125}; +# Point(8) = {-0.4884253336670712, 0.1078191282319664, 0, 0.125}; +# Point(9) = {-0.4868257975566199, 0.1106599068424483, 0, 0.125}; +# Point(10) = {-0.4846930063965668, 0.1135018003016681, 0, 0.125}; +# Point(11) = {-0.4820271400142729, 0.1163403835785654, 0, 0.125}; +# Point(12) = {-0.4788290988083472, 0.1191703902233889, 0, 0.125}; +# Point(13) = {-0.4751005105908123, 0.1219857416089041, 0, 0.125}; +# Point(14) = {-0.4708437376101668, 0.1247795819332056, 0, 0.125}; +# Point(15) = {-0.4660618835629463, 0.1275443187232316, 0, 0.125}; +# Point(16) = {-0.4607588003749649, 0.1302716685409717, 0, 0.125}; +# Point(17) = {-0.4549390945110529, 0.132952707559475, 0, 0.125}; +# Point(18) = {-0.448608132554204, 0.1355779266432996, 0, 0.125}; +# Point(19) = {-0.4417720457819508, 0.138137290538182, 0, 0.125}; +# Point(20) = {-0.4344377334597768, 0.140620300747629, 0, 0.125}; +# Point(21) = {-0.4266128645686593, 0.1430160616500159, 0, 0.125}; +# Point(22) = {-0.4183058776865576, 0.1453133493887722, 0, 0.125}; +# Point(23) = {-0.4095259787518715, 0.147500683050503, 0, 0.125}; +# Point(24) = {-0.4002831364505879, 0.1495663976315875, 0, 0.125}; +# Point(25) = {-0.3905880749878933, 0.1514987182830453, 0, 0.125}; +# Point(26) = {-0.3804522640292948, 0.1532858353164163, 0, 0.125}; +# Point(27) = {-0.3698879056254708, 0.1549159794501833, 0, 0.125}; +# Point(28) = {-0.3589079179688306, 0.1563774967770029, 0, 0.125}; +# Point(29) = {-0.3475259158676376, 0.1576589229368209, 0, 0.125}; +# Point(30) = {-0.3357561878650377, 0.158749055989923, 0, 0.125}; +# Point(31) = {-0.3236136699747923, 0.1596370274972017, 0, 0.125}; +# Point(32) = {-0.3111139160522804, 0.1603123713324616, 0, 0.125}; +# Point(33) = {-0.298273064867608, 0.160765089773461, 0, 0.125}; +# Point(34) = {-0.2851078039966239, 0.1609857164445887, 0, 0.125}; +# Point(35) = {-0.2716353306943914, 0.160965375714529, 0, 0.125}; +# Point(36) = {-0.2578733099632437, 0.1606958381868515, 0, 0.125}; +# Point(37) = {-0.2438398300730194, 0.1601695719599709, 0, 0.125}; +# Point(38) = {-0.2295533558334121, 0.1593797893750759, 0, 0.125}; +# Point(39) = {-0.2150326799566391, 0.1583204890160489, 0, 0.125}; +# Point(40) = {-0.2002968728818922, 0.1569864927736143, 0, 0.125}; +# Point(41) = {-0.18536523146042, 0.1553734778363979, 0, 0.125}; +# Point(42) = {-0.1702572269208345, 0.1534780035235666, 0, 0.125}; +# Point(43) = {-0.1549924525477129, 0.1512975329264932, 0, 0.125}; +# Point(44) = {-0.1395905715122586, 0.1488304493795921, 0, 0.125}; +# Point(45) = {-0.1240712652914332, 0.1460760678321895, 0, 0.125}; +# Point(46) = {-0.1084541831014299, 0.1430346412430583, 0, 0.125}; +# Point(47) = {-0.09275889275279087, 0.1397073621660917, 0, 0.125}; +# Point(48) = {-0.07700483330818747, 0.1360963597385416, 0, 0.125}; +# Point(49) = {-0.06151286635366404, 0.1323050298149023, 0, 0.125}; +# Point(50) = {-0.04602933219022032, 0.1283521764905442, 0, 0.125}; +# Point(51) = {-0.03051345534800332, 0.1242331665904082, 0, 0.125}; +# Point(52) = {-0.01498163190522334, 0.1199540932779839, 0, 0.125}; +# Point(53) = {0.0005498526140696458, 0.1155214539466913, 0, 0.125}; +# Point(54) = {0.01606484191716884, 0.1109421303284033, 0, 0.125}; +# Point(55) = {0.03154732664394777, 0.106223368423828, 0, 0.125}; +# Point(56) = {0.0469814611314705, 0.1013727584299359, 0, 0.125}; +# Point(57) = {0.06235157928986135, 0.09639821481480275, 0, 0.125}; +# Point(58) = {0.07764220964363855, 0.09130795666388933, 0, 0.125}; +# Point(59) = {0.09283808959671735, 0.08611048839446452, 0, 0.125}; +# Point(60) = {0.1079241789809607, 0.08081458090718853, 0, 0.125}; +# Point(61) = {0.1228856729475325, 0.07542925321638272, 0, 0.125}; +# Point(62) = {0.1377080142575372, 0.06996375457378261, 0, 0.125}; +# Point(63) = {0.1523769050236616, 0.06442754707512513, 0, 0.125}; +# Point(64) = {0.1668783179480157, 0.05883028871526293, 0, 0.125}; +# Point(65) = {0.1811985070933818, 0.05318181683604975, 0, 0.125}; +# Point(66) = {0.1953240182159306, 0.04749213189240609, 0, 0.125}; +# Point(67) = {0.2092416986775084, 0.04177138144606024, 0, 0.125}; +# Point(68) = {0.2229387069452062, 0.03602984428372727, 0, 0.125}; +# Point(69) = {0.2364025216754475, 0.03027791454712048, 0, 0.125}; +# Point(70) = {0.2496209503696738, 0.02452608575629232, 0, 0.125}; +# Point(71) = {0.2625821375791982, 0.01878493460541621, 0, 0.125}; +# Point(72) = {0.2752745726282818, 0.01306510441121807, 0, 0.125}; +# Point(73) = {0.28768709681727, 0.007377288098728577, 0, 0.125}; +# Point(74) = {0.2998089100619555, 0.001732210616722449, 0, 0.125}; +# Point(75) = {0.3116295769214332, -0.003859389314124759, 0, 0.125}; +# Point(76) = {0.3231390319647309, -0.009386778203927332, 0, 0.125}; +# Point(77) = {0.3343275844265582, -0.01483924761490708, 0, 0.125}; +# Point(78) = {0.3451859221046181, -0.02020613485126957, 0, 0.125}; +# Point(79) = {0.3557051144551212, -0.02547684454806881, 0, 0.125}; +# Point(80) = {0.3658766148492779, -0.03064087116872238, 0, 0.125}; +# Point(81) = {0.3756922619615632, -0.0356878223992288, 0, 0.125}; +# Point(82) = {0.3851442802702071, -0.0406074434050937, 0, 0.125}; +# Point(83) = {0.394225279661484, -0.04538964189492445, 0, 0.125}; +# Point(84) = {0.4029282541416501, -0.05002451391298904, 0, 0.125}; +# Point(85) = {0.4112465796735204, -0.05450237026215737, 0, 0.125}; +# Point(86) = {0.4191740111683733, -0.05881376343890812, 0, 0.125}; +# Point(87) = {0.4267046786777481, -0.06294951494382847, 0, 0.125}; +# Point(88) = {0.4338330828434404, -0.06690074281456823, 0, 0.125}; +# Point(89) = {0.4405540896772232, -0.07065888921378868, 0, 0.125}; +# Point(90) = {0.4468629247542237, -0.07421574789251445, 0, 0.125}; +# Point(91) = {0.4527551669150955, -0.0775634913396257, 0, 0.125}; +# Point(92) = {0.4582267415819197, -0.08069469742118066, 0, 0.125}; +# Point(93) = {0.4632739138007936, -0.08360237530891265, 0, 0.125}; +# Point(94) = {0.4678932811302005, -0.08627999049569551, 0, 0.125}; +# Point(95) = {0.4720817664982195, -0.08872148869699745, 0, 0.125}; +# Point(96) = {0.4758366111533843, -0.09092131844134463, 0, 0.125}; +# Point(97) = {0.4791553678333992, -0.09287445215953141, 0, 0.125}; +# Point(98) = {0.4820358942729613, -0.09457640559161551, 0, 0.125}; +# Point(99) = {0.4844763471666588, -0.09602325534252773, 0, 0.125}; +# Point(100) = {0.4864751766953637, -0.09721165443119822, 0, 0.125}; +# Point(101) = {0.4880311217148797, -0.09813884569428721, 0, 0.125}; +# Point(102) = {0.4891432056939881, -0.09880267292366274, 0, 0.125}; +# Point(103) = {0.4898107334756874, -0.09920158963645126, 0, 0.125}; +# Point(104) = {0.4900332889206208, -0.09933466539753058, 0, 0.125}; +# Point(105) = {0.4897824225031319, -0.09926905587549506, 0, 0.125}; +# Point(106) = {0.4890301110661922, -0.09907236506934192, 0, 0.125}; +# Point(107) = {0.4877772173496635, -0.09874500608402761, 0, 0.125}; +# Point(108) = {0.48602517690576, -0.09828766683852558, 0, 0.125}; +# Point(109) = {0.4837759946062035, -0.09770130916007558, 0, 0.125}; +# Point(110) = {0.4810322398085871, -0.09698716747297723, 0, 0.125}; +# Point(111) = {0.4777970402368822, -0.09614674703990023, 0, 0.125}; +# Point(112) = {0.4740740746447117, -0.09518182170326678, 0, 0.125}; +# Point(113) = {0.4698675643422793, -0.09409443106501386, 0, 0.125}; +# Point(114) = {0.4651822636784212, -0.09288687703518478, 0, 0.125}; +# Point(115) = {0.460023449577924, -0.09156171967354482, 0, 0.125}; +# Point(116) = {0.4543969102408585, -0.09012177224394632, 0, 0.125}; +# Point(117) = {0.4483089331151018, -0.08857009539864649, 0, 0.125}; +# Point(118) = {0.4417662922553667, -0.08690999040934186, 0, 0.125}; +# Point(119) = {0.4347762351819332, -0.0851449913634191, 0, 0.125}; +# Point(120) = {0.4273464693498908, -0.08327885624791403, 0, 0.125}; +# Point(121) = {0.419485148335155, -0.08131555684993674, 0, 0.125}; +# Point(122) = {0.411200857836944, -0.07925926741086739, 0, 0.125}; +# Point(123) = {0.4025026015879757, -0.07711435198240155, 0, 0.125}; +# Point(124) = {0.3933997872536054, -0.07488535044544484, 0, 0.125}; +# Point(125) = {0.3839022123897198, -0.07257696316779733, 0, 0.125}; +# Point(126) = {0.3740200505167618, -0.07019403429336624, 0, 0.125}; +# Point(127) = {0.3637638373540689, -0.06774153367408606, 0, 0.125}; +# Point(128) = {0.3531444572451353, -0.06522453747557577, 0, 0.125}; +# Point(129) = {0.3421731297908021, -0.06264820750853495, 0, 0.125}; +# Point(130) = {0.3308613966940724, -0.06001776935966011, 0, 0.125}; +# Point(131) = {0.3192211088076166, -0.05733848941811218, 0, 0.125}; +# Point(132) = {0.3072644133633567, -0.05461565091590426, 0, 0.125}; +# Point(133) = {0.2950037413531683, -0.05185452912263369, 0, 0.125}; +# Point(134) = {0.2824517950208982, -0.04906036585632723, 0, 0.125}; +# Point(135) = {0.2696215354188702, -0.04623834349241404, 0, 0.125}; +# Point(136) = {0.2565261699769623, -0.04339355867155523, 0, 0.125}; +# Point(137) = {0.2431791400293651, -0.04053099592384862, 0, 0.125}; +# Point(138) = {0.2295941082432855, -0.03765550144139543, 0, 0.125}; +# Point(139) = {0.2157849458952252, -0.03477175724299444, 0, 0.125}; +# Point(140) = {0.2017657199439165, -0.03188425598348005, 0, 0.125}; +# Point(141) = {0.187550679854507, -0.02899727666564914, 0, 0.125}; +# Point(142) = {0.1731542441359161, -0.02611486151457043, 0, 0.125}; +# Point(143) = {0.1585909865622793, -0.02324079427214604, 0, 0.125}; +# Point(144) = {0.1438756220597465, -0.02037858016395433, 0, 0.125}; +# Point(145) = {0.129022992251319, -0.0175314277805827, 0, 0.125}; +# Point(146) = {0.1140480506645569, -0.01470223310184333, 0, 0.125}; +# Point(147) = {0.09896584761949168, -0.01189356587453844, 0, 0.125}; +# Point(148) = {0.08379151482656089, -0.009107658532933174, 0, 0.125}; +# Point(149) = {0.06854024973648176, -0.006346397826038436, 0, 0.125}; +# Point(150) = {0.05322729969528361, -0.003611319287478529, 0, 0.125}; +# Point(151) = {0.03786794596792287, -0.00090360465249055, 0, 0.125}; +# Point(152) = {0.0224774877026287, 0.00177591770710904, 0, 0.125}; +# Point(153) = {0.007071225915134205, 0.004426769294862437, 0, 0.125}; +# Point(154) = {-0.00833555242305456, 0.007048814950562587, 0, 0.125}; +# Point(155) = {-0.02372759010533726, 0.009642253300220296, 0, 0.125}; +# Point(156) = {-0.03908967513210498, 0.01220760427359278, 0, 0.125}; +# Point(157) = {-0.05440665578848514, 0.01474569380579989, 0, 0.125}; +# Point(158) = {-0.06966345527617318, 0.01725763587663899, 0, 0.125}; +# Point(159) = {-0.08484508582421563, 0.01974481207672138, 0, 0.125}; +# Point(160) = {-0.09987987792382108, 0.02219618763023203, 0, 0.125}; +# Point(161) = {-0.1145078729404739, 0.02450371976411331, 0, 0.125}; +# Point(162) = {-0.1290321771824579, 0.0267015185742735, 0, 0.125}; +# Point(163) = {-0.143440065923266, 0.02879471001709845, 0, 0.125}; +# Point(164) = {-0.1577189448447794, 0.03078883518202784, 0, 0.125}; +# Point(165) = {-0.1718563428491159, 0.03268980457290044, 0, 0.125}; +# Point(166) = {-0.1858399037768357, 0.03450385196323842, 0, 0.125}; +# Point(167) = {-0.1996573773370766, 0.03623748825421298, 0, 0.125}; +# Point(168) = {-0.2132966095779342, 0.03789745574015834, 0, 0.125}; +# Point(169) = {-0.2267455332406906, 0.0394906831577609, 0, 0.125}; +# Point(170) = {-0.2399921583489679, 0.04102424186233269, 0, 0.125}; +# Point(171) = {-0.2530245633834605, 0.04250530343879837, 0, 0.125}; +# Point(172) = {-0.2658308873846617, 0.04394109901707172, 0, 0.125}; +# Point(173) = {-0.2783993233102972, 0.04533888052223981, 0, 0.125}; +# Point(174) = {-0.2907181129514687, 0.04670588405019788, 0, 0.125}; +# Point(175) = {-0.3027755436824813, 0.0480492955198111, 0, 0.125}; +# Point(176) = {-0.3145599472847223, 0.04937621871394801, 0, 0.125}; +# Point(177) = {-0.3260597010456697, 0.05069364578437131, 0, 0.125}; +# Point(178) = {-0.337263231291058, 0.05200843025992359, 0, 0.125}; +# Point(179) = {-0.3481590194623916, 0.05332726256406103, 0, 0.125}; +# Point(180) = {-0.3587356108043638, 0.05465664801682354, 0, 0.125}; +# Point(181) = {-0.3689816256782782, 0.0560028872679817, 0, 0.125}; +# Point(182) = {-0.3788857734692287, 0.05737205908247899, 0, 0.125}; +# Point(183) = {-0.3884368690074614, 0.05877000537646382, 0, 0.125}; +# Point(184) = {-0.3976238513788748, 0.06020231838219783, 0, 0.125}; +# Point(185) = {-0.40643580495675, 0.06167432980291591, 0, 0.125}; +# Point(186) = {-0.4148619824472646, 0.06319110180426264, 0, 0.125}; +# Point(187) = {-0.4228918297057104, 0.06475741967717524, 0, 0.125}; +# Point(188) = {-0.43051501204915, 0.06637778599795482, 0, 0.125}; +# Point(189) = {-0.4377214417649294, 0.06805641610468524, 0, 0.125}; +# Point(190) = {-0.4445013064933708, 0.06979723470503821, 0, 0.125}; +# Point(191) = {-0.4508450981473512, 0.07160387342876083, 0, 0.125}; +# Point(192) = {-0.4567436420215075, 0.073479669138689, 0, 0.125}; +# Point(193) = {-0.4621881257395756, 0.07542766281688272, 0, 0.125}; +# Point(194) = {-0.4671701276898881, 0.07745059884734995, 0, 0.125}; +# Point(195) = {-0.471681644606229, 0.07955092452372269, 0, 0.125}; +# Point(196) = {-0.4757151179639407, 0.0817307896190848, 0, 0.125}; +# Point(197) = {-0.4792634588791559, 0.0839920458658267, 0, 0.125}; +# Point(198) = {-0.4823200712220043, 0.08633624620581726, 0, 0.125}; +# Point(199) = {-0.4848788726822436, 0.08876464368523246, 0, 0.125}; +# Point(200) = {-0.4869343135575803, 0.09127818988394577, 0, 0.125}; +# Point(201) = {-0.4884813930704814, 0.09387753278635144, 0, 0.125}; +# Point(202) = {-0.4895156730580155, 0.09656301401871749, 0, 0.125}; +# +# // splines of the airfoil +# Spline(5) = {5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104}; +# Spline(6) = {104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,5}; +# +# // airfoil +# Line Loop(9) = {5, 6}; +# // complete domain +# Plane Surface(1) = {8, 9}; +# +# // labeling of the boundary parts +# Physical Line(1) = {4}; // inflow +# Physical Line(2) = {2}; // outflow +# Physical Line(3) = {1, 3}; // airfoil +# Physical Line(4) = {5, 6}; // upper/lower wall +# Physical Surface(1) = {10}; +# ``` +# From which we can construct a mesh like this: +# ![mesh_screenshot](https://github.com/trixi-framework/Trixi.jl/assets/75639095/67adfe3d-d403-4cd3-acaa-971a34df0709) +# +# The first four points define the bounding box = (near-field) domain: +# ```c++ +# // outer bounding box +# Point(1) = {-1.25, -0.5, 0, 1.0}; +# Point(2) = {1.25, -0.5, 0, 1.0}; +# Point(3) = {1.25, 0.5, 0, 1.0}; +# Point(4) = {-1.25, 0.5, 0, 1.0}; +# ``` +# which is constructed from connecting the points in lines: +# ```c++ +# // outer box +# Line(1) = {1, 2}; +# Line(2) = {2, 3}; +# Line(3) = {3, 4}; +# Line(4) = {4, 1}; +# // outer box +# Line Loop(8) = {1, 2, 3, 4}; +# ``` +# +# This is followed by a couple (in principle optional) settings where the most important one is +# ```c++ +# // Insist on quads instead of default triangles +# Mesh.RecombineAll = 1; +# ``` +# which forces `gmsh` to generate quadrilateral elements instead of the default triangles. +# This is strictly required to be able to use the mesh later with `p4est` with supports only straight-sided quads, +# i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. +# See for more details the (short) ]documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. +# In principle, it should also be possible to use the recombine function of `gmsh` to convert the traingles to quads, +# but this is observed to be less robust than enforcing quads from the beginning. +# +# Then the airfoil is defined by a set of points: +# ```c++ +# // points of the airfoil contour +# Point(5) = {-0.4900332889206208, 0.09933466539753061, 0, 0.125}; +# Point(6) = {-0.4900274857651495, 0.1021542752054094, 0, 0.125}; +# ... +# ``` +# which are connected by splines for the upper and lower part of the airfoil: +# ```c++ +# // splines of the airfoil +# Spline(5) = {5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, +# ... +# 96,97,98,99,100,101,102,103,104}; +# Spline(6) = {104,105,106,107,108,109,110,111,112,113,114,115, +# ... +# 200,201,202,5}; +# ``` +# which are then connected to form a single line loop for easy physical group assignment: +# ```c++ +# // airfoil +# Line Loop(9) = {5, 6}; +# ``` +# +# At the end of the file the physical groups are defined: +# ```c++ +# // labeling of the boundary parts +# Physical Line(1) = {4}; // inflow +# Physical Line(2) = {2}; // outflow +# Physical Line(3) = {1, 3}; // airfoil +# Physical Line(4) = {5, 6}; // upper/lower wall +# ``` +# which are crucial for the correct assignment of boundary conditions in `Trixi.jl`. +# +# After opening this file in `gmsh`, meshing the geometry and exporting to Abaqus `.inp` format, +# we can have a look at the input file: +# ``` +# *Heading +# +# *NODE +# 1, -1.25, -0.5, 0 +# 2, 1.25, -0.5, 0 +# 3, 1.25, 0.5, 0 +# 4, -1.25, 0.5, 0 +# ... +# ******* E L E M E N T S ************* +# *ELEMENT, type=T3D2, ELSET=Line1 +# 1, 1, 7 +# ... +# *ELEMENT, type=CPS4, ELSET=Surface10 +# 191, 272, 46, 263, 807 +# ... +# *NSET,NSET=PhysicalLine10 +# 1, 4, 52, 53, 54, 55, 56, 57, 58, +# *NSET,NSET=PhysicalLine20 +# 2, 3, 26, 27, 28, 29, 30, 31, 32, +# *NSET,NSET=PhysicalLine30 +# 1, 2, 3, 4, 7, 8, 9, 10, 11, 12, +# 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, +# 23, 24, 25, 33, 34, 35, 36, 37, 38, 39, +# 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, +# 50, 51, +# *NSET,NSET=PhysicalLine40 +# 5, 6, 59, 60, 61, 62, 63, 64, 65, 66, +# 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, +# 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, +# 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, +# 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, +# 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, +# 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, +# 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, +# 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, +# 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, +# 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, +# 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, +# 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, +# 187, 188, 189, 190, +# ``` +# +# First, the coordinates of the nodes are listed, followed by the elements. +# Note that `gmsh` exports also line elements of type `T3D2` which are ignored by `p4est`. +# The relevant elements in 2D which form the gridcells are of type `CPS4` which are defined by their four corner nodes. +# This is followed by the nodesets encoded via `*NSET` which are used to assign boundary conditions in `Trixi.jl`. +# `Trixi.jl` parses the `.inp` file and assigns the edges (in 2D, surfaces in 3D) of elements to the corresponding boundary condition based on +# the supplied `boundary_symbols` that have to be supplied to the `P4estMesh` constructor: +# ```julia +# # boundary symbols +# boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] +# mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) +# ``` +# The same boundary symbols have then also be supplied to the semidiscretization alongside the +# corresponding physical boundary conditions: +# ```julia +# # Supersonic inflow boundary condition. +# # Calculate the boundary flux entirely from the external solution state, i.e., set +# # external solution state values for everything entering the domain. +# @inline function boundary_condition_supersonic_inflow(u_inner, +# normal_direction::AbstractVector, +# x, t, surface_flux_function, +# equations::CompressibleEulerEquations2D) +# u_boundary = initial_condition_mach2_flow(x, t, equations) +# flux = Trixi.flux(u_boundary, normal_direction, equations) +# +# return flux +# end +# +# # Supersonic outflow boundary condition. +# # Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow +# # except all the solution state values are set from the internal solution as everything leaves the domain +# @inline function boundary_condition_supersonic_outflow(u_inner, +# normal_direction::AbstractVector, x, +# t, +# surface_flux_function, +# equations::CompressibleEulerEquations2D) +# flux = Trixi.flux(u_inner, normal_direction, equations) +# +# boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary +# :PhysicalLine20 => boundary_condition_supersonic_outflow, # Right boundary +# :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil +# :PhysicalLine40 => boundary_condition_supersonic_outflow) # Top and bottom boundary +# +# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, +# boundary_conditions = boundary_conditions) +# ``` +# Note that you **have to** supply the `boundary_symbols` keyword to the `P4estMesh` constructor +# to select the boundaries from the available nodesets in the `.inp` file. + +# ## Package versions + +# These results were obtained using the following versions. + +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["Trixi", "OrdinaryDiffEq", "Plots", "Download"], + mode=PKGMODE_MANIFEST) diff --git a/docs/make.jl b/docs/make.jl index df8ac04be12..dee87371bd1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -65,6 +65,7 @@ files = [ "Adaptive mesh refinement" => "adaptive_mesh_refinement.jl", "Structured mesh with curvilinear mapping" => "structured_mesh_mapping.jl", "Unstructured meshes with HOHQMesh.jl" => "hohqmesh_tutorial.jl", + "P4est mesh from gmsh" => "p4est_from_gmsh.jl", # Topic: other stuff "Explicit time stepping" => "time_stepping.jl", "Differentiable programming" => "differentiable_programming.jl", From 86bf70b84470d3474b4177d645adf23e3b870b03 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 10 Jan 2024 11:42:50 +0100 Subject: [PATCH 35/69] typos --- docs/literate/src/files/p4est_from_gmsh.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index f7925e68026..dbd3b5952c8 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -34,7 +34,7 @@ end #hide #md # There are plenty of possibilities for the user, see the [documentation](https://gmsh.info/doc/texinfo/gmsh.html) and [tutorials](https://gitlab.onelab.info/gmsh/gmsh/tree/master/tutorials). # To begin, we provide a complete geometry file in this tutorial. After this we give a breakdown -# of the most important parts required for succesfull mesh generation that can later be used by the `p4est` library +# of the most important parts required for successful mesh generation that can later be used by the `p4est` library # and `Trixi.jl`. # The associated `NACA6412.geo` file is given below: @@ -311,7 +311,7 @@ end #hide #md # This is strictly required to be able to use the mesh later with `p4est` with supports only straight-sided quads, # i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. # See for more details the (short) ]documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. -# In principle, it should also be possible to use the recombine function of `gmsh` to convert the traingles to quads, +# In principle, it should also be possible to use the recombine function of `gmsh` to convert the triangles to quads, # but this is observed to be less robust than enforcing quads from the beginning. # # Then the airfoil is defined by a set of points: From c0e6be8e3ddffde9ceb9f1ceaa1f38942d35a63f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:24:16 +0100 Subject: [PATCH 36/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 2a48da2a5b9..42d189f1e96 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -243,7 +243,7 @@ In addition to the list of [nodes](#nodes) and [elements](#elements) given above 1, 4, 52, 53, 54, 55, 56, 57, 58, ``` required which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. -By looping over every element and its associated edges consisting of two nodes, we query the read in `NSET`s if the current node pair is present. +By looping over every element and its associated edges, consisting of two nodes, we query the read in `NSET`s if the current node pair is present. To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user has to supply the `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: From 5c65d704605377cbebbf730a14c1510818b8f40f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:24:29 +0100 Subject: [PATCH 37/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index dbd3b5952c8..dd756aea087 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -311,7 +311,7 @@ end #hide #md # This is strictly required to be able to use the mesh later with `p4est` with supports only straight-sided quads, # i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. # See for more details the (short) ]documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. -# In principle, it should also be possible to use the recombine function of `gmsh` to convert the triangles to quads, +# In principle, it should also be possible to use the `recombine` function of `gmsh` to convert the triangles to quads, # but this is observed to be less robust than enforcing quads from the beginning. # # Then the airfoil is defined by a set of points: From c9ada6fd958799836d443ef8c51892ddecfeeabe Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:24:37 +0100 Subject: [PATCH 38/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index dd756aea087..aba5890f64d 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -310,7 +310,7 @@ end #hide #md # which forces `gmsh` to generate quadrilateral elements instead of the default triangles. # This is strictly required to be able to use the mesh later with `p4est` with supports only straight-sided quads, # i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. -# See for more details the (short) ]documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. +# See for more details the (short) [documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. # In principle, it should also be possible to use the `recombine` function of `gmsh` to convert the triangles to quads, # but this is observed to be less robust than enforcing quads from the beginning. # From d3e93876e28d7ec52a00669c8f394f7df00f412a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:24:52 +0100 Subject: [PATCH 39/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index aba5890f64d..05105d46699 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -308,7 +308,7 @@ end #hide #md # Mesh.RecombineAll = 1; # ``` # which forces `gmsh` to generate quadrilateral elements instead of the default triangles. -# This is strictly required to be able to use the mesh later with `p4est` with supports only straight-sided quads, +# This is strictly required to be able to use the mesh later with `p4est`, which supports only straight-sided quads, # i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. # See for more details the (short) [documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. # In principle, it should also be possible to use the `recombine` function of `gmsh` to convert the triangles to quads, From 7152cfa09055fa387aa9f14de1a20d58d0d45f07 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:25:06 +0100 Subject: [PATCH 40/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 42d189f1e96..48d9576152b 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -253,7 +253,7 @@ boundary_symbols = [:PhysicalLine1] mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) ``` By doing so, only nodesets with label in `boundary_symbols` are treated as boundaries and other nodesets which could be used for diagnostics are not treated as boundary. -Note that there is a leading double-colon `:` compared to the label in the `.inp` mesh file. +Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). A 2D example for this mesh generation is presented in From 38ca0199568a3d7441b351fa7bd059c536f62471 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:25:25 +0100 Subject: [PATCH 41/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 48d9576152b..a3380d39a2c 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -252,7 +252,8 @@ boundary_symbols = [:PhysicalLine1] mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) ``` -By doing so, only nodesets with label in `boundary_symbols` are treated as boundaries and other nodesets which could be used for diagnostics are not treated as boundary. +By doing so, only nodesets with a label present in `boundary_symbols` are treated as physical boundaries. +Other nodesets that could be used for diagnostics are not treated as external boundary. Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). From 6be55ee05733be51a6c220de9498c6ec43c1f667 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:25:41 +0100 Subject: [PATCH 42/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index a3380d39a2c..90cf6595409 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -245,7 +245,7 @@ In addition to the list of [nodes](#nodes) and [elements](#elements) given above required which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. By looping over every element and its associated edges, consisting of two nodes, we query the read in `NSET`s if the current node pair is present. -To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user has to supply the `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: +To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user must supply a `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: ```julia boundary_symbols = [:PhysicalLine1] From 22b7283673f452010a53bfb77fa48c79efe6192c Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:26:05 +0100 Subject: [PATCH 43/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 90cf6595409..a71b407cd98 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -257,7 +257,7 @@ Other nodesets that could be used for diagnostics are not treated as external bo Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). -A 2D example for this mesh generation is presented in +A 2D example for this mesh read-in for an unstructured mesh file created with `gmsh` is presented in `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`. ### Mesh in three spatial dimensions From 1d7f5a59bf51e5967f6180c6dd3b48f16802294d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:26:20 +0100 Subject: [PATCH 44/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index a71b407cd98..6cfbcdb90a0 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -387,7 +387,7 @@ The most notable difference is that boundaries are formed in 3D by faces defined A simple mesh file used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl` is given below: ``` *Heading - Something DIFFERENT FROM ** ***** HOHQMesh boundary information ***** ** + *NODE 1, -2, 0, 0 2, -1, 0, 0 From a3ca4e45808f7acb4a9b9ec11022c3b387fbe42f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:26:33 +0100 Subject: [PATCH 45/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Andrew Winters --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 05105d46699..25c0b76262c 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -20,7 +20,7 @@ redirect_stdio(stdout=devnull, stderr=devnull) do # code that prints annoying st trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_euler_NACA6412airfoil_mach2.jl"), tspan=(0.0, 0.5)) end #hide #md -# Conveniently, we can use the Plots package to have a first look at the results: +# Conveniently, we use the Plots package to have a first look at the results: # ```julia # using Plots # pd = PlotData2D(sol) From 0b99b49b1120b97993acd92a9d5a93296d541cfc Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:26:52 +0100 Subject: [PATCH 46/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Andrew Winters --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 25c0b76262c..735cef6588b 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -33,7 +33,7 @@ end #hide #md # The creation of an unstructured quadrilateral mesh using `gmsh` is driven by a **geometry file**. # There are plenty of possibilities for the user, see the [documentation](https://gmsh.info/doc/texinfo/gmsh.html) and [tutorials](https://gitlab.onelab.info/gmsh/gmsh/tree/master/tutorials). -# To begin, we provide a complete geometry file in this tutorial. After this we give a breakdown +# To begin, we provide a complete geometry file for the NACA6412 airfoil bounded by a rectangular box in this tutorial. After this we give a breakdown # of the most important parts required for successful mesh generation that can later be used by the `p4est` library # and `Trixi.jl`. From 013ff843452d95e9e3a88af8d3a22d518ec5d35e Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:30:13 +0100 Subject: [PATCH 47/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Andrew Winters --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 6cfbcdb90a0..e79ba8c58fc 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -242,7 +242,7 @@ In addition to the list of [nodes](#nodes) and [elements](#elements) given above *NSET,NSET=PhysicalLine1 1, 4, 52, 53, 54, 55, 56, 57, 58, ``` -required which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. +which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. By looping over every element and its associated edges, consisting of two nodes, we query the read in `NSET`s if the current node pair is present. To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user must supply a `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: From 0622b7a8365496b71323a95c021250235a226530 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:31:51 +0100 Subject: [PATCH 48/69] Update docs/src/meshes/p4est_mesh.md --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index e79ba8c58fc..53f365ca396 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -242,7 +242,7 @@ In addition to the list of [nodes](#nodes) and [elements](#elements) given above *NSET,NSET=PhysicalLine1 1, 4, 52, 53, 54, 55, 56, 57, 58, ``` -which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. +present which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. By looping over every element and its associated edges, consisting of two nodes, we query the read in `NSET`s if the current node pair is present. To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user must supply a `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: From dc2f232be46998fae76a74514722a3fd41262149 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:34:24 +0100 Subject: [PATCH 49/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 735cef6588b..898c6bad48c 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -340,10 +340,10 @@ end #hide #md # At the end of the file the physical groups are defined: # ```c++ # // labeling of the boundary parts -# Physical Line(1) = {4}; // inflow -# Physical Line(2) = {2}; // outflow -# Physical Line(3) = {1, 3}; // airfoil -# Physical Line(4) = {5, 6}; // upper/lower wall +# Physical Line(1) = {4}; // Inflow. Label in Abaqus .inp file: PhysicalLine1 +# Physical Line(2) = {2}; // Outflow. Label in Abaqus .inp file: PhysicalLine2 +# Physical Line(3) = {1, 3}; // Airfoil. Label in Abaqus .inp file: PhysicalLine3 +# Physical Line(4) = {5, 6}; //Upper and lower wall/farfield/... Label in Abaqus .inp file: PhysicalLine4 # ``` # which are crucial for the correct assignment of boundary conditions in `Trixi.jl`. # From 77f22a6706309acf1c3347cafd7dd24348e2acdd Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:37:19 +0100 Subject: [PATCH 50/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 898c6bad48c..0de705335e5 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -346,6 +346,7 @@ end #hide #md # Physical Line(4) = {5, 6}; //Upper and lower wall/farfield/... Label in Abaqus .inp file: PhysicalLine4 # ``` # which are crucial for the correct assignment of boundary conditions in `Trixi.jl`. +In particular, it is the users responsibility to keep track on the physical boundary names between the mesh generation and assignment of boundary condition functions in an elixir. # # After opening this file in `gmsh`, meshing the geometry and exporting to Abaqus `.inp` format, # we can have a look at the input file: From f9af6350d425e13e1a43e20e4355d410998984c2 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:40:56 +0100 Subject: [PATCH 51/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 0de705335e5..66201a4e5bd 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -65,6 +65,7 @@ end #hide #md # Mesh.ColorCarousel = 0; # # // points of the airfoil contour +# // Format: {x, y, z, DesiredCellSize}. See the [documentation](https://gmsh.info/doc/texinfo/gmsh.html#Points). # Point(5) = {-0.4900332889206208, 0.09933466539753061, 0, 0.125}; # Point(6) = {-0.4900274857651495, 0.1021542752054094, 0, 0.125}; # Point(7) = {-0.4894921489729144, 0.1049830248247787, 0, 0.125}; From 6405b0aa3b3758c1791e75936eb35abff152a13a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 17:45:39 +0100 Subject: [PATCH 52/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 66201a4e5bd..8e1a9ba8f62 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -65,7 +65,7 @@ end #hide #md # Mesh.ColorCarousel = 0; # # // points of the airfoil contour -# // Format: {x, y, z, DesiredCellSize}. See the [documentation](https://gmsh.info/doc/texinfo/gmsh.html#Points). +# // Format: {x, y, z, DesiredCellSize}. See the documentation: https://gmsh.info/doc/texinfo/gmsh.html#Points # Point(5) = {-0.4900332889206208, 0.09933466539753061, 0, 0.125}; # Point(6) = {-0.4900274857651495, 0.1021542752054094, 0, 0.125}; # Point(7) = {-0.4894921489729144, 0.1049830248247787, 0, 0.125}; From 5d43f8bea191a98fa722510747032a9b7ea5a299 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 13 Jan 2024 18:53:21 +0100 Subject: [PATCH 53/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 8e1a9ba8f62..33e4fb7cf58 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -347,7 +347,7 @@ end #hide #md # Physical Line(4) = {5, 6}; //Upper and lower wall/farfield/... Label in Abaqus .inp file: PhysicalLine4 # ``` # which are crucial for the correct assignment of boundary conditions in `Trixi.jl`. -In particular, it is the users responsibility to keep track on the physical boundary names between the mesh generation and assignment of boundary condition functions in an elixir. +# In particular, it is the users responsibility to keep track on the physical boundary names between the mesh generation and assignment of boundary condition functions in an elixir. # # After opening this file in `gmsh`, meshing the geometry and exporting to Abaqus `.inp` format, # we can have a look at the input file: From 8ff47e73c745d4260ae7a15d05a5756401320677 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sun, 14 Jan 2024 11:31:13 +0100 Subject: [PATCH 54/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 33e4fb7cf58..68f9cd2ab16 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -39,7 +39,8 @@ end #hide #md # The associated `NACA6412.geo` file is given below: # ```c++ -# // GMSH geometry script for a NACA 6412 airfoil in a box +# // GMSH geometry script for a NACA 6412 airfoil with 11 degree angle of attack +# // in a box (near-field mesh) # // see https://github.com/cfsengineering/GMSH-Airfoil-2D # # // outer bounding box From 246878cf79ae6d24e4448d7fda1d6df66156560e Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 16 Jan 2024 11:34:16 +0100 Subject: [PATCH 55/69] Apply suggestions from code review Polish --- docs/literate/src/files/p4est_from_gmsh.jl | 29 ++++++++++++---------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 68f9cd2ab16..52ffe1e90b4 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -11,7 +11,7 @@ # ## Running the simulation of a near-field flow around an airfoil -# Trixi.jl supports solving hyperbolic problems on several mesh types. +# Trixi.jl supports solving hyperbolic-parabolic problems on several mesh types. # A somewhat complex example that employs the `P4estMesh` is the near-field simulation of a # Mach 2 flow around the NACA6412 airfoil. @@ -33,15 +33,17 @@ end #hide #md # The creation of an unstructured quadrilateral mesh using `gmsh` is driven by a **geometry file**. # There are plenty of possibilities for the user, see the [documentation](https://gmsh.info/doc/texinfo/gmsh.html) and [tutorials](https://gitlab.onelab.info/gmsh/gmsh/tree/master/tutorials). -# To begin, we provide a complete geometry file for the NACA6412 airfoil bounded by a rectangular box in this tutorial. After this we give a breakdown +# To begin, we provide a complete geometry file for the NACA6412 airfoil bounded by a rectangular box. After this we give a breakdown # of the most important parts required for successful mesh generation that can later be used by the `p4est` library # and `Trixi.jl`. +# We emphasize that this near-field mesh should only be used for instructive purposes and not for actual production runs. # The associated `NACA6412.geo` file is given below: # ```c++ # // GMSH geometry script for a NACA 6412 airfoil with 11 degree angle of attack -# // in a box (near-field mesh) +# // in a box (near-field mesh). # // see https://github.com/cfsengineering/GMSH-Airfoil-2D +# // for software to generate gmsh `.geo` geometry files for NACA airfoils. # # // outer bounding box # Point(1) = {-1.25, -0.5, 0, 1.0}; @@ -365,20 +367,20 @@ end #hide #md # *ELEMENT, type=T3D2, ELSET=Line1 # 1, 1, 7 # ... -# *ELEMENT, type=CPS4, ELSET=Surface10 +# *ELEMENT, type=CPS4, ELSET=Surface1 # 191, 272, 46, 263, 807 # ... -# *NSET,NSET=PhysicalLine10 +# *NSET,NSET=PhysicalLine1 # 1, 4, 52, 53, 54, 55, 56, 57, 58, -# *NSET,NSET=PhysicalLine20 +# *NSET,NSET=PhysicalLine2 # 2, 3, 26, 27, 28, 29, 30, 31, 32, -# *NSET,NSET=PhysicalLine30 +# *NSET,NSET=PhysicalLine3 # 1, 2, 3, 4, 7, 8, 9, 10, 11, 12, # 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, # 23, 24, 25, 33, 34, 35, 36, 37, 38, 39, # 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, # 50, 51, -# *NSET,NSET=PhysicalLine40 +# *NSET,NSET=PhysicalLine4 # 5, 6, 59, 60, 61, 62, 63, 64, 65, 66, # 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, # 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, @@ -403,7 +405,7 @@ end #hide #md # the supplied `boundary_symbols` that have to be supplied to the `P4estMesh` constructor: # ```julia # # boundary symbols -# boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] +# boundary_symbols = [:PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4] # mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) # ``` # The same boundary symbols have then also be supplied to the semidiscretization alongside the @@ -432,16 +434,17 @@ end #hide #md # equations::CompressibleEulerEquations2D) # flux = Trixi.flux(u_inner, normal_direction, equations) # -# boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, # Left boundary -# :PhysicalLine20 => boundary_condition_supersonic_outflow, # Right boundary -# :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil -# :PhysicalLine40 => boundary_condition_supersonic_outflow) # Top and bottom boundary +# boundary_conditions = Dict(:PhysicalLine1 => boundary_condition_supersonic_inflow, # Left boundary +# :PhysicalLine2 => boundary_condition_supersonic_outflow, # Right boundary +# :PhysicalLine3 => boundary_condition_slip_wall, # Airfoil +# :PhysicalLine4 => boundary_condition_supersonic_outflow) # Top and bottom boundary # # semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, # boundary_conditions = boundary_conditions) # ``` # Note that you **have to** supply the `boundary_symbols` keyword to the `P4estMesh` constructor # to select the boundaries from the available nodesets in the `.inp` file. +# If the `boundary_symbols` keyword is not supplied, all boundaries will be assigned to the default set `:all`. # ## Package versions From 9de58c14cd4f648422695f333e353ddedfdf4da4 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 16 Jan 2024 22:35:20 +0100 Subject: [PATCH 56/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Andrew Winters --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 52ffe1e90b4..e64248cee8a 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -350,7 +350,7 @@ end #hide #md # Physical Line(4) = {5, 6}; //Upper and lower wall/farfield/... Label in Abaqus .inp file: PhysicalLine4 # ``` # which are crucial for the correct assignment of boundary conditions in `Trixi.jl`. -# In particular, it is the users responsibility to keep track on the physical boundary names between the mesh generation and assignment of boundary condition functions in an elixir. +# In particular, it is the responsibility of a user to keep track on the physical boundary names between the mesh generation and assignment of boundary condition functions in an elixir. # # After opening this file in `gmsh`, meshing the geometry and exporting to Abaqus `.inp` format, # we can have a look at the input file: From 30b2946fc7fb0a924a0a8cab011c3c6e54f169e2 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 17 Jan 2024 23:19:59 +0100 Subject: [PATCH 57/69] Comments --- src/meshes/p4est_mesh.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 25e4dbef566..5f6609cae6b 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -519,10 +519,10 @@ function parse_elements(meshfile, n_trees, n_dims) if length(content) == expected_content_length # Check that we still read in connectivity data content_int = parse.(Int64, content) # Add constituent nodes to the element_node_matrix. - # Important: Do not use idex form the Abaqus file, but the one from p4est + # Important: Do not use index form the Abaqus file, but the one from p4est element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id tree_id += 1 - else # Read all elements for this ELSET + else # Processed all elements for this ELSET el_list_follows = false end end @@ -568,7 +568,7 @@ function parse_node_sets(meshfile, boundary_symbols) end # This function assigns the edges of elements to boundaries by -# checking if the nodes of the edges are part of nodesets which correspond to boundaries. +# checking if the nodes that define the edges are part of nodesets which correspond to boundaries. function assign_boundaries_standard_abaqus!(boundary_names, n_trees, element_node_matrix, node_set_dict, ::Val{2}) # 2D version @@ -609,7 +609,7 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, end # This function assigns the edges of elements to boundaries by -# checking if the nodes of the edges are part of nodesets which correspond to boundaries. +# checking if the nodes that define the faces are part of nodesets which correspond to boundaries. function assign_boundaries_standard_abaqus!(boundary_names, n_trees, element_node_matrix, node_set_dict, ::Val{3}) # 3D version From 253c5025547c91e1ca3aa299fc0456d4472ea552 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 18 Jan 2024 13:26:20 +0100 Subject: [PATCH 58/69] Update p4est_mesh.jl --- src/meshes/p4est_mesh.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 5f6609cae6b..08ea0f6dd09 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -311,8 +311,9 @@ To create a curved unstructured mesh `P4estMesh` two strategies are available: straight-sided 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`, as distinguishing different physical boundaries is - non-trivial. + 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 `p4est_mesh_from_standard_abaqus` function. The `p4est_mesh_from_hohqmesh_abaqus` function obtains the mesh `polydeg` directly from the `meshfile` From d7292356f265216229512a2562585425deb0642a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 12:16:00 +0100 Subject: [PATCH 59/69] Update src/meshes/p4est_mesh.jl Co-authored-by: Johannes Markert <10619309+jmark@users.noreply.github.com> --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 08ea0f6dd09..25fd10da854 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -520,7 +520,7 @@ function parse_elements(meshfile, n_trees, n_dims) if length(content) == expected_content_length # Check that we still read in connectivity data content_int = parse.(Int64, content) # Add constituent nodes to the element_node_matrix. - # Important: Do not use index form the Abaqus file, but the one from p4est + # Important: Do not use index from the Abaqus file, but the one from p4est. element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id tree_id += 1 else # Processed all elements for this ELSET From b98172728504969c389ff4ce0f2c52765dfe3dff Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 12:16:33 +0100 Subject: [PATCH 60/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Johannes Markert <10619309+jmark@users.noreply.github.com> --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 53f365ca396..0c215c1befb 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -242,7 +242,7 @@ In addition to the list of [nodes](#nodes) and [elements](#elements) given above *NSET,NSET=PhysicalLine1 1, 4, 52, 53, 54, 55, 56, 57, 58, ``` -present which are used to associate the edges defined through their corner nodes with a label, in this case `PhysicalLine1`. +present which are used to associate the edges defined through their corner nodes with a label. In this case it is called `PhysicalLine1`. By looping over every element and its associated edges, consisting of two nodes, we query the read in `NSET`s if the current node pair is present. To prevent that every nodeset following `*NSET,NSET=` is treated as a boundary, the user must supply a `boundary_symbols` keyword to the [`P4estMesh`](@ref) constructor: From 9a116e57d4b75679e5c0f2608742c985410bc3a5 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 12:16:44 +0100 Subject: [PATCH 61/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Johannes Markert <10619309+jmark@users.noreply.github.com> --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 0c215c1befb..2aae1ae6651 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -253,7 +253,7 @@ boundary_symbols = [:PhysicalLine1] mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols) ``` By doing so, only nodesets with a label present in `boundary_symbols` are treated as physical boundaries. -Other nodesets that could be used for diagnostics are not treated as external boundary. +Other nodesets that could be used for diagnostics are not treated as external boundaries. Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). From 43d0ac917b7b28909b25b6cc478548b7090f7104 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 12:17:03 +0100 Subject: [PATCH 62/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Johannes Markert <10619309+jmark@users.noreply.github.com> --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 2aae1ae6651..b9bb3659bb4 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -257,7 +257,7 @@ Other nodesets that could be used for diagnostics are not treated as external bo Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). -A 2D example for this mesh read-in for an unstructured mesh file created with `gmsh` is presented in +A 2D example for this mesh, which is read-in for an unstructured mesh file created with `gmsh`, is presented in `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`. ### Mesh in three spatial dimensions From 282ba8241e81aea50d0a9bfc2e70c7b0b3bcc52b Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 12:17:40 +0100 Subject: [PATCH 63/69] Update docs/src/meshes/p4est_mesh.md Co-authored-by: Johannes Markert <10619309+jmark@users.noreply.github.com> --- docs/src/meshes/p4est_mesh.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index b9bb3659bb4..db75346cab3 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -384,7 +384,7 @@ transfinite map of the straight sided hexahedral element to find Also for a mesh in standard Abaqus format there are no qualitative changes when going from 2D to 3D. The most notable difference is that boundaries are formed in 3D by faces defined by four nodes while in 2D boundaries are edges consisting of two elements. -A simple mesh file used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl` is given below: +A simple mesh file, which is used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl`, is given below: ``` *Heading From 6ad175e9a51aa120371a5060ca7d71d96504d25f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 23 Jan 2024 12:48:21 +0100 Subject: [PATCH 64/69] more robust parsing --- src/meshes/p4est_mesh.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 25fd10da854..71bdcd91467 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -555,8 +555,14 @@ function parse_node_sets(meshfile, boundary_symbols) current_symbol = nothing end elseif current_symbol !== nothing # Read only if there was already a nodeset specified - # There is always a trailing comma, remove the corresponding empty string - append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)])) + try # Check if line contains nodes + # There is always a trailing comma, remove the corresponding empty string + append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)])) + catch # Something different, stop reading in nodes + # If parsing fails, set current_symbol to nothing + nodes_dict[current_symbol] = current_nodes + current_symbol = nothing + end end end # Safe the previous nodeset From 9676b90e55f1906bfe881b6dfb027114c2fd2b68 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 23 Jan 2024 13:02:57 +0100 Subject: [PATCH 65/69] print warning if there has been a bc name supplied for which no nodes have been found --- src/meshes/p4est_mesh.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 71bdcd91467..d8376eac29e 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -571,6 +571,12 @@ function parse_node_sets(meshfile, boundary_symbols) end end + for symbol in boundary_symbols + if !haskey(nodes_dict, symbol) + @warn "No nodes found for nodeset :" * "$symbol" * " !" + end + end + return nodes_dict end From 964c2f825e07a2b6fea07ca45e54d9d0fd445692 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 16:01:57 +0100 Subject: [PATCH 66/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Hendrik Ranocha --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index e64248cee8a..be3ceac2c63 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -35,7 +35,7 @@ end #hide #md # To begin, we provide a complete geometry file for the NACA6412 airfoil bounded by a rectangular box. After this we give a breakdown # of the most important parts required for successful mesh generation that can later be used by the `p4est` library -# and `Trixi.jl`. +# and Trixi.jl. # We emphasize that this near-field mesh should only be used for instructive purposes and not for actual production runs. # The associated `NACA6412.geo` file is given below: From 8ae228d4dc9472511ff2059c7f3788ab3bf0dcf0 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 16:02:08 +0100 Subject: [PATCH 67/69] Update docs/literate/src/files/p4est_from_gmsh.jl Co-authored-by: Hendrik Ranocha --- docs/literate/src/files/p4est_from_gmsh.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index be3ceac2c63..3079e401f33 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -400,8 +400,8 @@ end #hide #md # First, the coordinates of the nodes are listed, followed by the elements. # Note that `gmsh` exports also line elements of type `T3D2` which are ignored by `p4est`. # The relevant elements in 2D which form the gridcells are of type `CPS4` which are defined by their four corner nodes. -# This is followed by the nodesets encoded via `*NSET` which are used to assign boundary conditions in `Trixi.jl`. -# `Trixi.jl` parses the `.inp` file and assigns the edges (in 2D, surfaces in 3D) of elements to the corresponding boundary condition based on +# This is followed by the nodesets encoded via `*NSET` which are used to assign boundary conditions in Trixi.jl. +# Trixi.jl parses the `.inp` file and assigns the edges (in 2D, surfaces in 3D) of elements to the corresponding boundary condition based on # the supplied `boundary_symbols` that have to be supplied to the `P4estMesh` constructor: # ```julia # # boundary symbols From 16a7b25e08b28bca2371f279d431f450678bf941 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 16:02:26 +0100 Subject: [PATCH 68/69] Update src/meshes/p4est_mesh.jl Co-authored-by: Hendrik Ranocha --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index d8376eac29e..c5d39ef00c0 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -499,7 +499,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, end function parse_elements(meshfile, n_trees, n_dims) - @assert n_dims in [2, 3] "Only 2D and 3D meshes are supported" + @assert n_dims in (2, 3) "Only 2D and 3D meshes are supported" # Valid element types (that can be processed by p4est) based on dimension element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", From 85817af5a80640ee42b87480697032897e04cfa4 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 23 Jan 2024 16:05:39 +0100 Subject: [PATCH 69/69] Update docs/literate/src/files/p4est_from_gmsh.jl --- docs/literate/src/files/p4est_from_gmsh.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index 3079e401f33..356339cdd47 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -69,6 +69,7 @@ end #hide #md # # // points of the airfoil contour # // Format: {x, y, z, DesiredCellSize}. See the documentation: https://gmsh.info/doc/texinfo/gmsh.html#Points +# // These concrete points are generated using the tool from https://github.com/cfsengineering/GMSH-Airfoil-2D # Point(5) = {-0.4900332889206208, 0.09933466539753061, 0, 0.125}; # Point(6) = {-0.4900274857651495, 0.1021542752054094, 0, 0.125}; # Point(7) = {-0.4894921489729144, 0.1049830248247787, 0, 0.125};