Skip to content

Commit

Permalink
3D
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed Jan 5, 2024
1 parent 8f0e184 commit c783247
Show file tree
Hide file tree
Showing 2 changed files with 136 additions and 33 deletions.
60 changes: 60 additions & 0 deletions examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl
View file
Original file line number Diff line number Diff line change
@@ -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
109 changes: 76 additions & 33 deletions src/meshes/p4est_mesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit c783247

Please sign in to comment.