P4estMesh constructor for HOHQMesh Abaqus-type file

examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl

@@ -0,0 +1,97 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+equations = CompressibleEulerEquations2D(1.4)
+
+@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D)
+
+  # set the freestream flow parameters
+  rho_freestream = 1.0
+  u_freestream = 0.3
+  p_freestream = inv(equations.gamma)
+
+  theta = pi / 90.0 # analogous with a two degree angle of attack
+  si, co = sincos(theta)
+  v1 = u_freestream * co
+  v2 = u_freestream * si
+
+  prim = SVector(rho_freestream, v1, v2, p_freestream)
+  return prim2cons(prim, equations)
+end
+
+initial_condition = uniform_flow_state
+
+boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
+boundary_conditions = Dict( :Body    => boundary_condition_uniform_flow,
+                            :Button1 => boundary_condition_slip_wall,
+                            :Button2 => boundary_condition_slip_wall,
+                            :Eye1    => boundary_condition_slip_wall,
+                            :Eye2    => boundary_condition_slip_wall,
+                            :Smile   => boundary_condition_slip_wall,
+                            :Bowtie  => boundary_condition_slip_wall )
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=5, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "abaqus_gingerbread_man.inp")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/0e9e990a04b5105d1d2e3096a6e41272/raw/67c2f7f889e637292f551284777c368590d2f12b/abaqus_gingerbread_man.inp",
+                                      default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = P4estMesh{2}(mesh_file, true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+amr_indicator = IndicatorLöhner(semi, variable=density)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level=0,
+                                      med_level=1, med_threshold=0.05,
+                                      max_level=3, max_threshold=0.1)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        amr_callback)
+
+
+###############################################################################
+# run the simulation
+sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7,
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/meshes/p4est_mesh.jl

