P4estMesh constructor for HOHQMesh Abaqus-type file

docs/src/meshes/p4est_mesh.md

@@ -13,7 +13,7 @@ wrap numerical fluxes implemented only for Cartesian meshes. This simplifies
 the re-use of existing functionality for the [`TreeMesh`](@ref) but is usually
 less efficient, cf. [PR #550](https://github.com/trixi-framework/Trixi.jl/pull/550).
 
-## Construction of P4estMesh from an Abaqus file
+## Construction of 2D P4estMesh from an Abaqus file
 
 One available option to construct a [`P4estMesh`](@ref) is to read in an Abaqus (`.inp`) mesh file.
 We briefly describe the structure of this file, the conventions it uses, and how the mesh file
@@ -230,7 +230,7 @@ For completeness, we provide the entire Abaqus mesh file for the example mesh in
 
 ## Short discussion on 3D
 
-The Abaqus file format with high-order boundary information from `HOHQMesh` is very similar to the
+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:
 
 1. The element connectivity would be given in terms of the eight corners that define a hexahedron.

examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl

@@ -4,7 +4,7 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# semidiscretization of the linear advection equation
+# semidiscretization of the compressible Euler equations
 
 equations = CompressibleEulerEquations3D(1.4)
 
@@ -20,7 +20,7 @@ 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/4643ecb2fa452fe431268b85d8ab1a65426c5719/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
 

src/meshes/face_interpolant.jl

@@ -36,13 +36,14 @@ 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, fvals, wbary)
+function lagrange_interpolation_2d(x, nodes, function_values, barycentric_weights)
 
-  f_intermediate = zeros(eltype(fvals), length(nodes))
+  f_intermediate = zeros(eltype(function_values), length(nodes))
   for j in eachindex(nodes)
-    f_intermediate[j] = lagrange_interpolation(x[2], nodes, view(fvals, j, :), wbary)
+    f_intermediate[j] = lagrange_interpolation(x[2], nodes, view(function_values, j, :),
+                                               barycentric_weights)
   end
-  point_value = lagrange_interpolation(x[1], nodes, f_intermediate, wbary)
+  point_value = lagrange_interpolation(x[1], nodes, f_intermediate, barycentric_weights)
 
   return point_value
 end

src/meshes/transfinite_mappings_3d.jl

@@ -53,14 +53,14 @@ function straight_side_hex_map(xi, eta, zeta, corner_points)
 
   coordinate = zeros(eltype(xi), 3)
   for j in 1:3
-    coordinate[j] += (0.125 * ( corner_points[j, 1] * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta)
-                              + corner_points[j, 2] * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta)
-                              + corner_points[j, 3] * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta)
-                              + corner_points[j, 4] * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta)
-                              + corner_points[j, 5] * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta)
-                              + corner_points[j, 6] * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta)
-                              + corner_points[j, 7] * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta)
-                              + corner_points[j, 8] * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta) ) )
+    coordinate[j] += (0.125 * ( corner_points[j, 1] * (1 - xi) * (1 - eta) * (1 - zeta)
+                              + corner_points[j, 2] * (1 + xi) * (1 - eta) * (1 - zeta)
+                              + corner_points[j, 3] * (1 + xi) * (1 + eta) * (1 - zeta)
+                              + corner_points[j, 4] * (1 - xi) * (1 + eta) * (1 - zeta)
+                              + corner_points[j, 5] * (1 - xi) * (1 - eta) * (1 + zeta)
+                              + corner_points[j, 6] * (1 + xi) * (1 - eta) * (1 + zeta)
+                              + corner_points[j, 7] * (1 + xi) * (1 + eta) * (1 + zeta)
+                              + corner_points[j, 8] * (1 - xi) * (1 + eta) * (1 + zeta) ) )
   end
 
   return coordinate
