P4estMesh constructor for HOHQMesh Abaqus-type file

examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl

@@ -0,0 +1,93 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler
+
+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/0d924b1d7e7d3cc1070a6cc22fe1d501687aa6dd/abaqus_gingerbread_man.inp",
+                                      default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = P4estMesh{2}(mesh_file)
+
+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_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_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

examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl

@@ -0,0 +1,69 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict( :Bottom => boundary_condition,
+                            :Top    => boundary_condition,
+                            :Circle => boundary_condition,
+                            :Cut    => boundary_condition )
+
+solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs)
+
+# Unstructured 3D half circle mesh from HOHQMesh
+default_mesh_file = joinpath(@__DIR__, "abaqus_half_circle_3d.inp")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/11461efbfb02c42e06aca338b3d0b645/raw/81deeb1ebc4945952c30af5bb75fe222a18d975c/abaqus_half_circle_3d.inp",
+                                      default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = P4estMesh{3}(mesh_file, initial_refinement_level=0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test,
+                                    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,
+                                     extra_analysis_integrals=(entropy,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=50,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback);
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1, # 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/auxiliary/p4est.jl

@@ -8,16 +8,16 @@
 """
     init_p4est()
 
-Initialize p4est by calling `p4est_init` and setting the log level to `SC_LP_ERROR`.
-This function will check if p4est is already initialized
+Initialize `p4est` by calling `p4est_init` and setting the log level to `SC_LP_ERROR`.
+This function will check if `p4est` is already initialized
 and if yes, do nothing, thus it is safe to call it multiple times.
 """
 function init_p4est()
   if p4est_package_id()[] >= 0
     return nothing
   end
 
-  # Initialize p4est with log level ERROR to prevent a lot of output in AMR simulations
+  # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations
   p4est_init(C_NULL, SC_LP_ERROR)
 
   return nothing
@@ -42,7 +42,7 @@ function unsafe_load_sc(::Type{T}, sc_array, i=1) where T
 end
 
 
-# Create new p4est from a p4est_connectivity
+# Create new `p4est` from a p4est_connectivity
 # 2D
 function new_p4est(conn::Ptr{p4est_connectivity_t}, initial_refinement_level)
   p4est_new_ext(0, # No MPI communicator
@@ -61,7 +61,7 @@ function new_p4est(conn::Ptr{p8est_connectivity_t}, initial_refinement_level)
 end
 
 
-# Save p4est data to file
+# Save `p4est` data to file
 # 2D
 function save_p4est!(file, p4est::Ptr{p4est_t})
   # Don't save user data of the quads
@@ -75,7 +75,7 @@ function save_p4est!(file, p8est::Ptr{p8est_t})
 end
 
 
-# Load p4est from file
+# Load `p4est` from file
 # 2D
 function load_p4est(file, ::Val{2})
   conn_vec = Vector{Ptr{p4est_connectivity_t}}(undef, 1)
@@ -89,28 +89,28 @@ function load_p4est(file, ::Val{3})
 end
 
 
-# Read p4est connectivity from Abaqus mesh file (.inp)
+# Read `p4est` connectivity from Abaqus mesh file (.inp)
 # 2D
 read_inp_p4est(meshfile, ::Val{2}) = p4est_connectivity_read_inp(meshfile)
 # 3D
 read_inp_p4est(meshfile, ::Val{3}) = p8est_connectivity_read_inp(meshfile)
 
 
-# Refine p4est if refine_fn_c returns 1
+# Refine `p4est` if refine_fn_c returns 1
 # 2D
 refine_p4est!(p4est::Ptr{p4est_t}, recursive, refine_fn_c, init_fn_c) = p4est_refine(p4est, recursive, refine_fn_c, init_fn_c)
 # 3D
 refine_p4est!(p8est::Ptr{p8est_t}, recursive, refine_fn_c, init_fn_c) = p8est_refine(p8est, recursive, refine_fn_c, init_fn_c)
 
 
-# Refine p4est if coarsen_fn_c returns 1
+# Refine `p4est` if coarsen_fn_c returns 1
 # 2D
 coarsen_p4est!(p4est::Ptr{p4est_t}, recursive, coarsen_fn_c, init_fn_c) = p4est_coarsen(p4est, recursive, coarsen_fn_c, init_fn_c)
 # 3D
 coarsen_p4est!(p8est::Ptr{p8est_t}, recursive, coarsen_fn_c, init_fn_c) = p8est_coarsen(p8est, recursive, coarsen_fn_c, init_fn_c)
 
 
-# Let p4est iterate over each cell volume and cell face.
+# Let `p4est` iterate over each cell volume and cell face.
 # Call iter_volume_c for each cell and iter_face_c for each face.
 # 2D
 function iterate_p4est(p4est::Ptr{p4est_t}, user_data;

src/meshes/face_interpolant.jl

@@ -0,0 +1,52 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+
+
+#     CurvedFace{RealT<:Real}
+#
+# Contains the data needed to represent a curved face with data points (x,y,z) as a Lagrange polynomial
+# interpolant written in barycentric form at a given set of nodes.
+struct CurvedFace{RealT<:Real}
+  nodes               ::Vector{RealT}
+  barycentric_weights ::Vector{RealT}
+  coordinates         ::Array{RealT, 3} #[ndims, nnodes, nnodes]
+end
+
+
+# evalute the Gamma face interpolant at a particular point s = (s_1, s_2) and return the (x,y,z) coordinate
+function evaluate_at(s, boundary_face::CurvedFace)
+
+   @unpack nodes, barycentric_weights, coordinates = boundary_face
+
+   x_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 1, :, :),
+                                                                  barycentric_weights)
+   y_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 2, :, :),
+                                                                  barycentric_weights)
+   z_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 3, :, :),
+                                                                  barycentric_weights)
+
+   return x_coordinate_at_s_on_boundary_face,
+          y_coordinate_at_s_on_boundary_face,
+          z_coordinate_at_s_on_boundary_face
+end
+
+
+# Calculate a 2D Lagrange interpolating polynomial in barycentric 2 form
+# of a function f(x,y) at a given coordinate (x,y) for a given node distribution.
+function lagrange_interpolation_2d(x, nodes, function_values, barycentric_weights)
+
+  f_intermediate = zeros(eltype(function_values), length(nodes))
+  for j in eachindex(nodes)
+    f_intermediate[j] = lagrange_interpolation(x[2], nodes, view(function_values, j, :),
+                                               barycentric_weights)
+  end
+  point_value = lagrange_interpolation(x[1], nodes, f_intermediate, barycentric_weights)
+
+  return point_value
+end
+
+
+end # @muladd

src/meshes/mesh_io.jl

@@ -157,7 +157,7 @@ function save_mesh_file(mesh::P4estMesh, output_directory, timestep=0)
 
   p4est_file = joinpath(output_directory, p4est_filename)
 
-  # Save the complete connectivity/p4est data to disk.
+  # Save the complete connectivity and `p4est` data to disk.
   save_p4est!(p4est_file, mesh.p4est)
 
   # Open file (clobber existing content)
@@ -264,7 +264,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
     boundary_names = boundary_names_ .|> Symbol
 
     p4est_file = joinpath(dirname(mesh_file), p4est_filename)
-    # Prevent Julia crashes when p4est can't find the file
+    # Prevent Julia crashes when `p4est` can't find the file
     @assert isfile(p4est_file)
 
     p4est = load_p4est(p4est_file, Val(ndims))

src/meshes/meshes.jl

@@ -9,6 +9,8 @@ include("tree_mesh.jl")
 include("structured_mesh.jl")
 include("surface_interpolant.jl")
 include("unstructured_mesh.jl")
+include("face_interpolant.jl")
+include("transfinite_mappings_3d.jl")
 include("p4est_mesh.jl")
 include("mesh_io.jl")
 include("dgmulti_meshes.jl")