Different Boundary Conditions for Standard Abaqus

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
 

examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl

@@ -0,0 +1,104 @@
+
+using Trixi
+using OrdinaryDiffEq
+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
+volume_flux = flux_ranocha
+
+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, # Left boundary
+                           :PhysicalLine20 => boundary_condition_outflow, # Right boundary
+                           :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil
+                           :PhysicalLine40 => boundary_condition_outflow) # Top and bottom boundary
+
+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
+###############################################################################
+sol = solve(ode, SSPRK104(; thread = OrdinaryDiffEq.True());
+            dt = 42.0, # overwritten
+            callback = callbacks);
+
+summary_callback() # print the timer summary

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

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)
 
@@ -373,7 +376,8 @@ function P4estMesh{NDIMS}(meshfile::String;
                                                                                               polydeg,
                                                                                               initial_refinement_level,
                                                                                               NDIMS,
-                                                                                              RealT)
+                                                                                              RealT,
+                                                                                              boundary_symbols)
     end
 
     return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes,
@@ -444,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)
@@ -469,12 +474,186 @@ 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, n_dimensions)
+
+        # 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
+    end
 
     return p4est, tree_node_coordinates, nodes, boundary_names
 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"]
+    # 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)
+                el_list_follows = true
+            elseif el_list_follows
+                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
+                    element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id
+                    tree_id += 1
+                else # Read all elements for this ELSET
+                    el_list_follows = false
+                end
+            end
+        end
+    end
+
+    return element_node_matrix
+end
+
+function parse_node_sets(meshfile, boundary_symbols)
+    nodes_dict = Dict{Symbol, Vector{Int64}}()
+    current_symbol = nothing
+    current_nodes = Int64[]
+
+    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, "*NSET,NSET=")
+                # Safe the previous nodeset
+                if current_symbol !== nothing
+                    nodes_dict[current_symbol] = current_nodes
+                end
+
+                current_symbol = Symbol(split(line, "=")[2])
+                if current_symbol in boundary_symbols
+                    # New nodeset
+                    current_nodes = Int64[]
+                else # Read only boundary node sets
+                    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)]))
+            end
+        end
+        # Safe the previous nodeset
+        if current_symbol !== nothing
+            nodes_dict[current_symbol] = current_nodes
+        end
+    end
+
+    return nodes_dict
+end
+
 """
     P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness;
                          polydeg, RealT=Float64,