@@ -93,20 +93,20 @@ function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:Curved
   corners = zeros(eltype(xi), (3, 8))
 
   # Compute values along the face edges
-  edge_values[:, 1] .= evaluate_at(SVector(xi  , -1.0), face_curves[1])
-  edge_values[:, 2] .= evaluate_at(SVector(1.0 , zeta), face_curves[1])
-  edge_values[:, 3] .= evaluate_at(SVector(xi  ,  1.0), face_curves[1])
-  edge_values[:, 4] .= evaluate_at(SVector(-1.0, zeta), face_curves[1])
+  edge_values[:, 1] .= evaluate_at(SVector(xi,   -1), face_curves[1])
+  edge_values[:, 2] .= evaluate_at(SVector( 1, zeta), face_curves[1])
+  edge_values[:, 3] .= evaluate_at(SVector(xi,    1), face_curves[1])
+  edge_values[:, 4] .= evaluate_at(SVector(-1, zeta), face_curves[1])
 
-  edge_values[:, 5] .= evaluate_at(SVector(xi  , -1.0), face_curves[2])
-  edge_values[:, 6] .= evaluate_at(SVector(1.0 , zeta), face_curves[2])
-  edge_values[:, 7] .= evaluate_at(SVector(xi  ,  1.0), face_curves[2])
-  edge_values[:, 8] .= evaluate_at(SVector(-1.0, zeta), face_curves[2])
+  edge_values[:, 5] .= evaluate_at(SVector(xi,   -1), face_curves[2])
+  edge_values[:, 6] .= evaluate_at(SVector( 1, zeta), face_curves[2])
+  edge_values[:, 7] .= evaluate_at(SVector(xi,    1), face_curves[2])
+  edge_values[:, 8] .= evaluate_at(SVector(-1, zeta), face_curves[2])
 
-  edge_values[:, 9]  .= evaluate_at(SVector(eta, -1.0), face_curves[6])
-  edge_values[:, 10] .= evaluate_at(SVector(eta, -1.0), face_curves[4])
-  edge_values[:, 11] .= evaluate_at(SVector(eta,  1.0), face_curves[4])
-  edge_values[:, 12] .= evaluate_at(SVector(eta,  1.0), face_curves[6])
+  edge_values[:, 9]  .= evaluate_at(SVector(eta, -1), face_curves[6])
+  edge_values[:, 10] .= evaluate_at(SVector(eta, -1), face_curves[4])
+  edge_values[:, 11] .= evaluate_at(SVector(eta,  1), face_curves[4])
+  edge_values[:, 12] .= evaluate_at(SVector(eta,  1), face_curves[6])
 
   # Compute values on the face
   face_values[:, 1] .= evaluate_at(SVector( xi, zeta), face_curves[1])
@@ -131,23 +131,23 @@ function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:Curved
   # Compute the transfinite mapping
   for j in 1:3
     # Linear interpolation between opposite faces
-    coordinate[j] = ( 0.5 * ( face_values[j, 6] * (1.0 - xi  ) + face_values[j, 4] * (1.0 + xi  )
-                            + face_values[j, 1] * (1.0 - eta ) + face_values[j, 2] * (1.0 + eta )
-                            + face_values[j, 3] * (1.0 - zeta) + face_values[j, 5] * (1.0 + zeta) ) )
+    coordinate[j] = ( 0.5 * ( face_values[j, 6] * (1 - xi  ) + face_values[j, 4] * (1 + xi  )
+                            + face_values[j, 1] * (1 - eta ) + face_values[j, 2] * (1 + eta )
+                            + face_values[j, 3] * (1 - zeta) + face_values[j, 5] * (1 + zeta) ) )
 
     # Edge corrections to ensure faces match