@@ -300,6 +300,85 @@ function P4estMesh{NDIMS}(meshfile::String;
 end
 
 
+# TODO: FIX ME, there is probably a better way to dispatch correctly on this mesh constructor routine
+"""
+    P4estMesh{NDIMS}(meshfile::String, fromHOHQMesh::Bool; RealT=Float64,
+                     initial_refinement_level=0, unsaved_changes=true)
+
+Import an unstructured conforming mesh from an Abaqus-style mesh file (`.inp`). High-order, curved
+boundary information as well as boundary names are provided by [`HOHQMesh`](@ref). This creates
+an unstructured, curved `P4estMesh` from the curved mesh.
+
+The polynomial degree of the mesh boundaries is provided from the `meshfile`. The computation of
+the mapped element coordinates is done with transfinite interpolation with linear blending similar
+to [`UnstructuredMesh2D`](@ref).
+
+# Arguments
+- `meshfile::String`: an uncurved Abaqus mesh file that was generated using [`HOHQMesh`](@ref)
+                      to be imported by p4est. Additional curved boundary information and boundary
+                      names are provided in the second part of the meshfile.
+- `RealT::Type`: the type that should be used for coordinates.
+- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts.
+- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
+"""
+function P4estMesh{NDIMS}(meshfile::String, fromHOHQMesh::Bool;
+                          RealT=Float64,
+                          initial_refinement_level=0,
+                          unsaved_changes=true) where NDIMS
+  # Prevent p4est from crashing Julia if the file doesn't exist
+  @assert isfile(meshfile)
+
+  # Have p4est create the mesh connectivity from the Abaqus portion of the meshfile
+  conn = read_inp_p4est(meshfile, Val(NDIMS))
+
+  # These need to be of the type Int for unsafe_wrap below to work
+  n_trees::Int = conn.num_trees
+  n_vertices::Int = conn.num_vertices
+
+  # Extract a copy of the element corners to compute the tree node coordinates
+  vertices = unsafe_wrap(Array, conn.vertices, (3, n_vertices))
+
+  # readin all the information from the mesh file into a string array
+  file_lines = readlines(open(meshfile))
+
+  # Set the starting file index to read in the additional information provided by HOHQMesh.
+  # The six is to account for the other headers in the meshfile.
+  file_idx = 6 + n_trees + n_vertices
+
+  # Get the polynomial order of the mesh boundary information
+  current_line = split(file_lines[file_idx])
+  mesh_polydeg = parse(Int, current_line[5])
+  mesh_nnodes = mesh_polydeg + 1
+
+  # Create the Chebyshev-Gauss-Lobatto nodes used by HOHQMesh to represent the boundaries
+  cheby_nodes_, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes)
+  nodes = SVector{mesh_nnodes}(cheby_nodes_)
+
+  # Allocate the memory for the tree node coordinates
+  tree_node_coordinates = Array{RealT, NDIMS+2}(undef, NDIMS,
+                                                ntuple(_ -> length(nodes), NDIMS)...,
+                                                n_trees)
+
+  # Compute the tree node coordinates and return the updated file index
+  file_idx = calc_tree_node_coordinates!(tree_node_coordinates, file_lines, nodes, vertices, RealT)
+
+  # Allocate the memory for the boundary labels
+  boundary_names = Array{Symbol}(undef, (4, n_trees))
+
+  # Read in the boundary names from the last portion of the meshfile
+  # Note here the boundary names where "---" means an internal connection
+  for tree in 1:n_trees
+    boundary_names[:, tree] = map(Symbol, split(file_lines[file_idx]))
+    file_idx  += 1
+  end
+
+  p4est = new_p4est(conn, initial_refinement_level)
+
+  return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes,
+                          boundary_names, "", unsaved_changes)
+end
+
+
 """
     P4estMesh(trees_per_face_dimension, layers, inner_radius, thickness;
               polydeg, RealT=Float64,
@@ -1021,6 +1100,110 @@ function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction)
 end
 
 
+# Calculate physical coordinates of each element of an unstructured mesh read
+# in from a HOHQMesh file. This calculation is done with the transfinite interpolation
+# routines found in `mappings_geometry_curved_2d.jl` or `mappings_geometry_straight_2d.jl`
+function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4},
+                                     file_lines::Vector{String}, nodes, vertices, RealT)
+  # Get the number of trees and the number of interpolation nodes
+  n_trees = last(size(node_coordinates))
+  nnodes = length(nodes)
+
+  # Move the file index forward to start reading in the element information
+  # Set the starting file index to start reading in the element information
+  # The seven account for the other headers in the mesh file.
+  file_idx = 7 + n_trees + last(size(vertices))
+
+  # Create a work set of Gamma curves to create the node coordinates
+  CurvedSurfaceT    = CurvedSurface{RealT}
+  surface_curves    = Array{CurvedSurfaceT}(undef, 4)
+
+  # Create other work arrays to perform the mesh construction
+  element_node_ids  = Array{Int}(undef, 4)
+  curved_check      = Vector{Int}(undef, 4)
+  cornerNodeVals    = Array{RealT}(undef, (4, 2))
+  tempNodes         = Array{RealT}(undef, (4, 2))
+  curve_vals        = Array{RealT}(undef, (nnodes, 2))
+
+  # Create the barycentric weights used for the surface interpolations
+  bary_weights_ = barycentric_weights(nodes)
+  bary_weights = SVector{nnodes}(bary_weights_)
+
+  # Loop through all the trees, i.e., the elements generated by HOHQMesh and create the node coordinates
+  for tree in 1:n_trees
+    # pull the corner node IDs
+    current_line        = split(file_lines[file_idx])
+    element_node_ids[1] = parse(Int, current_line[1])
+    element_node_ids[2] = parse(Int, current_line[2])
+    element_node_ids[3] = parse(Int, current_line[3])
+    element_node_ids[4] = parse(Int, current_line[4])
+
+    # pull the (x,y) values of the four corners of the current tree out of the vertices array
+    for i in 1:4
+      cornerNodeVals[i, :] .= vertices[1:2, element_node_ids[i]]
+    end
+    # pull the information to check if boundary is curved in order to read in additional data
+    file_idx += 1
+    current_line    = split(file_lines[file_idx])
+    curved_check[1] = parse(Int, current_line[1])
+    curved_check[2] = parse(Int, current_line[2])
+    curved_check[3] = parse(Int, current_line[3])
+    curved_check[4] = parse(Int, current_line[4])
+    if sum(curved_check) == 0
+      # Create the node coordinates on this particular element
+      calc_node_coordinates!(node_coordinates, tree, nodes, cornerNodeVals)
+      # quadrilateral element is straight sided so we are ready to increment the file index
+      file_idx  += 1
+    else
+      # quadrilateral element has at least one curved side
+      # flip node ordering to make sure the element is right-handed for the interpolations
+      m1 = 1
+      m2 = 2
+      @views tempNodes[1, :] .= cornerNodeVals[4, :]
+      @views tempNodes[2, :] .= cornerNodeVals[2, :]
+      @views tempNodes[3, :] .= cornerNodeVals[3, :]
+      @views tempNodes[4, :] .= cornerNodeVals[1, :]
+      for i in 1:4
+        if curved_check[i] == 0
+          # when curved_check[i] is 0 then the "curve" from cornerNode(i) to cornerNode(i+1) is a
+          # straight line. So we must construct the interpolant for this line
+          for k in 1:nnodes
+            curve_vals[k, 1] = tempNodes[m1, 1] + 0.5 * (nodes[k] + 1.0) * ( tempNodes[m2, 1]
+                                                                            - tempNodes[m1, 1] )
+            curve_vals[k, 2] = tempNodes[m1, 2] + 0.5 * (nodes[k] + 1.0) * ( tempNodes[m2, 2]
+                                                                            - tempNodes[m1, 2] )
+          end
+        else
+          # when curved_check[i] is 1 this curved boundary information is supplied by the mesh
+          # generator. So we just read it into a work array
+          for k in 1:nnodes
+            file_idx += 1
+            current_line = split(file_lines[file_idx])
+            curve_vals[k, 1] = parse(RealT,current_line[1])
+            curve_vals[k, 2] = parse(RealT,current_line[2])
+          end
+        end
+        # construct the curve interpolant for the current side
+        surface_curves[i] = CurvedSurfaceT(nodes, bary_weights, copy(curve_vals))
+        # indexing update that contains a "flip" to ensure correct element orientation
+        # if we need to construct the straight line "curves" when curved_check[i] == 0
+        m1 += 1
+        if i == 3
+          m2 = 1
+        else
+          m2 += 1
+        end
+      end
+      # Create the node coordinates on this particular element
+      calc_node_coordinates!(node_coordinates, tree, nodes, surface_curves)
+      file_idx += 1
+    end
+  end
+
+  return file_idx
+end
+
+
 function balance!(mesh::P4estMesh{2}, init_fn=C_NULL)
   p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn)
   # Due to a bug in p4est, the forest needs to be rebalanced twice sometimes

test/test_p4est_2d.jl

@@ -85,6 +85,13 @@ isdir(outdir) && rm(outdir, recursive=true)
       tspan = (0.0, 0.3))
   end
 
+  @trixi_testset "elixir_euler_wall_bc_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc_amr.jl"),
+      l2   = [0.020291447969983396, 0.017479614254319948, 0.011387644425613437, 0.0514420126021293],
+      linf = [0.3582779022370579, 0.32073537890751663, 0.221818049107692, 0.9209559420400415],
+      tspan = (0.0, 0.15))
+  end
+
   @trixi_testset "elixir_eulergravity_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
       l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],