-    coordinate[j] -= ( 0.25 * ( edge_values[j, 1 ] * (1.0 - eta) * (1.0 - zeta)
-                              + edge_values[j, 2 ] * (1.0 + xi ) * (1.0 - eta )
-                              + edge_values[j, 3 ] * (1.0 - eta) * (1.0 + zeta)
-                              + edge_values[j, 4 ] * (1.0 - xi ) * (1.0 - eta )
-                              + edge_values[j, 5 ] * (1.0 + eta) * (1.0 - zeta)
-                              + edge_values[j, 6 ] * (1.0 + xi ) * (1.0 + eta )
-                              + edge_values[j, 7 ] * (1.0 + eta) * (1.0 + zeta)
-                              + edge_values[j, 8 ] * (1.0 - xi ) * (1.0 + eta )
-                              + edge_values[j, 9 ] * (1.0 - xi ) * (1.0 - zeta)
-                              + edge_values[j, 10] * (1.0 + xi ) * (1.0 - zeta)
-                              + edge_values[j, 11] * (1.0 + xi ) * (1.0 + zeta)
-                              + edge_values[j, 12] * (1.0 - xi ) * (1.0 + zeta) ) )
+    coordinate[j] -= ( 0.25 * ( edge_values[j, 1 ] * (1 - eta) * (1 - zeta)
+                              + edge_values[j, 2 ] * (1 + xi ) * (1 - eta )
+                              + edge_values[j, 3 ] * (1 - eta) * (1 + zeta)
+                              + edge_values[j, 4 ] * (1 - xi ) * (1 - eta )
+                              + edge_values[j, 5 ] * (1 + eta) * (1 - zeta)
+                              + edge_values[j, 6 ] * (1 + xi ) * (1 + eta )
+                              + edge_values[j, 7 ] * (1 + eta) * (1 + zeta)
+                              + edge_values[j, 8 ] * (1 - xi ) * (1 + eta )
+                              + edge_values[j, 9 ] * (1 - xi ) * (1 - zeta)
+                              + edge_values[j, 10] * (1 + xi ) * (1 - zeta)
+                              + edge_values[j, 11] * (1 + xi ) * (1 + zeta)
+                              + edge_values[j, 12] * (1 - xi ) * (1 + zeta) ) )
 
     # Subtracted interior twice, so add back the straight-sided hexahedral mapping
     coordinate[j] += coordinate_straight[j]

src/meshes/unstructured_mesh.jl

@@ -45,9 +45,9 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_
 
   # readin the number of nodes, number of interfaces, number of elements and local polynomial degree
   current_line = split(file_lines[2])
-  n_corners    = parse(Int, current_line[1])
-  n_surfaces   = parse(Int, current_line[2])
-  n_elements   = parse(Int, current_line[3])
+  n_corners = parse(Int, current_line[1])
+  n_surfaces = parse(Int, current_line[2])
+  n_elements = parse(Int, current_line[3])
   mesh_polydeg = parse(Int, current_line[4])
 
   mesh_nnodes = mesh_polydeg + 1
@@ -56,17 +56,17 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_
   # the mesh file. Thus, this cannot be type stable at all. Hence, we allocate
   # the memory now and introduce a function barrier before continuing to read
   # data from the file.
-  corner_nodes      = Array{RealT}(undef, (2, n_corners))
-  interface_info    = Array{Int}(undef, (6, n_surfaces))
-  element_node_ids  = Array{Int}(undef, (4, n_elements))
-  curved_check      = Vector{Int}(undef, 4)
-  cornerNodeVals    = Array{RealT}(undef, (4, 2))
-  tempNodes         = Array{RealT}(undef, (4, 2))
-  curve_vals        = Array{RealT}(undef, (mesh_nnodes, 2))
+  corner_nodes = Array{RealT}(undef, (2, n_corners))
+  interface_info = Array{Int}(undef, (6, n_surfaces))
+  element_node_ids = Array{Int}(undef, (4, n_elements))
+  curved_check = Vector{Int}(undef, 4)
+  quad_corners = Array{RealT}(undef, (4, 2))
+  quad_corners_flipped = Array{RealT}(undef, (4, 2))
+  curve_values = Array{RealT}(undef, (mesh_nnodes, 2))
   element_is_curved = Array{Bool}(undef, n_elements)
-  CurvedSurfaceT    = CurvedSurface{RealT}
-  surface_curves    = Array{CurvedSurfaceT}(undef, (4, n_elements))
-  boundary_names    = Array{Symbol}(undef, (4, n_elements))
+  CurvedSurfaceT = CurvedSurface{RealT}
+  surface_curves = Array{CurvedSurfaceT}(undef, (4, n_elements))
+  boundary_names = Array{Symbol}(undef, (4, n_elements))
 
   # create the Chebyshev-Gauss-Lobatto nodes used to represent any curved boundaries that are
   # required to construct the sides
@@ -76,7 +76,7 @@ function UnstructuredMesh2D(filename; RealT=Float64, periodicity=false, unsaved_
   bary_weights = SVector{mesh_nnodes}(bary_weights_)
 
   arrays = (; corner_nodes, interface_info, element_node_ids, curved_check,
-              cornerNodeVals, tempNodes, curve_vals,
+              quad_corners, quad_corners_flipped, curve_values,
               element_is_curved, surface_curves, boundary_names)
   counters = (; n_corners, n_surfaces, n_elements)
 
@@ -99,7 +99,7 @@ end
 
 function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights)
   @unpack ( corner_nodes, interface_info, element_node_ids, curved_check,
-            cornerNodeVals, tempNodes, curve_vals,
+            quad_corners, quad_corners_flipped, curve_values,
             element_is_curved, surface_curves, boundary_names ) = arrays
   @unpack n_corners, n_surfaces, n_elements = counters
   mesh_nnodes = length(cheby_nodes)
@@ -110,7 +110,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
   # readin an store the nodes that dictate the corners of the elements needed to construct the
   # element geometry terms
   for j in 1:n_corners
-    current_line       = split(file_lines[file_idx])
+    current_line = split(file_lines[file_idx])
     corner_nodes[1, j] = parse(RealT, current_line[1])
     corner_nodes[2, j] = parse(RealT, current_line[2])
     file_idx += 1
@@ -127,7 +127,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
   # container to for the interface neighbour information and connectivity
   n_boundaries = 0
   for j in 1:n_surfaces
-    current_line   = split(file_lines[file_idx])
+    current_line = split(file_lines[file_idx])
     interface_info[1, j] = parse(Int, current_line[1])
     interface_info[2, j] = parse(Int, current_line[2])
     interface_info[3, j] = parse(Int, current_line[3])
@@ -149,18 +149,18 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
 
   for j in 1:n_elements
     # pull the corner node IDs
-    current_line           = split(file_lines[file_idx])
+    current_line = split(file_lines[file_idx])
     element_node_ids[1, j] = parse(Int, current_line[1])
     element_node_ids[2, j] = parse(Int, current_line[2])
     element_node_ids[3, j] = parse(Int, current_line[3])
     element_node_ids[4, j] = parse(Int, current_line[4])
     for i in 1:4
       # pull the (x,y) values of these corners out of the nodes array
-      cornerNodeVals[i, :] .= corner_nodes[:, element_node_ids[i, j]]
+      quad_corners[i, :] .= corner_nodes[:, element_node_ids[i, j]]
     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])
+    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])
@@ -178,32 +178,30 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
       # 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, :]
+      @views quad_corners_flipped[1, :] .= quad_corners[4, :]
+      @views quad_corners_flipped[2, :] .= quad_corners[2, :]
+      @views quad_corners_flipped[3, :] .= quad_corners[3, :]
+      @views quad_corners_flipped[4, :] .= quad_corners[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
+          # when curved_check[i] is 0 then the "curve" from corner `i` to corner `i+1` is a
           # straight line. So we must construct the interpolant for this line
           for k in 1:mesh_nnodes
-            curve_vals[k, 1] = tempNodes[m1, 1] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 1]
-                                                                                 - tempNodes[m1, 1])
-            curve_vals[k, 2] = tempNodes[m1, 2] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 2]
-                                                                                 - tempNodes[m1, 2])
+            curve_values[k, 1] = linear_interpolate(cheby_nodes[k], quad_corners_flipped[m1, 1], quad_corners_flipped[m2, 1])
+            curve_values[k, 2] = linear_interpolate(cheby_nodes[k], quad_corners_flipped[m1, 2], quad_corners_flipped[m2, 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:mesh_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])
+            file_idx += 1
+            current_line = split(file_lines[file_idx])
+            curve_values[k, 1] = parse(RealT,current_line[1])
+            curve_values[k, 2] = parse(RealT,current_line[2])
           end
         end
         # construct the curve interpolant for the current side
-        surface_curves[i, j] = CurvedSurfaceT(cheby_nodes, bary_weights, copy(curve_vals))
+        surface_curves[i, j] = CurvedSurfaceT(cheby_nodes, bary_weights, copy(curve_values))
         # 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

test/test_p4est_3d.jl

@@ -123,8 +123,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_hohqmesh.jl"),
-      l2   = [0.004206309690620563, 0.004126097809243564, 0.004249654337625898, 0.0036399062057524295, 0.007401447142652835],
-      linf = [0.04530413742813599, 0.027657282796417437, 0.056698206703512, 0.04839054846610935, 0.11544778343499651])
+      l2   = [0.0042023406458005464, 0.004122532789279737, 0.0042448149597303616, 0.0036361316700401765, 0.007389845952982495],
+      linf = [0.04530610539892499, 0.02765695110527666, 0.05670295599308606, 0.048396544302230504, 0.1154589758186293])
   end
 end
 

test/test_unstructured_2d.jl

@@ -21,8 +21,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_euler_free_stream.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
-      l2   = [3.357431396258123e-14, 2.439943089578555e-13, 1.4655386790023588e-13, 4.670410488845425e-13],
-      linf = [1.0169198816356584e-11, 6.838458965763294e-11, 4.456759961080081e-11, 1.4207657272891083e-10],
+      l2   = [3.3937971107485363e-14, 2.447586447887882e-13, 1.4585205789296455e-13, 4.716993468962946e-13],
+      linf = [8.804734719092266e-12, 6.261270668606045e-11, 2.93670088247211e-11, 1.205400224080222e-10],
       tspan = (0.0, 0.1),
       atol = 3.0e-13)
   end
@@ -107,8 +107,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_shallowwater_well_balanced.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
-      l2   = [1.2164292510839083, 2.8049631232564753e-12, 1.5664980749498454e-12, 1.216429251083908],
-      linf = [1.5138512282315868, 5.0263880749562876e-11, 3.5805418805003075e-11, 1.513851228231574],
+      l2   = [1.2164292510839076, 2.6118925543469468e-12, 1.1636046671473883e-12, 1.2164292510839079],
+      linf = [1.5138512282315846, 4.998482888288039e-11, 2.0246214978154587e-11, 1.513851228231574],
       tspan = (0.0, 0.25))
   end
 
@@ -121,8 +121,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_shallowwater_dirichlet.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"),
-      l2   = [1.1577518608963063e-5, 4.371173779970303e-13, 4.2152984234036224e-13, 1.1577518608935235e-5],
-      linf = [8.394063878491842e-5, 9.632644747754261e-11, 9.54898901893391e-11, 8.394063879602065e-5],
+      l2   = [1.1577518608940115e-5, 4.867189932537344e-13, 4.647273240470541e-13, 1.1577518608933468e-5],
+      linf = [8.394063878602864e-5, 1.1469760027632646e-10, 1.1146619484429974e-10, 8.394063879602065e-5],
       tspan = (0.0, 2.0))
   end