From e0e9aa4baf081d789a85f7af7415cda7d0851b11 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Dec 2024 15:52:52 +0100 Subject: [PATCH 01/32] First steps towards quadratic/curved elements form standard Abaqus --- src/meshes/p4est_mesh.jl | 184 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 177 insertions(+), 7 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 65c0a431b29..b491ce952e7 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -452,10 +452,12 @@ function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_le connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_pw.num_trees[] + n_trees::Int = connectivity_pw.num_trees[] # = number elements n_vertices::Int = connectivity_pw.num_vertices[] # Extract a copy of the element vertices to compute the tree node coordinates + # `vertices` store coordinates of all three dimensions (even for the 2D case) + # since the Abaqus `.inp` format always stores 3D coordinates. vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) # Readin all the information from the mesh file into a string array @@ -497,6 +499,119 @@ function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_le return connectivity, tree_node_coordinates, nodes, boundary_names end +# p4est can handle only linear elements. This function checks the `meshfile` +# for quadratic elements (highest order supported by Standard Abaqus) and +# replaces them with linear elements. The higher-order (quadratic) boundaries are handled +# "internally" by Trixi as for the HOHQMesh-Abaqus case. +function preprocess_standard_abaqus_for_p4est(meshfile, dim, + linear_trusses, + linear_quads, + linear_hexes, + quadratic_trusses, + quadratic_quads, + quadratic_hexes) + meshfile_p4est_rdy = replace(meshfile, ".inp" => "_p4est_ready.inp") + order = 1 # Highest order of elements present in the mesh. + + elements_begin_index = 0 # Line number where the element section begins + open(meshfile, "r") do infile + open(meshfile_p4est_rdy, "w") do outfile + for (line_index, line) in enumerate(eachline(infile)) + # Copy header and node data to p4est_ready file + println(outfile, line) + if occursin("******* E L E M E N T S *************", line) + elements_begin_index = line_index + 1 + break + end + end + end + end + + current_element_type = "" + open(meshfile, "r") do infile + open(meshfile_p4est_rdy, "a") do outfile + for (line_index, line) in enumerate(eachline(infile)) + # Act only in the element section + if line_index >= elements_begin_index + # Check if a new element type/element set is defined + if startswith(line, "*") + # Retrieve element type + current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] + + # Check for linear elements - then we just copy the line + if occursin(linear_trusses, current_element_type) || + occursin(linear_quads, current_element_type) || + occursin(linear_hexes, current_element_type) + println(outfile, line) + else # Quadratic element - replace with linear + order = 2 + if occursin(quadratic_trusses, current_element_type) + # Distinction between 2D/3D should not make a difference here, + # Done for clarity & consistency here. + linear_truss_type = dim == 2 ? "T2D2" : "T3D2" + new_line = replace(line, + current_element_type => linear_truss_type) + elseif occursin(quadratic_quads, current_element_type) + linear_quad_type = "CPS4" + new_line = replace(line, + current_element_type => linear_quad_type) + elseif occursin(quadratic_hexes, current_element_type) + linear_hex_type = "C3D8" + new_line = replace(line, + current_element_type => linear_hex_type) + end + println(outfile, new_line) + end + else # Element data in line + # Check for linear elements - then we just copy the line + if occursin(linear_trusses, current_element_type) || + occursin(linear_quads, current_element_type) + println(outfile, line) + else + parts = split(line, ',') + if occursin(quadratic_trusses, current_element_type) + # Print the first (element), second (vertex 1), and fourth/last (vertex 2) indices to file + # For node order of quadratic trusses, check e.g. + # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch29s02ael13.html + new_line = join([parts[1], parts[2], parts[4]], ',') + elseif occursin(quadratic_quads, current_element_type) + # Print the first (element), second to fifth (vertices 1-4) indices to file + # For node order of quadratic quads, check e.g. + # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael02.html + new_line = join([ + parts[1], + parts[2], + parts[3], + parts[4], + parts[5] + ], ',') + elseif occursin(quadratic_hexes, current_element_type) + # Print the first (element), second to ninth (vertices 1-8) indices to file + # The node order is fortunately the same for hexes/bricks of type "C3D20", "C3D27", see + # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael03.html + new_line = join([ + parts[1], + parts[2], + parts[3], + parts[4], + parts[5], + parts[6], + parts[7], + parts[8], + parts[9] + ], ',') + end + println(outfile, new_line) + end + end # Line start if + end # Check if in element section of `.inp` file + end # Loop over lines in `.inp` files + end # Open output file in append mode + end # Open meshfile in read mode + + return meshfile_p4est_rdy, order +end + # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] # and a list of boundary names for the `P4estMesh`. The tree node coordinates are computed according to # the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. All boundary @@ -505,14 +620,54 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, initial_refinement_level, n_dimensions, RealT, boundary_symbols) + + ### Regular expressions for Abaqus element types ### + + # Linear trusses begin with T2D2 (2D) or T3D2 (3D), + # and may be followed by a bunch of Abaqus specifics ("E", "H", "T", etc.). + # These are, however, not relevant for the p4est mesh connectivity. + # For a full list of elements, see http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt10eli01.html + linear_trusses = r"^(T2D2|T3D2).*$" + + # These are the standard Abaqus linear quads. + # Note that there are many(!) more variants designed for specific purposes in the Abaqus solver. + # To keep it simple we support only basic quads, membranes, and shells here. + linear_quads = r"^(CPE4|CPEG4|CPS4|M3D4|S4|SFM3D4).*$" + + # Quadratic trusses begin with T2D3 (2D) or T3D3 (3D), + # and may be followed by a bunch of Abaqus specifics ("E", "H", "T", etc.). + # These are, however, not relevant for the p4est mesh connectivity. + quadratic_trusses = r"^(T2D3|T3D3).*$" + + # Same logic as for linear quads: + # Support only basic quads, membranes, and shells here. + quadratic_quads = r"^(CPE8|CPS8|CAX8|S8|M3D8|M3D9).*$" + + # In 3D only standard hexahedra/bricks are supported. + linear_hexes = r"^(C3D8).*$" + quadratic_hexes = r"^(C3D20|C3D27).*$" + + # Preprocess the meshfile to replace quadratic elements with linear elements + meshfile_p4est_rdy, mesh_polydeg = preprocess_standard_abaqus_for_p4est(meshfile, + n_dimensions, + linear_trusses, + linear_quads, + linear_hexes, + quadratic_trusses, + quadratic_quads, + quadratic_hexes) + # Create the mesh connectivity using `p4est` - connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) + connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_pw.num_trees[] + n_trees::Int = connectivity_pw.num_trees[] # = number elements n_vertices::Int = connectivity_pw.num_vertices[] + # Extract a copy of the element vertices to compute the tree node coordinates + # `vertices` store coordinates of all three dimensions (even for the 2D case) + # since the Abaqus `.inp` format always stores 3D coordinates. vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, (2^n_dimensions, n_trees)) @@ -524,14 +679,27 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, ntuple(_ -> length(nodes), n_dimensions)..., n_trees) - calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, - tree_to_vertex) + + # No nonlinearity in the mesh. Rely on the mapping function to realize the curvature. + if mesh_polydeg == 1 + calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, + tree_to_vertex) + else # mesh_polydeg = 2 => Nonlinearity in supplied mesh + mesh_nnodes = mesh_polydeg + 1 + + # Note: We ASSUME that the additional node between the end-vertices lies + # on the center on that line, such that we can use Chebyshev-Gauss-Lobatto nodes! + # For polydeg = 2, we have the 3 nodes [-1, 0, 1] (within the reference element). + cheby_nodes, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) + nodes = SVector{mesh_nnodes}(cheby_nodes) + end if boundary_symbols === nothing # There's no simple and generic way to distinguish boundaries without any information given. # Name all of them :all. boundary_names = fill(:all, 2 * n_dimensions, n_trees) else # Boundary information given + # TODO: Not sure which meshfile is required here - linear or quadratic? # 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 @@ -552,6 +720,7 @@ end function parse_elements(meshfile, n_trees, n_dims) @assert n_dims in (2, 3) "Only 2D and 3D meshes are supported" # Valid element types (that can be processed by p4est) based on dimension + # TODO: Make this more general: Use the linear element types defined above! element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : ["*ELEMENT, type=C3D8"] @@ -1492,7 +1661,7 @@ end # 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 + # Get the number of trees (elements) and the number of interpolation nodes n_trees = last(size(node_coordinates)) nnodes = length(nodes) @@ -1528,11 +1697,12 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, # Pull the (x,y) values of the four vertices of the current tree out of the global vertices array for i in 1:4 - quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] + quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] # 2D => 1:2 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]) + # Note: This strategy is HOHQMesh-Abaqus.inp specific! curved_check[1] = parse(Int, current_line[2]) curved_check[2] = parse(Int, current_line[3]) curved_check[3] = parse(Int, current_line[4]) From 101448ddb406b25bc4c7dc42b7f7cb02c8dc9bd1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Dec 2024 15:59:08 +0100 Subject: [PATCH 02/32] idea --- src/meshes/p4est_mesh.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index b491ce952e7..fdbd4062444 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -538,6 +538,7 @@ function preprocess_standard_abaqus_for_p4est(meshfile, dim, # Retrieve element type current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] + # TODO: Do I want to keep trusses at all? Same for quads in 3D. # Check for linear elements - then we just copy the line if occursin(linear_trusses, current_element_type) || occursin(linear_quads, current_element_type) || From 0b7d6468c8b902cd6bbdadf9e476d56c131dcb0c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Dec 2024 16:00:41 +0100 Subject: [PATCH 03/32] Notes --- src/meshes/p4est_mesh.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index fdbd4062444..adc1ca69683 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -539,6 +539,9 @@ function preprocess_standard_abaqus_for_p4est(meshfile, dim, current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] # TODO: Do I want to keep trusses at all? Same for quads in 3D. + # TODO: Should also probably start labeling the elements starting from 1 + # to have agreement with element & p4est tree index. + # Check for linear elements - then we just copy the line if occursin(linear_trusses, current_element_type) || occursin(linear_quads, current_element_type) || From 1e9af3070a6ff810b7009e8ba395ce20482f308e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Dec 2024 18:10:39 +0100 Subject: [PATCH 04/32] preproc abaqus --- src/meshes/p4est_mesh.jl | 197 ++++++++++++++++++++++++++------------- 1 file changed, 133 insertions(+), 64 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index adc1ca69683..f27edd90416 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -499,63 +499,140 @@ function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_le return connectivity, tree_node_coordinates, nodes, boundary_names end -# p4est can handle only linear elements. This function checks the `meshfile` -# for quadratic elements (highest order supported by Standard Abaqus) and -# replaces them with linear elements. The higher-order (quadratic) boundaries are handled -# "internally" by Trixi as for the HOHQMesh-Abaqus case. -function preprocess_standard_abaqus_for_p4est(meshfile, dim, - linear_trusses, - linear_quads, - linear_hexes, - quadratic_trusses, - quadratic_quads, - quadratic_hexes) - meshfile_p4est_rdy = replace(meshfile, ".inp" => "_p4est_ready.inp") - order = 1 # Highest order of elements present in the mesh. - - elements_begin_index = 0 # Line number where the element section begins +function preprocess_standard_abaqus(meshfile, + linear_quads, + linear_hexes, + quadratic_quads, + quadratic_hexes, + n_dimensions) + meshfile_preproc = replace(meshfile, ".inp" => "_preproc.inp") + + # Line number where the element section begins + elements_begin_idx = 0 + # Line number where the node and element sets begin (if present) + sets_begin_idx = 0 + + total_lines = 0 open(meshfile, "r") do infile - open(meshfile_p4est_rdy, "w") do outfile + # Copy header and node data to pre-processed file + open(meshfile_preproc, "w") do outfile for (line_index, line) in enumerate(eachline(infile)) - # Copy header and node data to p4est_ready file println(outfile, line) if occursin("******* E L E M E N T S *************", line) - elements_begin_index = line_index + 1 + elements_begin_idx = line_index + 1 break end + + total_lines += 1 end end + + # Find the line number where the node and element sets begin + for (line_index, line) in enumerate(eachline(infile)) + if line_index >= elements_begin_idx + if startswith(line, "*ELSET") || startswith(line, "*NSET") + sets_begin_idx = line_index + break + end + + total_lines += 1 + end + end + + # Catch the case where there are no node or element sets + if sets_begin_idx == 0 + sets_begin_idx = total_lines + 1 + end end + # Need to add `elements_begin_idx - 1` since the loop above starts at `elements_begin_idx` + sets_begin_idx += elements_begin_idx - 1 - current_element_type = "" open(meshfile, "r") do infile - open(meshfile_p4est_rdy, "a") do outfile + open(meshfile_preproc, "a") do outfile + print_following_lines = false + element_index = 1 + for (line_index, line) in enumerate(eachline(infile)) # Act only in the element section - if line_index >= elements_begin_index + if elements_begin_idx <= line_index < sets_begin_idx # Check if a new element type/element set is defined - if startswith(line, "*") + if startswith(line, "*ELEMENT") # Retrieve element type current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] - # TODO: Do I want to keep trusses at all? Same for quads in 3D. - # TODO: Should also probably start labeling the elements starting from 1 - # to have agreement with element & p4est tree index. + if n_dimensions == 2 && + (occursin(linear_quads, current_element_type) || + occursin(quadratic_quads, current_element_type)) + print_following_lines = true + println(outfile, line) + elseif n_dimensions == 3 && + (occursin(linear_hexes, current_element_type) || + occursin(quadratic_hexes, current_element_type)) + print_following_lines = true + println(outfile, line) + else + print_following_lines = false + end + else # Element data in line + # Check for linear elements - then we just copy the line + if print_following_lines + parts = split(line, ',') + parts[1] = string(element_index) + println(outfile, join(parts, ',')) + element_index += 1 + end + end + end + + # Print node and element sets as they are + if line_index >= sets_begin_idx + println(outfile, line) + end + end # Loop over lines in `.inp` files + end # Open output file in append mode + end # Open meshfile in read mode + + # TODO: Return + return meshfile_preproc, elements_begin_idx, sets_begin_idx +end + +# p4est can handle only linear elements. This function checks the `meshfile` +# for quadratic elements (highest order supported by Standard Abaqus) and +# replaces them with linear elements. The higher-order (quadratic) boundaries are handled +# "internally" by Trixi as for the HOHQMesh-Abaqus case. +function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, + linear_quads, + linear_hexes, + quadratic_quads, + quadratic_hexes, + elements_begin_idx, + sets_begin_idx) + meshfile_p4est_rdy = replace(meshfile_pre_proc, + "_preproc.inp" => "_p4est_ready.inp") + order = 1 # Highest order of elements present in the mesh. + + current_element_type = "" + new_line = "" + + open(meshfile_pre_proc, "r") do infile + open(meshfile_p4est_rdy, "w") do outfile + for (line_index, line) in enumerate(eachline(infile)) + # Copy header and node data + if line_index < elements_begin_idx || line_index >= sets_begin_idx + println(outfile, line) + else # Element data + # Check if a new element type/element set is defined + if startswith(line, "*ELEMENT") + # Retrieve element type + current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] # Check for linear elements - then we just copy the line - if occursin(linear_trusses, current_element_type) || - occursin(linear_quads, current_element_type) || + if occursin(linear_quads, current_element_type) || occursin(linear_hexes, current_element_type) println(outfile, line) else # Quadratic element - replace with linear order = 2 - if occursin(quadratic_trusses, current_element_type) - # Distinction between 2D/3D should not make a difference here, - # Done for clarity & consistency here. - linear_truss_type = dim == 2 ? "T2D2" : "T3D2" - new_line = replace(line, - current_element_type => linear_truss_type) - elseif occursin(quadratic_quads, current_element_type) + if occursin(quadratic_quads, current_element_type) linear_quad_type = "CPS4" new_line = replace(line, current_element_type => linear_quad_type) @@ -568,17 +645,12 @@ function preprocess_standard_abaqus_for_p4est(meshfile, dim, end else # Element data in line # Check for linear elements - then we just copy the line - if occursin(linear_trusses, current_element_type) || - occursin(linear_quads, current_element_type) + if occursin(linear_quads, current_element_type) || + occursin(linear_hexes, current_element_type) println(outfile, line) else parts = split(line, ',') - if occursin(quadratic_trusses, current_element_type) - # Print the first (element), second (vertex 1), and fourth/last (vertex 2) indices to file - # For node order of quadratic trusses, check e.g. - # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch29s02ael13.html - new_line = join([parts[1], parts[2], parts[4]], ',') - elseif occursin(quadratic_quads, current_element_type) + if occursin(quadratic_quads, current_element_type) # Print the first (element), second to fifth (vertices 1-4) indices to file # For node order of quadratic quads, check e.g. # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael02.html @@ -607,11 +679,11 @@ function preprocess_standard_abaqus_for_p4est(meshfile, dim, end println(outfile, new_line) end - end # Line start if - end # Check if in element section of `.inp` file - end # Loop over lines in `.inp` files - end # Open output file in append mode - end # Open meshfile in read mode + end + end + end + end + end return meshfile_p4est_rdy, order end @@ -627,22 +699,11 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, ### Regular expressions for Abaqus element types ### - # Linear trusses begin with T2D2 (2D) or T3D2 (3D), - # and may be followed by a bunch of Abaqus specifics ("E", "H", "T", etc.). - # These are, however, not relevant for the p4est mesh connectivity. - # For a full list of elements, see http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt10eli01.html - linear_trusses = r"^(T2D2|T3D2).*$" - # These are the standard Abaqus linear quads. # Note that there are many(!) more variants designed for specific purposes in the Abaqus solver. # To keep it simple we support only basic quads, membranes, and shells here. linear_quads = r"^(CPE4|CPEG4|CPS4|M3D4|S4|SFM3D4).*$" - # Quadratic trusses begin with T2D3 (2D) or T3D3 (3D), - # and may be followed by a bunch of Abaqus specifics ("E", "H", "T", etc.). - # These are, however, not relevant for the p4est mesh connectivity. - quadratic_trusses = r"^(T2D3|T3D3).*$" - # Same logic as for linear quads: # Support only basic quads, membranes, and shells here. quadratic_quads = r"^(CPE8|CPS8|CAX8|S8|M3D8|M3D9).*$" @@ -651,15 +712,22 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, linear_hexes = r"^(C3D8).*$" quadratic_hexes = r"^(C3D20|C3D27).*$" + meshfile_preproc, elements_begin_idx, sets_begin_idx = preprocess_standard_abaqus(meshfile, + linear_quads, + linear_hexes, + quadratic_quads, + quadratic_hexes, + n_dimensions) + # Preprocess the meshfile to replace quadratic elements with linear elements - meshfile_p4est_rdy, mesh_polydeg = preprocess_standard_abaqus_for_p4est(meshfile, - n_dimensions, - linear_trusses, + meshfile_p4est_rdy, mesh_polydeg = preprocess_standard_abaqus_for_p4est(meshfile_preproc, linear_quads, linear_hexes, - quadratic_trusses, quadratic_quads, - quadratic_hexes) + quadratic_hexes, + elements_begin_idx, + sets_begin_idx) + mesh_polydeg = 1 # TODO # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) @@ -703,7 +771,8 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # Name all of them :all. boundary_names = fill(:all, 2 * n_dimensions, n_trees) else # Boundary information given - # TODO: Not sure which meshfile is required here - linear or quadratic? + # TODO: Not sure which meshfile is required here! + # 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 From 49d79c4e76978b71e4c5bdbc9da06011b2d410ea Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Dec 2024 18:18:07 +0100 Subject: [PATCH 05/32] shorten --- src/meshes/p4est_mesh.jl | 38 ++++++++++---------------------------- 1 file changed, 10 insertions(+), 28 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index f27edd90416..c642c8b74ec 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -522,7 +522,6 @@ function preprocess_standard_abaqus(meshfile, elements_begin_idx = line_index + 1 break end - total_lines += 1 end end @@ -534,7 +533,6 @@ function preprocess_standard_abaqus(meshfile, sets_begin_idx = line_index break end - total_lines += 1 end end @@ -588,11 +586,10 @@ function preprocess_standard_abaqus(meshfile, if line_index >= sets_begin_idx println(outfile, line) end - end # Loop over lines in `.inp` files - end # Open output file in append mode - end # Open meshfile in read mode + end + end + end - # TODO: Return return meshfile_preproc, elements_begin_idx, sets_begin_idx end @@ -609,8 +606,9 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, sets_begin_idx) meshfile_p4est_rdy = replace(meshfile_pre_proc, "_preproc.inp" => "_p4est_ready.inp") - order = 1 # Highest order of elements present in the mesh. + order = 1 # Assume linear elements by default + # Some useful function-wide variables current_element_type = "" new_line = "" @@ -651,31 +649,15 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, else parts = split(line, ',') if occursin(quadratic_quads, current_element_type) - # Print the first (element), second to fifth (vertices 1-4) indices to file - # For node order of quadratic quads, check e.g. + # Print the first (element), second to fifth (vertices 1-4) indices to file. + # For node order of quadratic (secod-order) quads, see # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael02.html - new_line = join([ - parts[1], - parts[2], - parts[3], - parts[4], - parts[5] - ], ',') + new_line = join(parts[1:5], ',') elseif occursin(quadratic_hexes, current_element_type) - # Print the first (element), second to ninth (vertices 1-8) indices to file + # Print the first (element), second to ninth (vertices 1-8) indices to file. # The node order is fortunately the same for hexes/bricks of type "C3D20", "C3D27", see # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael03.html - new_line = join([ - parts[1], - parts[2], - parts[3], - parts[4], - parts[5], - parts[6], - parts[7], - parts[8], - parts[9] - ], ',') + new_line = join(parts[1:9], ',') end println(outfile, new_line) end From 2c597a418381bd5a2d9554869d3cc9f6fbed9bb4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 13 Dec 2024 11:39:42 +0100 Subject: [PATCH 06/32] work on quad mesches --- src/meshes/p4est_mesh.jl | 217 +++++++++++++++++++++++++++++++-------- 1 file changed, 177 insertions(+), 40 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index c642c8b74ec..5f841925784 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -512,7 +512,6 @@ function preprocess_standard_abaqus(meshfile, # Line number where the node and element sets begin (if present) sets_begin_idx = 0 - total_lines = 0 open(meshfile, "r") do infile # Copy header and node data to pre-processed file open(meshfile_preproc, "w") do outfile @@ -522,7 +521,6 @@ function preprocess_standard_abaqus(meshfile, elements_begin_idx = line_index + 1 break end - total_lines += 1 end end @@ -533,22 +531,22 @@ function preprocess_standard_abaqus(meshfile, sets_begin_idx = line_index break end - total_lines += 1 end end - - # Catch the case where there are no node or element sets - if sets_begin_idx == 0 - sets_begin_idx = total_lines + 1 - end end - # Need to add `elements_begin_idx - 1` since the loop above starts at `elements_begin_idx` - sets_begin_idx += elements_begin_idx - 1 + # Catch the case where there are no node or element sets + if sets_begin_idx == 0 + sets_begin_idx = length(readlines(meshfile)) + 1 + else + # Need to add `elements_begin_idx - 1` since the loop above starts at `elements_begin_idx` + sets_begin_idx += elements_begin_idx - 1 + end + element_index = 0 + preproc_element_section_lines = 0 open(meshfile, "r") do infile open(meshfile_preproc, "a") do outfile print_following_lines = false - element_index = 1 for (line_index, line) in enumerate(eachline(infile)) # Act only in the element section @@ -563,21 +561,24 @@ function preprocess_standard_abaqus(meshfile, occursin(quadratic_quads, current_element_type)) print_following_lines = true println(outfile, line) + preproc_element_section_lines += 1 elseif n_dimensions == 3 && (occursin(linear_hexes, current_element_type) || occursin(quadratic_hexes, current_element_type)) print_following_lines = true println(outfile, line) + preproc_element_section_lines += 1 else print_following_lines = false end else # Element data in line - # Check for linear elements - then we just copy the line if print_following_lines + element_index += 1 parts = split(line, ',') + # Exchange element index parts[1] = string(element_index) println(outfile, join(parts, ',')) - element_index += 1 + preproc_element_section_lines += 1 end end end @@ -589,6 +590,8 @@ function preprocess_standard_abaqus(meshfile, end end end + # Adjust `sets_begin_idx` to correct line number after removing unnecessary elements + sets_begin_idx = elements_begin_idx + preproc_element_section_lines return meshfile_preproc, elements_begin_idx, sets_begin_idx end @@ -670,6 +673,34 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, return meshfile_p4est_rdy, order end +function read_nodes_standard_abaqus(meshfile, n_dimensions, + elements_begin_idx, RealT) + mesh_nodes = Dict{Int, SVector{n_dimensions, RealT}}() + + node_section = false + open(meshfile, "r") do infile + for (line_index, line) in enumerate(eachline(infile)) + if line_index < elements_begin_idx + if startswith(line, "*NODE") + node_section = true + end + if node_section + parts = split(line, ',') + if length(parts) >= n_dimensions + 1 + node_id = parse(Int, parts[1]) + coordinates = SVector{n_dimensions, RealT}([parse(RealT, + parts[i]) + for i in 2:(n_dimensions + 1)]) + mesh_nodes[node_id] = coordinates + end + end + end + end + end + + return mesh_nodes +end + # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] # and a list of boundary names for the `P4estMesh`. The tree node coordinates are computed according to # the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. All boundary @@ -709,7 +740,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - mesh_polydeg = 1 # TODO + #mesh_polydeg = 1 # TODO # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) @@ -723,10 +754,11 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # `vertices` store coordinates of all three dimensions (even for the 2D case) # since the Abaqus `.inp` format always stores 3D coordinates. vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) + tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, (2^n_dimensions, n_trees)) - basis = LobattoLegendreBasis(RealT, polydeg) + basis = LobattoLegendreBasis(RealT, mesh_polydeg) nodes = basis.nodes tree_node_coordinates = Array{RealT, n_dimensions + 2}(undef, n_dimensions, @@ -744,8 +776,25 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # Note: We ASSUME that the additional node between the end-vertices lies # on the center on that line, such that we can use Chebyshev-Gauss-Lobatto nodes! # For polydeg = 2, we have the 3 nodes [-1, 0, 1] (within the reference element). - cheby_nodes, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) - nodes = SVector{mesh_nnodes}(cheby_nodes) + cheby_nodes_, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) + cheby_nodes = SVector{mesh_nnodes}(cheby_nodes_) + + println("cheby_nodes: ", cheby_nodes) + + mesh_nodes = read_nodes_standard_abaqus(meshfile_preproc, n_dimensions, + elements_begin_idx, RealT) + + element_lines = readlines(open(meshfile_preproc))[elements_begin_idx:(sets_begin_idx - 1)] + + if n_dimensions == 2 + calc_tree_node_coordinates!(tree_node_coordinates, element_lines, cheby_nodes, + vertices, RealT, + linear_quads, mesh_nodes) + else # n_dimensions == 3 + calc_tree_node_coordinates!(tree_node_coordinates, element_lines, cheby_nodes, + vertices, RealT, + linear_hexes, mesh_nodes) + end end if boundary_symbols === nothing @@ -756,9 +805,10 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # TODO: Not sure which meshfile is required here! # Read in nodes belonging to boundaries - node_set_dict = parse_node_sets(meshfile, boundary_symbols) + node_set_dict = parse_node_sets(meshfile_p4est_rdy, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries - element_node_matrix = parse_elements(meshfile, n_trees, n_dimensions) + element_node_matrix = parse_elements(meshfile_p4est_rdy, n_trees, n_dimensions, + elements_begin_idx, sets_begin_idx) # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) @@ -772,13 +822,9 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, return connectivity, tree_node_coordinates, nodes, boundary_names end -function parse_elements(meshfile, n_trees, n_dims) - @assert n_dims in (2, 3) "Only 2D and 3D meshes are supported" - # Valid element types (that can be processed by p4est) based on dimension - # TODO: Make this more general: Use the linear element types defined above! - element_types = n_dims == 2 ? - ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", - "*ELEMENT, type=S4"] : ["*ELEMENT, type=C3D8"] +function parse_elements(meshfile, n_trees, n_dims, + elements_begin_idx, + sets_begin_idx) # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 @@ -786,20 +832,22 @@ function parse_elements(meshfile, n_trees, n_dims) 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. - # Important: Do not use index from the Abaqus file, but the one from p4est. - element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id - tree_id += 1 - else # Processed all elements for this ELSET - el_list_follows = false + open(meshfile, "r") do infile + for (line_index, line) in enumerate(eachline(infile)) + if elements_begin_idx <= line_index < sets_begin_idx + if startswith(line, "*ELEMENT") + 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. + # Important: Do not use index from the Abaqus file, but the one from p4est. + element_node_matrix[tree_id, :] = content_int[2:end] # First entry is element id + tree_id += 1 + else # Processed all elements for this `ELSET` + el_list_follows = false + end end end end @@ -1822,6 +1870,95 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, return file_idx end +function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, + element_lines, nodes, vertices, RealT, + linear_quads, mesh_nodes) + nnodes = length(nodes) + + # 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) + quad_vertices = Array{RealT}(undef, (4, 2)) + curve_values = 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_) + + element_set_order = 1 + tree = 0 + for line_idx in 1:length(element_lines) + line = element_lines[line_idx] + + # Check if a new element type/element set is defined + if startswith(line, "*ELEMENT") + # Retrieve element type + current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] + + # Check if these are linear elements + if occursin(linear_quads, current_element_type) + element_set_order = 1 + else + element_set_order = 2 + end + else # Element data + tree += 1 + + # Pull the vertex node IDs + line_split = split(line, r",\s+") + element_nodes = parse.(Int, line_split) + + if element_set_order == 1 + element_node_ids[1] = element_nodes[2] + element_node_ids[2] = element_nodes[3] + element_node_ids[3] = element_nodes[4] + element_node_ids[4] = element_nodes[5] + + # Create the node coordinates on this particular element + # Pull the (x,y) values of the four vertices of the current tree out of the global vertices array + for i in 1:4 + quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] # 2D => 1:2 + end + + # CARE: Not tested if this works or if we need to interpolate + # linear elements to have full higher-order elements + calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) + else # element_set_order == 2 + for edge in 1:4 + node1 = element_nodes[edge + 1] # "Left" node + node2 = element_nodes[edge + 5] # "Middle" node + if edge < 4 + node3 = element_nodes[edge + 2] # "Right" node + else + node3 = element_nodes[2] # "Right" node + end + + node1_coords = mesh_nodes[node1] + node2_coords = mesh_nodes[node2] + node3_coords = mesh_nodes[node3] + + curve_values[1, 1] = node1_coords[1] + curve_values[1, 2] = node1_coords[2] + curve_values[2, 1] = node2_coords[1] + curve_values[2, 2] = node2_coords[2] + curve_values[3, 1] = node3_coords[1] + curve_values[3, 2] = node3_coords[2] + + # Construct the curve interpolant for the current side + surface_curves[edge] = CurvedSurfaceT(nodes, bary_weights, + copy(curve_values)) + end + + # Create the node coordinates on this particular element + calc_node_coordinates!(node_coordinates, tree, nodes, surface_curves) + end + end + end +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 `transfinite_mappings_3d.jl` From 70d46ca37a922da5990ffeb2ba80f0264d963afd Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 13 Dec 2024 14:33:26 +0100 Subject: [PATCH 07/32] Seemingly working version 2nd order Abaqus 2D --- src/meshes/p4est_mesh.jl | 103 ++++++++++++++++++++++++--------------- 1 file changed, 64 insertions(+), 39 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 5f841925784..6b7bcffbbfa 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -526,11 +526,9 @@ function preprocess_standard_abaqus(meshfile, # Find the line number where the node and element sets begin for (line_index, line) in enumerate(eachline(infile)) - if line_index >= elements_begin_idx - if startswith(line, "*ELSET") || startswith(line, "*NSET") - sets_begin_idx = line_index - break - end + if startswith(line, "*ELSET") || startswith(line, "*NSET") + sets_begin_idx = line_index + break end end end @@ -740,7 +738,6 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - #mesh_polydeg = 1 # TODO # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) @@ -758,9 +755,18 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, (2^n_dimensions, n_trees)) - basis = LobattoLegendreBasis(RealT, mesh_polydeg) + basis = LobattoLegendreBasis(RealT, polydeg) nodes = basis.nodes + if mesh_polydeg == 2 + mesh_nnodes = mesh_polydeg + 1 + # Note: We ASSUME that the additional node between the end-vertices lies + # on the center on that line, such that we can use Chebyshev-Gauss-Lobatto nodes! + # For polydeg = 2, we have the 3 nodes [-1, 0, 1] (within the reference element). + cheby_nodes, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) + nodes = SVector{mesh_nnodes}(cheby_nodes) + end + tree_node_coordinates = Array{RealT, n_dimensions + 2}(undef, n_dimensions, ntuple(_ -> length(nodes), n_dimensions)..., @@ -768,31 +774,21 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # No nonlinearity in the mesh. Rely on the mapping function to realize the curvature. if mesh_polydeg == 1 - calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, vertices, - tree_to_vertex) + calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, + vertices, tree_to_vertex) else # mesh_polydeg = 2 => Nonlinearity in supplied mesh - mesh_nnodes = mesh_polydeg + 1 - - # Note: We ASSUME that the additional node between the end-vertices lies - # on the center on that line, such that we can use Chebyshev-Gauss-Lobatto nodes! - # For polydeg = 2, we have the 3 nodes [-1, 0, 1] (within the reference element). - cheby_nodes_, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes) - cheby_nodes = SVector{mesh_nnodes}(cheby_nodes_) - - println("cheby_nodes: ", cheby_nodes) - mesh_nodes = read_nodes_standard_abaqus(meshfile_preproc, n_dimensions, elements_begin_idx, RealT) element_lines = readlines(open(meshfile_preproc))[elements_begin_idx:(sets_begin_idx - 1)] if n_dimensions == 2 - calc_tree_node_coordinates!(tree_node_coordinates, element_lines, cheby_nodes, - vertices, RealT, + calc_tree_node_coordinates!(tree_node_coordinates, element_lines, + nodes, vertices, RealT, linear_quads, mesh_nodes) else # n_dimensions == 3 - calc_tree_node_coordinates!(tree_node_coordinates, element_lines, cheby_nodes, - vertices, RealT, + calc_tree_node_coordinates!(tree_node_coordinates, element_lines, + nodes, vertices, RealT, linear_hexes, mesh_nodes) end end @@ -802,8 +798,6 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # Name all of them :all. boundary_names = fill(:all, 2 * n_dimensions, n_trees) else # Boundary information given - # TODO: Not sure which meshfile is required here! - # Read in nodes belonging to boundaries node_set_dict = parse_node_sets(meshfile_p4est_rdy, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries @@ -1922,30 +1916,61 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, for i in 1:4 quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] # 2D => 1:2 end - - # CARE: Not tested if this works or if we need to interpolate - # linear elements to have full higher-order elements calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) else # element_set_order == 2 for edge in 1:4 node1 = element_nodes[edge + 1] # "Left" node node2 = element_nodes[edge + 5] # "Middle" node - if edge < 4 - node3 = element_nodes[edge + 2] # "Right" node - else - node3 = element_nodes[2] # "Right" node - end + node3 = edge == 4 ? element_nodes[2] : element_nodes[edge + 2] # "Right" node node1_coords = mesh_nodes[node1] node2_coords = mesh_nodes[node2] node3_coords = mesh_nodes[node3] - curve_values[1, 1] = node1_coords[1] - curve_values[1, 2] = node1_coords[2] - curve_values[2, 1] = node2_coords[1] - curve_values[2, 2] = node2_coords[2] - curve_values[3, 1] = node3_coords[1] - curve_values[3, 2] = node3_coords[2] + # The nodes for an Abaqus element are labeled following a closed path + # around the element: + # + # <---- + # *-----*-----* + # | | + # | | | ^ + # | * * | + # v | | | + # | | + # *-----*-----* + # ----> + # `node_coordinates`, however, requires to sort the nodes into a + # unique coordinate system, + # + # *-----*-----* + # | | + # | | + # * * + # | | + # ^η | | + # | *-----*-----* + # |----> ξ + # thus we need to flip the nodes for the second xi and eta edges met. + + if edge in [1, 2] + curve_values[1, 1] = node1_coords[1] + curve_values[1, 2] = node1_coords[2] + + curve_values[2, 1] = node2_coords[1] + curve_values[2, 2] = node2_coords[2] + + curve_values[3, 1] = node3_coords[1] + curve_values[3, 2] = node3_coords[2] + else # Flip "left" and "right" nodes + curve_values[1, 1] = node3_coords[1] + curve_values[1, 2] = node3_coords[2] + + curve_values[2, 1] = node2_coords[1] + curve_values[2, 2] = node2_coords[2] + + curve_values[3, 1] = node1_coords[1] + curve_values[3, 2] = node1_coords[2] + end # Construct the curve interpolant for the current side surface_curves[edge] = CurvedSurfaceT(nodes, bary_weights, From 7838c24c8cbb28365496e389a0144fbede864263 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 13 Dec 2024 14:50:58 +0100 Subject: [PATCH 08/32] Comments --- src/meshes/p4est_mesh.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 6b7bcffbbfa..1d388c3aa6a 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -554,6 +554,8 @@ function preprocess_standard_abaqus(meshfile, # Retrieve element type current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] + # Keep only quads (2D) or hexes (3D), i.e., eliminate all other elements + # like trusses (2D/3D) or quads in 3D. if n_dimensions == 2 && (occursin(linear_quads, current_element_type) || occursin(quadratic_quads, current_element_type)) @@ -671,6 +673,7 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, return meshfile_p4est_rdy, order end +# Read all nodes (not only vertices, i.e., endpoints of elements) into a dict function read_nodes_standard_abaqus(meshfile, n_dimensions, elements_begin_idx, RealT) mesh_nodes = Dict{Int, SVector{n_dimensions, RealT}}() @@ -1947,7 +1950,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, # | | # * * # | | - # ^η | | + # ^ η | | # | *-----*-----* # |----> ξ # thus we need to flip the nodes for the second xi and eta edges met. From b85e41b21ea1077d84e4b89e69c386e5a058c4a6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 14 Dec 2024 23:09:00 +0100 Subject: [PATCH 09/32] experiment with viz --- .../elixir_advection_amr_unstructured_flag.jl | 19 +++++++++++++------ src/meshes/p4est_mesh.jl | 2 +- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index 4bfb2d3e375..067340809a2 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -8,13 +8,14 @@ using Trixi advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) -initial_condition = initial_condition_gauss +initial_condition = initial_condition_constant boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = Dict(:all => boundary_condition) -solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) +solver = DGSEM(polydeg = 2, surface_flux = flux_lax_friedrichs) +#= # Deformed rectangle that looks like a waving flag, # lower and upper faces are sinus curves, left and right are vertical lines. f1(s) = SVector(-5.0, 5 * s - 5.0) @@ -36,6 +37,10 @@ mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2e mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping_flag, initial_refinement_level = 1) +=# + +#mesh = P4estMesh{2}("./sd7003_laminar.inp") +mesh = P4estMesh{2}("./Test_Curved_Abaqus.inp") semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) @@ -43,7 +48,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 10.0) +tspan = (0.0, 0.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -57,7 +62,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_restart = SaveRestartCallback(interval = 100, save_final_restart = true) -save_solution = SaveSolutionCallback(interval = 100, +save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) @@ -75,8 +80,10 @@ stepsize_callback = StepsizeCallback(cfl = 0.7) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_restart, save_solution, - amr_callback, stepsize_callback); + #save_restart, + save_solution, + #amr_callback, + stepsize_callback); ############################################################################### # run the simulation diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 1d388c3aa6a..bfe60ae23d2 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -762,7 +762,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, nodes = basis.nodes if mesh_polydeg == 2 - mesh_nnodes = mesh_polydeg + 1 + mesh_nnodes = mesh_polydeg + 1 # = 3 # Note: We ASSUME that the additional node between the end-vertices lies # on the center on that line, such that we can use Chebyshev-Gauss-Lobatto nodes! # For polydeg = 2, we have the 3 nodes [-1, 0, 1] (within the reference element). From 82c506a82cb75036b13aa492bbf7fb76b5a2eee7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Dec 2024 15:08:54 +0100 Subject: [PATCH 10/32] test hohqmesh --- .../p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index 067340809a2..aabef736b46 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -40,7 +40,10 @@ mesh = P4estMesh{2}(mesh_file, polydeg = 3, =# #mesh = P4estMesh{2}("./sd7003_laminar.inp") -mesh = P4estMesh{2}("./Test_Curved_Abaqus.inp") + +#mesh = P4estMesh{2}("./Test_Curved_Abaqus.inp") + +mesh = P4estMesh{2}("./Test_Abaqus_HOHQMesh_Style.inp") semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) From c1f0d6c3db514ffecb93e0b9f63e615f98e9ae45 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 16 Dec 2024 09:33:22 +0100 Subject: [PATCH 11/32] revert --- .../elixir_advection_amr_unstructured_flag.jl | 22 +++++-------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index aabef736b46..4bfb2d3e375 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -8,14 +8,13 @@ using Trixi advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) -initial_condition = initial_condition_constant +initial_condition = initial_condition_gauss boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = Dict(:all => boundary_condition) -solver = DGSEM(polydeg = 2, surface_flux = flux_lax_friedrichs) +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) -#= # Deformed rectangle that looks like a waving flag, # lower and upper faces are sinus curves, left and right are vertical lines. f1(s) = SVector(-5.0, 5 * s - 5.0) @@ -37,13 +36,6 @@ mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2e mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping_flag, initial_refinement_level = 1) -=# - -#mesh = P4estMesh{2}("./sd7003_laminar.inp") - -#mesh = P4estMesh{2}("./Test_Curved_Abaqus.inp") - -mesh = P4estMesh{2}("./Test_Abaqus_HOHQMesh_Style.inp") semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) @@ -51,7 +43,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 0.0) +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -65,7 +57,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_restart = SaveRestartCallback(interval = 100, save_final_restart = true) -save_solution = SaveSolutionCallback(interval = 1000, +save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) @@ -83,10 +75,8 @@ stepsize_callback = StepsizeCallback(cfl = 0.7) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - #save_restart, - save_solution, - #amr_callback, - stepsize_callback); + save_restart, save_solution, + amr_callback, stepsize_callback); ############################################################################### # run the simulation From 3bf13639070703d5624869d93d9d042d8c275879 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 16 Dec 2024 09:41:07 +0100 Subject: [PATCH 12/32] comments --- src/meshes/p4est_mesh.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index bfe60ae23d2..4802238f0c9 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -499,6 +499,9 @@ function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_le return connectivity, tree_node_coordinates, nodes, boundary_names end +# This removes irrelevant elements from the meshfile, i.e., trusses/beams which are 1D elements. +# Thus, if `n_dimensions == 2`, only quads are kept, and if `n_dimensions == 3`, +# only hexes are kept, i.e., in that case also quads are removed. function preprocess_standard_abaqus(meshfile, linear_quads, linear_hexes, @@ -673,7 +676,8 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, return meshfile_p4est_rdy, order end -# Read all nodes (not only vertices, i.e., endpoints of elements) into a dict +# Read all nodes (not only vertices, i.e., endpoints of elements) into a dict. +# Those are required to enable higher-order boundaries for quadratic elements. function read_nodes_standard_abaqus(meshfile, n_dimensions, elements_begin_idx, RealT) mesh_nodes = Dict{Int, SVector{n_dimensions, RealT}}() @@ -1867,6 +1871,8 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, return file_idx end +# Version for quadratic 2D elements, i.e., second-order quads. +# TODO: 3D version function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element_lines, nodes, vertices, RealT, linear_quads, mesh_nodes) @@ -1942,8 +1948,8 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, # | | # *-----*-----* # ----> - # `node_coordinates`, however, requires to sort the nodes into a - # unique coordinate system, + # `curve_values`, however, requires to sort the nodes into a + # valid coordinate system, # # *-----*-----* # | | @@ -1953,7 +1959,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, # ^ η | | # | *-----*-----* # |----> ξ - # thus we need to flip the nodes for the second xi and eta edges met. + # thus we need to flip the node order for the second xi and eta edges met. if edge in [1, 2] curve_values[1, 1] = node1_coords[1] From fde8369cb8187f24796de962128083f4bc6cbca4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 16 Dec 2024 17:47:53 +0100 Subject: [PATCH 13/32] SD7003 --- .../elixir_navierstokes_SD70003airfoil.jl | 151 ++++++++++++++++++ 1 file changed, 151 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl new file mode 100644 index 00000000000..f1d94d30258 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -0,0 +1,151 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +U_inf = 0.2 +c_inf = 1.0 + +rho_inf = 1.4 # with gamma = 1.4 => p_inf = 1.0 + +Re = 10000.0 +airfoil_cord_length = 1.0 + +t_c = airfoil_cord_length / U_inf + +aoa = 4 * pi / 180 +u_x = U_inf * cos(aoa) +u_y = U_inf * sin(aoa) + +gamma = 1.4 +prandtl_number() = 0.72 +mu() = rho_inf * U_inf * airfoil_cord_length / Re + +equations = CompressibleEulerEquations2D(gamma) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), + Prandtl = prandtl_number(), + gradient_variables = GradientVariablesPrimitive()) + +@inline function initial_condition_mach02_flow(x, t, equations) + # set the freestream flow parameters + rho_freestream = 1.4 + + v1 = 0.19951281005196486 # 0.2 * cos(aoa) + v2 = 0.01395129474882506 # 0.2 * sin(aoa) + + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach02_flow + +# Boundary conditions for free-stream testing +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) + +velocity_bc_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_airfoil = BoundaryConditionNavierStokesWall(velocity_bc_airfoil, heat_bc) + +polydeg = 3 + +surf_flux = flux_hllc +vol_flux = flux_chandrashekar +solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, + volume_integral = VolumeIntegralFluxDifferencing(vol_flux)) + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +#path = "/storage/home/daniel/PERK4/SD7003/" + +path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Laminar/" + +#mesh_file = path * "sd7003_laminar_straight_sided_Trixi.inp" + +mesh_file = path * "sd7003_laminar.inp" + +boundary_symbols = [:Airfoil, :FarField] + +mesh = P4estMesh{2}(mesh_file, boundary_symbols = boundary_symbols, + initial_refinement_level = 0) + +boundary_conditions = Dict(:FarField => boundary_condition_free_stream, + :Airfoil => boundary_condition_slip_wall) + +boundary_conditions_parabolic = Dict(:FarField => boundary_condition_free_stream, + :Airfoil => boundary_condition_airfoil) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.005 * t_c) # Try to get into a state where initial pressure wave is gone +tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +# Choose analysis interval such that roughly every dt_c = 0.005 a record is taken +analysis_interval = 25 # PERK4_Multi, PERKSingle + +f_aoa() = aoa +f_rho_inf() = rho_inf +f_U_inf() = U_inf +f_linf() = airfoil_cord_length + +drag_coefficient = AnalysisSurfaceIntegral((:Airfoil,), + DragCoefficientPressure(f_aoa(), f_rho_inf(), + f_U_inf(), f_linf())) + +drag_coefficient_shear_force = AnalysisSurfaceIntegral((:Airfoil,), + DragCoefficientShearStress(f_aoa(), + f_rho_inf(), + f_U_inf(), + f_linf())) + +lift_coefficient = AnalysisSurfaceIntegral((:Airfoil,), + LiftCoefficientPressure(f_aoa(), f_rho_inf(), + f_U_inf(), f_linf())) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_errors = Symbol[], + analysis_integrals = (drag_coefficient, + drag_coefficient_shear_force, + lift_coefficient)) + +stepsize_callback = StepsizeCallback(cfl = 2.0) # PERK_4 Multi E = 5, ..., 14 + +# For plots etc +save_solution = SaveSolutionCallback(interval = 1_000_000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim, + output_directory = "out") + +alive_callback = AliveCallback(alive_interval = 10) + +save_restart = SaveRestartCallback(interval = Int(10^7), # Only at end + save_final_restart = true) + +callbacks = CallbackSet(stepsize_callback, # For measurements: Fixed timestep (do not use this) + alive_callback, # Not needed for measurement run + save_solution, # For plotting during measurement run + #save_restart, # For restart with measurements + summary_callback); + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); + dt = 1.0, save_everystep = false, callback = callbacks) + +summary_callback() # print the timer summary From 5592a391a6d2af514dab9cf5e326c38cb8c569fc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 17 Dec 2024 17:30:11 +0100 Subject: [PATCH 14/32] start 3d --- Project.toml | 4 + .../elixir_navierstokes_SD70003airfoil.jl | 3 +- .../elixir_navierstokes_SD70003airfoil.jl | 164 ++++++++++++++++++ src/meshes/p4est_mesh.jl | 147 +++++++++++++++- 4 files changed, 314 insertions(+), 4 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl diff --git a/Project.toml b/Project.toml index 296883f5420..e200188003e 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -24,11 +25,13 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -51,6 +54,7 @@ TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [weakdeps] Convex = "f65535da-76fb-5f13-bab9-19810c17039a" diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl index f1d94d30258..b93a24101f6 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -64,8 +64,7 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Laminar/" #mesh_file = path * "sd7003_laminar_straight_sided_Trixi.inp" - -mesh_file = path * "sd7003_laminar.inp" +#mesh_file = path * "sd7003_laminar.inp" boundary_symbols = [:Airfoil, :FarField] diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl new file mode 100644 index 00000000000..7a1d9491e42 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -0,0 +1,164 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +U_inf = 0.2 +c_inf = 1.0 + +rho_inf = 1.4 # with gamma = 1.4 => p_inf = 1.0 + +Re = 10000.0 +wall_cord_length = 1.0 + +t_c = wall_cord_length / U_inf + +aoa = 4 * pi / 180 +u_x = U_inf * cos(aoa) +u_y = U_inf * sin(aoa) + +gamma = 1.4 +prandtl_number() = 0.72 +mu() = rho_inf * U_inf * wall_cord_length / Re + +equations = CompressibleEulerEquations3D(gamma) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), + Prandtl = prandtl_number(), + gradient_variables = GradientVariablesPrimitive()) + +@inline function initial_condition_mach02_flow(x, t, equations) + # set the freestream flow parameters + rho_freestream = 1.4 + + v1 = 0.19951281005196486 # 0.2 * cos(aoa) + v2 = 0.01395129474882506 # 0.2 * sin(aoa) + v3 = 0.0 + + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, v3, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach02_flow + +# Boundary conditions for free-stream testing +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) + +velocity_bc_wall = NoSlip((x, t, equations) -> SVector(0.0, 0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_wall = BoundaryConditionNavierStokesWall(velocity_bc_wall, heat_bc) + +polydeg = 3 + +surf_flux = flux_hllc +vol_flux = flux_chandrashekar +solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, + volume_integral = VolumeIntegralFluxDifferencing(vol_flux)) + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +#path = "/storage/home/daniel/PERK4/SD7003/" + +path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Turbulent/" +mesh_file = path * "sd7003_turbulent.inp" + +boundary_symbols = [:wall, :rieminv] + +mesh = P4estMesh{3}(mesh_file, boundary_symbols = boundary_symbols) + +boundary_conditions = Dict(:rieminv => boundary_condition_free_stream, + :wall => boundary_condition_slip_wall) + +velocity_bc = NoSlip() do x, t, equations_parabolic + Trixi.velocity(initial_condition_mach02_flow(x, + t, + equations_parabolic), + equations_parabolic) +end + +heat_bc = Isothermal() do x, t, equations_parabolic + Trixi.temperature(initial_condition_mach02_flow(x, + t, + equations_parabolic), + equations_parabolic) +end + +boundary_condition_rieminv = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) + +boundary_conditions_parabolic = Dict(:rieminv => boundary_condition_rieminv, + :wall => boundary_condition_wall) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.0 * t_c) # Try to get into a state where initial pressure wave is gone +#tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +# Choose analysis interval such that roughly every dt_c = 0.005 a record is taken +analysis_interval = 25 # PERK4_Multi, PERKSingle + +f_aoa() = aoa +f_rho_inf() = rho_inf +f_U_inf() = U_inf +f_linf() = wall_cord_length + +drag_coefficient = AnalysisSurfaceIntegral((:wall,), + DragCoefficientPressure(f_aoa(), f_rho_inf(), + f_U_inf(), f_linf())) + +drag_coefficient_shear_force = AnalysisSurfaceIntegral((:wall,), + DragCoefficientShearStress(f_aoa(), + f_rho_inf(), + f_U_inf(), + f_linf())) + +lift_coefficient = AnalysisSurfaceIntegral((:wall,), + LiftCoefficientPressure(f_aoa(), f_rho_inf(), + f_U_inf(), f_linf())) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_errors = Symbol[], + analysis_integrals = (drag_coefficient, + drag_coefficient_shear_force, + lift_coefficient)) + +stepsize_callback = StepsizeCallback(cfl = 2.0) # PERK_4 Multi E = 5, ..., 14 + +# For plots etc +save_solution = SaveSolutionCallback(interval = 1_000_000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim, + output_directory = "out") + +alive_callback = AliveCallback(alive_interval = 10) + +save_restart = SaveRestartCallback(interval = Int(10^7), # Only at end + save_final_restart = true) + +callbacks = CallbackSet(stepsize_callback, # For measurements: Fixed timestep (do not use this) + alive_callback, # Not needed for measurement run + save_solution, # For plotting during measurement run + #save_restart, # For restart with measurements + summary_callback); + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); + dt = 1.0, save_everystep = false, callback = callbacks) + +summary_callback() # print the timer summary diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 4802238f0c9..da4183ca2bc 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -745,7 +745,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - + mesh_polydeg = 1 # TODO: Only for testing # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) connectivity_pw = PointerWrapper(connectivity) @@ -1872,7 +1872,6 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, end # Version for quadratic 2D elements, i.e., second-order quads. -# TODO: 3D version function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element_lines, nodes, vertices, RealT, linear_quads, mesh_nodes) @@ -2090,6 +2089,150 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, return file_idx end +# TODO: 3D version +function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, + element_lines, nodes, vertices, RealT, + linear_hexes, mesh_nodes) + nnodes = length(nodes) + + # Create a work set of Gamma curves to create the node coordinates + CurvedFaceT = CurvedFace{RealT} + face_curves = Array{CurvedFaceT}(undef, 6) + + # Create other work arrays to perform the mesh construction + element_node_ids = Array{Int}(undef, 8) + hex_vertices = Array{RealT}(undef, (3, 8)) + curve_values = Array{RealT}(undef, (3, nnodes, nnodes)) + + # Create the barycentric weights used for the surface interpolations + bary_weights_ = barycentric_weights(nodes) + bary_weights = SVector{nnodes}(bary_weights_) + + element_set_order = 1 + tree = 0 + for line_idx in 1:length(element_lines) + line = element_lines[line_idx] + + # Check if a new element type/element set is defined + if startswith(line, "*ELEMENT") + # Retrieve element type + current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] + + # Check if these are linear elements + if occursin(linear_hexes, current_element_type) + element_set_order = 1 + else + element_set_order = 2 + end + else # Element data + tree += 1 + + # Pull the vertex node IDs + line_split = split(line, r",\s+") + element_nodes = parse.(Int, line_split) + + if element_set_order == 1 + element_node_ids[1] = element_nodes[2] + element_node_ids[2] = element_nodes[3] + element_node_ids[3] = element_nodes[4] + element_node_ids[4] = element_nodes[5] + element_node_ids[5] = element_nodes[6] + element_node_ids[6] = element_nodes[7] + element_node_ids[7] = element_nodes[8] + element_node_ids[8] = element_nodes[9] + + # Create the node coordinates on this particular element + # Pull the (x,y, z) values of the four vertices of the current tree out of the global vertices array + for i in 1:8 + hex_vertices[i, :] .= vertices[:, element_node_ids[i]] # 3D => 1:3 = : + end + calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) + else # element_set_order == 2 + for face in 1:6 + # The second-order face has the following node distribution: + # NW N NE + # *-----*-----* + # | | + # | | + # W * * E + # | | + # | | + # *-----*-----* + # SW S SE + + # TODO: Care for special cases! + node1 = element_nodes[face + 1] # "SW" node + node2 = element_nodes[face + 8] # "S" node + node3 = element_nodes[face + 2] # "SE" node + node4 = element_nodes[face + 18] # "E" node + node5 = element_nodes[face + 6] # "NE" node + node6 = element_nodes[face + 12] # "N" node + node7 = element_nodes[face + 5] # "NW" node + node8 = element_nodes[face + 17] # "W" node + + #node3 = edge == 4 ? element_nodes[2] : element_nodes[edge + 2] # "Right" node + + node1_coords = mesh_nodes[node1] + node2_coords = mesh_nodes[node2] + node3_coords = mesh_nodes[node3] + + # The nodes for an Abaqus element are labeled following a closed path + # around the element: + # + # <---- + # *-----*-----* + # | | + # | | | ^ + # | * * | + # v | | | + # | | + # *-----*-----* + # ----> + # `curve_values`, however, requires to sort the nodes into a + # valid coordinate system, + # + # *-----*-----* + # | | + # | | + # * * + # | | + # ^ η | | + # | *-----*-----* + # |----> ξ + # thus we need to flip the node order for the second xi and eta edges met. + + if edge in [1, 2] + curve_values[1, 1] = node1_coords[1] + curve_values[1, 2] = node1_coords[2] + + curve_values[2, 1] = node2_coords[1] + curve_values[2, 2] = node2_coords[2] + + curve_values[3, 1] = node3_coords[1] + curve_values[3, 2] = node3_coords[2] + else # Flip "left" and "right" nodes + curve_values[1, 1] = node3_coords[1] + curve_values[1, 2] = node3_coords[2] + + curve_values[2, 1] = node2_coords[1] + curve_values[2, 2] = node2_coords[2] + + curve_values[3, 1] = node1_coords[1] + curve_values[3, 2] = node1_coords[2] + end + + # Construct the curve interpolant for the current side + surface_curves[edge] = CurvedSurfaceT(nodes, bary_weights, + copy(curve_values)) + end + + # Create the node coordinates on this particular element + calc_node_coordinates!(node_coordinates, tree, nodes, surface_curves) + end + end + end +end + # Given the eight `hex_vertices` for a hexahedral element extract # the four `face_vertices` for a particular `face_index`. function get_vertices_for_bilinear_interpolant!(face_vertices, face_index, hex_vertices) From 1c51559c9edce96ee4bb2371b82b73c0518cd7f9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Dec 2024 10:16:53 +0100 Subject: [PATCH 15/32] continue --- .../elixir_navierstokes_SD70003airfoil.jl | 2 +- src/meshes/p4est_mesh.jl | 122 +++++++++++------- 2 files changed, 76 insertions(+), 48 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl index 7a1d9491e42..0fd658d476d 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -63,7 +63,7 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, #path = "/storage/home/daniel/PERK4/SD7003/" path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Turbulent/" -mesh_file = path * "sd7003_turbulent.inp" +mesh_file = path * "sd7003_turbulent_fix_lines.inp" boundary_symbols = [:wall, :rieminv] diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index da4183ca2bc..1b6aca7802c 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -728,7 +728,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, # In 3D only standard hexahedra/bricks are supported. linear_hexes = r"^(C3D8).*$" - quadratic_hexes = r"^(C3D20|C3D27).*$" + quadratic_hexes = r"^(C3D27).*$" meshfile_preproc, elements_begin_idx, sets_begin_idx = preprocess_standard_abaqus(meshfile, linear_quads, @@ -745,7 +745,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - mesh_polydeg = 1 # TODO: Only for testing + # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) connectivity_pw = PointerWrapper(connectivity) @@ -2103,6 +2103,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, element_node_ids = Array{Int}(undef, 8) hex_vertices = Array{RealT}(undef, (3, 8)) curve_values = Array{RealT}(undef, (3, nnodes, nnodes)) + face_nodes = Array{Int}(undef, 9) # Create the barycentric weights used for the surface interpolations bary_weights_ = barycentric_weights(nodes) @@ -2154,80 +2155,107 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, # *-----*-----* # | | # | | - # W * * E + # W * *C * E # | | # | | # *-----*-----* # SW S SE - # TODO: Care for special cases! - node1 = element_nodes[face + 1] # "SW" node - node2 = element_nodes[face + 8] # "S" node - node3 = element_nodes[face + 2] # "SE" node - node4 = element_nodes[face + 18] # "E" node - node5 = element_nodes[face + 6] # "NE" node - node6 = element_nodes[face + 12] # "N" node - node7 = element_nodes[face + 5] # "NW" node - node8 = element_nodes[face + 17] # "W" node + face_nodes[1] = element_nodes[face + 1] # "SW" node + face_nodes[2] = element_nodes[face + 9] # "S" node + face_nodes[3] = element_nodes[face + 2] # "SE" node + face_nodes[4] = element_nodes[face + 18] # "E" node + face_nodes[5] = element_nodes[face + 6] # "NE" node + face_nodes[6] = element_nodes[face + 13] # "N" node + face_nodes[7] = element_nodes[face + 5] # "NW" node + face_nodes[8] = element_nodes[face + 17] # "W" node + + if face < 4 + face_nodes[9] = element_nodes[face + 24] # "C" node + elseif face == 4 + face_nodes[3] = element_nodes[2] + face_nodes[4] = element_nodes[18] + face_nodes[5] = element_nodes[6] + elseif face == 5 + face_nodes[1] = element_nodes[2] + face_nodes[2] = element_nodes[10] + face_nodes[3] = element_nodes[3] + face_nodes[4] = element_nodes[11] + face_nodes[5] = element_nodes[4] + face_nodes[6] = element_nodes[12] + face_nodes[7] = element_nodes[5] + face_nodes[8] = element_nodes[13] + face_nodes[9] = element_nodes[23] + elseif face == 6 + face_nodes[1] = element_nodes[6] + face_nodes[2] = element_nodes[14] + face_nodes[3] = element_nodes[7] + face_nodes[4] = element_nodes[15] + face_nodes[5] = element_nodes[8] + face_nodes[6] = element_nodes[16] + face_nodes[7] = element_nodes[9] + face_nodes[8] = element_nodes[17] + face_nodes[9] = element_nodes[24] + end - #node3 = edge == 4 ? element_nodes[2] : element_nodes[edge + 2] # "Right" node - - node1_coords = mesh_nodes[node1] - node2_coords = mesh_nodes[node2] - node3_coords = mesh_nodes[node3] + node1_coords = mesh_nodes[face_nodes[1]] + node2_coords = mesh_nodes[face_nodes[2]] + node3_coords = mesh_nodes[face_nodes[3]] + node4_coords = mesh_nodes[face_nodes[4]] + node5_coords = mesh_nodes[face_nodes[5]] + node6_coords = mesh_nodes[face_nodes[6]] + node7_coords = mesh_nodes[face_nodes[7]] + node8_coords = mesh_nodes[face_nodes[8]] + node9_coords = mesh_nodes[face_nodes[9]] # The nodes for an Abaqus element are labeled following a closed path # around the element: # - # <---- - # *-----*-----* - # | | - # | | | ^ - # | * * | - # v | | | - # | | - # *-----*-----* - # ----> + # <---- + # 7 6 5 + # *-----*-----* + # | | + # | | | ^ + # | 8* 9* *4 | + # v | | | + # | | + # *-----*-----* + # 1 2 3 + # ----> # `curve_values`, however, requires to sort the nodes into a # valid coordinate system, # # *-----*-----* # | | # | | - # * * + # * * * # | | # ^ η | | # | *-----*-----* # |----> ξ # thus we need to flip the node order for the second xi and eta edges met. - if edge in [1, 2] - curve_values[1, 1] = node1_coords[1] - curve_values[1, 2] = node1_coords[2] - - curve_values[2, 1] = node2_coords[1] - curve_values[2, 2] = node2_coords[2] - - curve_values[3, 1] = node3_coords[1] - curve_values[3, 2] = node3_coords[2] - else # Flip "left" and "right" nodes - curve_values[1, 1] = node3_coords[1] - curve_values[1, 2] = node3_coords[2] + # Proceed along bottom edge + curve_values[:, 1, 1] = node1_coords + curve_values[:, 2, 1] = node2_coords + curve_values[:, 3, 1] = node3_coords - curve_values[2, 1] = node2_coords[1] - curve_values[2, 2] = node2_coords[2] + # Proceed along middle line + curve_values[:, 1, 2] = node8_coords + curve_values[:, 2, 2] = node9_coords + curve_values[:, 3, 2] = node4_coords - curve_values[3, 1] = node1_coords[1] - curve_values[3, 2] = node1_coords[2] - end + # Proceed along top edge + curve_values[:, 1, 3] = node7_coords + curve_values[:, 2, 3] = node6_coords + curve_values[:, 3, 3] = node5_coords # Construct the curve interpolant for the current side - surface_curves[edge] = CurvedSurfaceT(nodes, bary_weights, - copy(curve_values)) + face_curves[face] = CurvedFaceT(nodes, bary_weights, copy(curve_values)) end # Create the node coordinates on this particular element - calc_node_coordinates!(node_coordinates, tree, nodes, surface_curves) + calc_node_coordinates!(node_coordinates, tree, nodes, face_curves) end end end From 1ddc717591393fda1d1252df9814e6f2f6c8a18d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Dec 2024 13:45:04 +0100 Subject: [PATCH 16/32] Continue --- .../elixir_navierstokes_SD70003airfoil.jl | 2 +- src/meshes/p4est_mesh.jl | 23 ++++++++++++------- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl index 0fd658d476d..c00fd14371e 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -63,7 +63,7 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, #path = "/storage/home/daniel/PERK4/SD7003/" path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Turbulent/" -mesh_file = path * "sd7003_turbulent_fix_lines.inp" +mesh_file = path * "sd7003_reduced.inp" boundary_symbols = [:wall, :rieminv] diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 1b6aca7802c..6b57fd091c0 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -745,7 +745,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - + #mesh_polydeg = 1 # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) connectivity_pw = PointerWrapper(connectivity) @@ -2145,7 +2145,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, # Create the node coordinates on this particular element # Pull the (x,y, z) values of the four vertices of the current tree out of the global vertices array for i in 1:8 - hex_vertices[i, :] .= vertices[:, element_node_ids[i]] # 3D => 1:3 = : + hex_vertices[:, i] .= vertices[:, element_node_ids[i]] # 3D => 1:3 = : end calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) else # element_set_order == 2 @@ -2170,9 +2170,11 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, face_nodes[7] = element_nodes[face + 5] # "NW" node face_nodes[8] = element_nodes[face + 17] # "W" node - if face < 4 - face_nodes[9] = element_nodes[face + 24] # "C" node - elseif face == 4 + if face <= 4 + face_nodes[9] = element_nodes[face + 24] # "C" node + end + + if face == 4 face_nodes[3] = element_nodes[2] face_nodes[4] = element_nodes[18] face_nodes[5] = element_nodes[6] @@ -2235,19 +2237,24 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, # |----> ξ # thus we need to flip the node order for the second xi and eta edges met. + # TODO: Need probably to do some flips for the different faces, as for the 2D case # Proceed along bottom edge curve_values[:, 1, 1] = node1_coords - curve_values[:, 2, 1] = node2_coords + #curve_values[:, 2, 1] = node2_coords + curve_values[:, 2, 1] = 0.5 .* (node1_coords .+ node3_coords) + curve_values[:, 3, 1] = node3_coords # Proceed along middle line curve_values[:, 1, 2] = node8_coords - curve_values[:, 2, 2] = node9_coords + #curve_values[:, 2, 2] = node9_coords + curve_values[:, 2, 2] = 0.5 .* (node8_coords .+ node4_coords) curve_values[:, 3, 2] = node4_coords # Proceed along top edge curve_values[:, 1, 3] = node7_coords - curve_values[:, 2, 3] = node6_coords + #curve_values[:, 2, 3] = node6_coords + curve_values[:, 2, 3] = 0.5 .* (node7_coords .+ node5_coords) curve_values[:, 3, 3] = node5_coords # Construct the curve interpolant for the current side From 0bbafa01549cb6cc9c07efb8bf827c25e251f691 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Dec 2024 14:27:22 +0100 Subject: [PATCH 17/32] Reduce to 2D --- .../elixir_navierstokes_SD70003airfoil.jl | 164 -------------- src/meshes/p4est_mesh.jl | 202 +----------------- src/meshes/t8code_mesh.jl | 2 - 3 files changed, 11 insertions(+), 357 deletions(-) delete mode 100644 examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl deleted file mode 100644 index c00fd14371e..00000000000 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ /dev/null @@ -1,164 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the compressible Euler equations - -U_inf = 0.2 -c_inf = 1.0 - -rho_inf = 1.4 # with gamma = 1.4 => p_inf = 1.0 - -Re = 10000.0 -wall_cord_length = 1.0 - -t_c = wall_cord_length / U_inf - -aoa = 4 * pi / 180 -u_x = U_inf * cos(aoa) -u_y = U_inf * sin(aoa) - -gamma = 1.4 -prandtl_number() = 0.72 -mu() = rho_inf * U_inf * wall_cord_length / Re - -equations = CompressibleEulerEquations3D(gamma) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), - Prandtl = prandtl_number(), - gradient_variables = GradientVariablesPrimitive()) - -@inline function initial_condition_mach02_flow(x, t, equations) - # set the freestream flow parameters - rho_freestream = 1.4 - - v1 = 0.19951281005196486 # 0.2 * cos(aoa) - v2 = 0.01395129474882506 # 0.2 * sin(aoa) - v3 = 0.0 - - p_freestream = 1.0 - - prim = SVector(rho_freestream, v1, v2, v3, p_freestream) - return prim2cons(prim, equations) -end - -initial_condition = initial_condition_mach02_flow - -# Boundary conditions for free-stream testing -boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) - -velocity_bc_wall = NoSlip((x, t, equations) -> SVector(0.0, 0.0, 0.0)) -heat_bc = Adiabatic((x, t, equations) -> 0.0) -boundary_condition_wall = BoundaryConditionNavierStokesWall(velocity_bc_wall, heat_bc) - -polydeg = 3 - -surf_flux = flux_hllc -vol_flux = flux_chandrashekar -solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, - volume_integral = VolumeIntegralFluxDifferencing(vol_flux)) - -############################################################################### -# Get the uncurved mesh from a file (downloads the file if not available locally) - -#path = "/storage/home/daniel/PERK4/SD7003/" - -path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Turbulent/" -mesh_file = path * "sd7003_reduced.inp" - -boundary_symbols = [:wall, :rieminv] - -mesh = P4estMesh{3}(mesh_file, boundary_symbols = boundary_symbols) - -boundary_conditions = Dict(:rieminv => boundary_condition_free_stream, - :wall => boundary_condition_slip_wall) - -velocity_bc = NoSlip() do x, t, equations_parabolic - Trixi.velocity(initial_condition_mach02_flow(x, - t, - equations_parabolic), - equations_parabolic) -end - -heat_bc = Isothermal() do x, t, equations_parabolic - Trixi.temperature(initial_condition_mach02_flow(x, - t, - equations_parabolic), - equations_parabolic) -end - -boundary_condition_rieminv = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc) - -boundary_conditions_parabolic = Dict(:rieminv => boundary_condition_rieminv, - :wall => boundary_condition_wall) - -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, solver; - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 0.0 * t_c) # Try to get into a state where initial pressure wave is gone -#tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() - -# Choose analysis interval such that roughly every dt_c = 0.005 a record is taken -analysis_interval = 25 # PERK4_Multi, PERKSingle - -f_aoa() = aoa -f_rho_inf() = rho_inf -f_U_inf() = U_inf -f_linf() = wall_cord_length - -drag_coefficient = AnalysisSurfaceIntegral((:wall,), - DragCoefficientPressure(f_aoa(), f_rho_inf(), - f_U_inf(), f_linf())) - -drag_coefficient_shear_force = AnalysisSurfaceIntegral((:wall,), - DragCoefficientShearStress(f_aoa(), - f_rho_inf(), - f_U_inf(), - f_linf())) - -lift_coefficient = AnalysisSurfaceIntegral((:wall,), - LiftCoefficientPressure(f_aoa(), f_rho_inf(), - f_U_inf(), f_linf())) - -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - output_directory = "out", - save_analysis = true, - analysis_errors = Symbol[], - analysis_integrals = (drag_coefficient, - drag_coefficient_shear_force, - lift_coefficient)) - -stepsize_callback = StepsizeCallback(cfl = 2.0) # PERK_4 Multi E = 5, ..., 14 - -# For plots etc -save_solution = SaveSolutionCallback(interval = 1_000_000, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim, - output_directory = "out") - -alive_callback = AliveCallback(alive_interval = 10) - -save_restart = SaveRestartCallback(interval = Int(10^7), # Only at end - save_final_restart = true) - -callbacks = CallbackSet(stepsize_callback, # For measurements: Fixed timestep (do not use this) - alive_callback, # Not needed for measurement run - save_solution, # For plotting during measurement run - #save_restart, # For restart with measurements - summary_callback); - -############################################################################### -# run the simulation - -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); - dt = 1.0, save_everystep = false, callback = callbacks) - -summary_callback() # print the timer summary diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 6b57fd091c0..41a492e1ce2 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -414,7 +414,6 @@ end function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level, n_dimensions, RealT) connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_hohqmesh_abaqus(meshfile, - initial_refinement_level, n_dimensions, RealT) @@ -431,7 +430,6 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, - initial_refinement_level, n_dimensions, RealT, boundary_symbols) @@ -445,7 +443,7 @@ end # and a list of boundary names for the `P4estMesh`. High-order boundary curve information as well as # the boundary names on each tree are provided by the `meshfile` created by # [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl). -function p4est_connectivity_from_hohqmesh_abaqus(meshfile, initial_refinement_level, +function p4est_connectivity_from_hohqmesh_abaqus(meshfile, n_dimensions, RealT) # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) @@ -707,18 +705,20 @@ function read_nodes_standard_abaqus(meshfile, n_dimensions, end # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] -# and a list of boundary names for the `P4estMesh`. The tree node coordinates are computed according to -# the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. All boundary -# names are given the name `:all`. +# and a list of boundary names for the `P4estMesh`. For linear meshes, the tree node coordinates are +# computed according to the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. +# For quadratic (second-order) meshes, the tree node coordinates are read from the file, similar as for +# `p4est_connectivity_from_hohqmesh_abaqus`. function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, - initial_refinement_level, n_dimensions, + n_dimensions, RealT, boundary_symbols) ### Regular expressions for Abaqus element types ### # These are the standard Abaqus linear quads. - # Note that there are many(!) more variants designed for specific purposes in the Abaqus solver. + # Note that there are many(!) more variants designed for specific purposes in the Abaqus solver, see + # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt10eli01.html # To keep it simple we support only basic quads, membranes, and shells here. linear_quads = r"^(CPE4|CPEG4|CPS4|M3D4|S4|SFM3D4).*$" @@ -745,7 +745,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, elements_begin_idx, sets_begin_idx) - #mesh_polydeg = 1 + # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile_p4est_rdy, Val(n_dimensions)) connectivity_pw = PointerWrapper(connectivity) @@ -794,9 +794,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, nodes, vertices, RealT, linear_quads, mesh_nodes) else # n_dimensions == 3 - calc_tree_node_coordinates!(tree_node_coordinates, element_lines, - nodes, vertices, RealT, - linear_hexes, mesh_nodes) + @error "Three dimensional quadratic elements are not supported yet!" end end @@ -1871,6 +1869,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, return file_idx end +# TODO: 3D version # Version for quadratic 2D elements, i.e., second-order quads. function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element_lines, nodes, vertices, RealT, @@ -2089,185 +2088,6 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, return file_idx end -# TODO: 3D version -function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, - element_lines, nodes, vertices, RealT, - linear_hexes, mesh_nodes) - nnodes = length(nodes) - - # Create a work set of Gamma curves to create the node coordinates - CurvedFaceT = CurvedFace{RealT} - face_curves = Array{CurvedFaceT}(undef, 6) - - # Create other work arrays to perform the mesh construction - element_node_ids = Array{Int}(undef, 8) - hex_vertices = Array{RealT}(undef, (3, 8)) - curve_values = Array{RealT}(undef, (3, nnodes, nnodes)) - face_nodes = Array{Int}(undef, 9) - - # Create the barycentric weights used for the surface interpolations - bary_weights_ = barycentric_weights(nodes) - bary_weights = SVector{nnodes}(bary_weights_) - - element_set_order = 1 - tree = 0 - for line_idx in 1:length(element_lines) - line = element_lines[line_idx] - - # Check if a new element type/element set is defined - if startswith(line, "*ELEMENT") - # Retrieve element type - current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] - - # Check if these are linear elements - if occursin(linear_hexes, current_element_type) - element_set_order = 1 - else - element_set_order = 2 - end - else # Element data - tree += 1 - - # Pull the vertex node IDs - line_split = split(line, r",\s+") - element_nodes = parse.(Int, line_split) - - if element_set_order == 1 - element_node_ids[1] = element_nodes[2] - element_node_ids[2] = element_nodes[3] - element_node_ids[3] = element_nodes[4] - element_node_ids[4] = element_nodes[5] - element_node_ids[5] = element_nodes[6] - element_node_ids[6] = element_nodes[7] - element_node_ids[7] = element_nodes[8] - element_node_ids[8] = element_nodes[9] - - # Create the node coordinates on this particular element - # Pull the (x,y, z) values of the four vertices of the current tree out of the global vertices array - for i in 1:8 - hex_vertices[:, i] .= vertices[:, element_node_ids[i]] # 3D => 1:3 = : - end - calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) - else # element_set_order == 2 - for face in 1:6 - # The second-order face has the following node distribution: - # NW N NE - # *-----*-----* - # | | - # | | - # W * *C * E - # | | - # | | - # *-----*-----* - # SW S SE - - face_nodes[1] = element_nodes[face + 1] # "SW" node - face_nodes[2] = element_nodes[face + 9] # "S" node - face_nodes[3] = element_nodes[face + 2] # "SE" node - face_nodes[4] = element_nodes[face + 18] # "E" node - face_nodes[5] = element_nodes[face + 6] # "NE" node - face_nodes[6] = element_nodes[face + 13] # "N" node - face_nodes[7] = element_nodes[face + 5] # "NW" node - face_nodes[8] = element_nodes[face + 17] # "W" node - - if face <= 4 - face_nodes[9] = element_nodes[face + 24] # "C" node - end - - if face == 4 - face_nodes[3] = element_nodes[2] - face_nodes[4] = element_nodes[18] - face_nodes[5] = element_nodes[6] - elseif face == 5 - face_nodes[1] = element_nodes[2] - face_nodes[2] = element_nodes[10] - face_nodes[3] = element_nodes[3] - face_nodes[4] = element_nodes[11] - face_nodes[5] = element_nodes[4] - face_nodes[6] = element_nodes[12] - face_nodes[7] = element_nodes[5] - face_nodes[8] = element_nodes[13] - face_nodes[9] = element_nodes[23] - elseif face == 6 - face_nodes[1] = element_nodes[6] - face_nodes[2] = element_nodes[14] - face_nodes[3] = element_nodes[7] - face_nodes[4] = element_nodes[15] - face_nodes[5] = element_nodes[8] - face_nodes[6] = element_nodes[16] - face_nodes[7] = element_nodes[9] - face_nodes[8] = element_nodes[17] - face_nodes[9] = element_nodes[24] - end - - node1_coords = mesh_nodes[face_nodes[1]] - node2_coords = mesh_nodes[face_nodes[2]] - node3_coords = mesh_nodes[face_nodes[3]] - node4_coords = mesh_nodes[face_nodes[4]] - node5_coords = mesh_nodes[face_nodes[5]] - node6_coords = mesh_nodes[face_nodes[6]] - node7_coords = mesh_nodes[face_nodes[7]] - node8_coords = mesh_nodes[face_nodes[8]] - node9_coords = mesh_nodes[face_nodes[9]] - - # The nodes for an Abaqus element are labeled following a closed path - # around the element: - # - # <---- - # 7 6 5 - # *-----*-----* - # | | - # | | | ^ - # | 8* 9* *4 | - # v | | | - # | | - # *-----*-----* - # 1 2 3 - # ----> - # `curve_values`, however, requires to sort the nodes into a - # valid coordinate system, - # - # *-----*-----* - # | | - # | | - # * * * - # | | - # ^ η | | - # | *-----*-----* - # |----> ξ - # thus we need to flip the node order for the second xi and eta edges met. - - # TODO: Need probably to do some flips for the different faces, as for the 2D case - # Proceed along bottom edge - curve_values[:, 1, 1] = node1_coords - #curve_values[:, 2, 1] = node2_coords - curve_values[:, 2, 1] = 0.5 .* (node1_coords .+ node3_coords) - - curve_values[:, 3, 1] = node3_coords - - # Proceed along middle line - curve_values[:, 1, 2] = node8_coords - #curve_values[:, 2, 2] = node9_coords - curve_values[:, 2, 2] = 0.5 .* (node8_coords .+ node4_coords) - curve_values[:, 3, 2] = node4_coords - - # Proceed along top edge - curve_values[:, 1, 3] = node7_coords - #curve_values[:, 2, 3] = node6_coords - curve_values[:, 2, 3] = 0.5 .* (node7_coords .+ node5_coords) - curve_values[:, 3, 3] = node5_coords - - # Construct the curve interpolant for the current side - face_curves[face] = CurvedFaceT(nodes, bary_weights, copy(curve_values)) - end - - # Create the node coordinates on this particular element - calc_node_coordinates!(node_coordinates, tree, nodes, face_curves) - end - end - end -end - # Given the eight `hex_vertices` for a hexahedral element extract # the four `face_vertices` for a particular `face_index`. function get_vertices_for_bilinear_interpolant!(face_vertices, face_index, hex_vertices) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index fa5b1d8f81a..76419a70899 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -747,7 +747,6 @@ function T8codeMesh(meshfile::AbaqusFile{NDIMS}; if header == " File created by HOHQMesh" # Mesh curvature and boundary naming is handled with additional information available in meshfile connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_hohqmesh_abaqus(meshfile.path, - initial_refinement_level, NDIMS, RealT) # Apply user defined mapping. @@ -757,7 +756,6 @@ function T8codeMesh(meshfile::AbaqusFile{NDIMS}; connectivity, tree_node_coordinates, nodes, boundary_names = p4est_connectivity_from_standard_abaqus(meshfile.path, mapping, polydeg, - initial_refinement_level, NDIMS, RealT, boundary_symbols) From 32a78f12e214e479c00217c2eb47bc291c02345d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Dec 2024 15:08:21 +0100 Subject: [PATCH 18/32] example --- .../elixir_navierstokes_SD70003airfoil.jl | 88 +++++++++---------- src/meshes/p4est_mesh.jl | 6 +- 2 files changed, 46 insertions(+), 48 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl index b93a24101f6..f448351e65c 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -4,9 +4,10 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -U_inf = 0.2 -c_inf = 1.0 +gamma = 1.4 +U_inf = 0.2 +aoa = 4 * pi / 180 rho_inf = 1.4 # with gamma = 1.4 => p_inf = 1.0 Re = 10000.0 @@ -14,11 +15,6 @@ airfoil_cord_length = 1.0 t_c = airfoil_cord_length / U_inf -aoa = 4 * pi / 180 -u_x = U_inf * cos(aoa) -u_y = U_inf * sin(aoa) - -gamma = 1.4 prandtl_number() = 0.72 mu() = rho_inf * U_inf * airfoil_cord_length / Re @@ -28,7 +24,7 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), gradient_variables = GradientVariablesPrimitive()) @inline function initial_condition_mach02_flow(x, t, equations) - # set the freestream flow parameters + # set the freestream flow parameters such that c_inf = 1.0 => Mach 0.2 rho_freestream = 1.4 v1 = 0.19951281005196486 # 0.2 * cos(aoa) @@ -39,61 +35,54 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), prim = SVector(rho_freestream, v1, v2, p_freestream) return prim2cons(prim, equations) end - initial_condition = initial_condition_mach02_flow -# Boundary conditions for free-stream testing -boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) - -velocity_bc_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) -heat_bc = Adiabatic((x, t, equations) -> 0.0) -boundary_condition_airfoil = BoundaryConditionNavierStokesWall(velocity_bc_airfoil, heat_bc) - -polydeg = 3 - surf_flux = flux_hllc vol_flux = flux_chandrashekar -solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, +solver = DGSEM(polydeg = 3, surface_flux = surf_flux, volume_integral = VolumeIntegralFluxDifferencing(vol_flux)) ############################################################################### # Get the uncurved mesh from a file (downloads the file if not available locally) -#path = "/storage/home/daniel/PERK4/SD7003/" +# Get quadratic meshfile +mesh_file_name = "SD7003_2D_Quadratic.inp" -path = "/home/daniel/ownCloud - Döhring, Daniel (1MH1D4@rwth-aachen.de)@rwth-aachen.sciebo.de/Job/Doktorand/Content/Meshes/PERK_mesh/SD7003Laminar/" +mesh_file = Trixi.download("https://gist.githubusercontent.com/DanielDoehring/bd2aa20f7e6839848476a0e87ede9f69/raw/1bc8078b4a57634819dc27010f716e68a225c9c6/SD7003_2D_Quadratic.inp", + joinpath(@__DIR__, mesh_file_name)) -#mesh_file = path * "sd7003_laminar_straight_sided_Trixi.inp" -#mesh_file = path * "sd7003_laminar.inp" +# There is also a linear mesh file available at +# https://gist.githubusercontent.com/DanielDoehring/375df933da8a2081f58588529bed21f0/raw/18592aa90f1c86287b4f742fd405baf55c3cf133/SD7003_2D_Linear.inp boundary_symbols = [:Airfoil, :FarField] +mesh = P4estMesh{2}(mesh_file, boundary_symbols = boundary_symbols) -mesh = P4estMesh{2}(mesh_file, boundary_symbols = boundary_symbols, - initial_refinement_level = 0) +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = Dict(:FarField => boundary_condition_free_stream, - :Airfoil => boundary_condition_slip_wall) +velocity_bc_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) +heat_bc = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_airfoil = BoundaryConditionNavierStokesWall(velocity_bc_airfoil, heat_bc) + +boundary_conditions_hyp = Dict(:FarField => boundary_condition_free_stream, + :Airfoil => boundary_condition_slip_wall) -boundary_conditions_parabolic = Dict(:FarField => boundary_condition_free_stream, - :Airfoil => boundary_condition_airfoil) +boundary_conditions_para = Dict(:FarField => boundary_condition_free_stream, + :Airfoil => boundary_condition_airfoil) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) + boundary_conditions = (boundary_conditions_hyp, + boundary_conditions_para)) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 0.005 * t_c) # Try to get into a state where initial pressure wave is gone -tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +tspan = (0.0, 0.01 * t_c) +#tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -# Choose analysis interval such that roughly every dt_c = 0.005 a record is taken -analysis_interval = 25 # PERK4_Multi, PERKSingle - f_aoa() = aoa f_rho_inf() = rho_inf f_U_inf() = U_inf @@ -113,15 +102,21 @@ lift_coefficient = AnalysisSurfaceIntegral((:Airfoil,), LiftCoefficientPressure(f_aoa(), f_rho_inf(), f_U_inf(), f_linf())) +# For long simulation run, use a large interval. +# For measurements once the simulation has settled in, one should use a +# significantly smaller interval, e.g. 50 to record the drag/lift coefficients. +analysis_interval = 10_000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", save_analysis = true, - analysis_errors = Symbol[], + analysis_errors = Symbol[], # Turn off standard errors analysis_integrals = (drag_coefficient, drag_coefficient_shear_force, lift_coefficient)) -stepsize_callback = StepsizeCallback(cfl = 2.0) # PERK_4 Multi E = 5, ..., 14 +stepsize_callback = StepsizeCallback(cfl = 2.0) + +alive_callback = AliveCallback(alive_interval = 50) # For plots etc save_solution = SaveSolutionCallback(interval = 1_000_000, @@ -130,21 +125,22 @@ save_solution = SaveSolutionCallback(interval = 1_000_000, solution_variables = cons2prim, output_directory = "out") -alive_callback = AliveCallback(alive_interval = 10) - save_restart = SaveRestartCallback(interval = Int(10^7), # Only at end save_final_restart = true) -callbacks = CallbackSet(stepsize_callback, # For measurements: Fixed timestep (do not use this) - alive_callback, # Not needed for measurement run - save_solution, # For plotting during measurement run - #save_restart, # For restart with measurements - summary_callback); +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback, + save_solution, + save_restart); ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); +sol = solve(ode, + CarpenterKennedy2N54(williamson_condition = false, + thread = OrdinaryDiffEq.True()); dt = 1.0, save_everystep = false, callback = callbacks) summary_callback() # print the timer summary diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 41a492e1ce2..57968bfa50e 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -598,7 +598,7 @@ function preprocess_standard_abaqus(meshfile, end # p4est can handle only linear elements. This function checks the `meshfile` -# for quadratic elements (highest order supported by Standard Abaqus) and +# for quadratic elements (highest order supported by standard Abaqus) and # replaces them with linear elements. The higher-order (quadratic) boundaries are handled # "internally" by Trixi as for the HOHQMesh-Abaqus case. function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, @@ -730,6 +730,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, linear_hexes = r"^(C3D8).*$" quadratic_hexes = r"^(C3D27).*$" + # Preprocess the meshfile to remove lower-dimensional elements meshfile_preproc, elements_begin_idx, sets_begin_idx = preprocess_standard_abaqus(meshfile, linear_quads, linear_hexes, @@ -737,7 +738,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, quadratic_hexes, n_dimensions) - # Preprocess the meshfile to replace quadratic elements with linear elements + # Copy of mesh for p4est with linear elements only meshfile_p4est_rdy, mesh_polydeg = preprocess_standard_abaqus_for_p4est(meshfile_preproc, linear_quads, linear_hexes, @@ -787,6 +788,7 @@ function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, mesh_nodes = read_nodes_standard_abaqus(meshfile_preproc, n_dimensions, elements_begin_idx, RealT) + # Extract element section from pre-processed Abaqus meshfile element_lines = readlines(open(meshfile_preproc))[elements_begin_idx:(sets_begin_idx - 1)] if n_dimensions == 2 From 3acd22af73dbb4d2500ea6171a25eb80f5d8746e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Dec 2024 15:29:23 +0100 Subject: [PATCH 19/32] todo --- .../p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl index f448351e65c..554fd4fc4b7 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl @@ -119,13 +119,13 @@ stepsize_callback = StepsizeCallback(cfl = 2.0) alive_callback = AliveCallback(alive_interval = 50) # For plots etc -save_solution = SaveSolutionCallback(interval = 1_000_000, +save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim, output_directory = "out") -save_restart = SaveRestartCallback(interval = Int(10^7), # Only at end +save_restart = SaveRestartCallback(interval = analysis_interval, save_final_restart = true) callbacks = CallbackSet(summary_callback, @@ -138,6 +138,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation +# TODO: PERK ? sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); From 0d0c07d14a596835290ad4b1ad642c9b9a2bab0e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 19 Dec 2024 09:40:41 +0100 Subject: [PATCH 20/32] Mention that boundary is in general only higher-order, not curved --- src/meshes/p4est_mesh.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 57968bfa50e..6fd0b6a0b0c 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -142,7 +142,7 @@ end RealT=Float64, initial_refinement_level=0, periodicity=true, unsaved_changes=true, p4est_partition_allow_for_coarsening=true) -Create a structured curved `P4estMesh` of the specified size. +Create a structured curved/higher-order `P4estMesh` of the specified size. There are three ways to map the mesh to the physical domain. 1. Define a `mapping` that maps the hypercube `[-1, 1]^n`. @@ -319,9 +319,9 @@ mesh from an Abaqus mesh file (`.inp`). Each element of the conforming mesh pars from the `meshfile` is created as a [`p4est`](https://github.com/cburstedde/p4est) tree datatype. -To create a curved unstructured mesh `P4estMesh` two strategies are available: +To create a curved/higher-order unstructured mesh `P4estMesh` two strategies are available: -- `p4est_mesh_from_hohqmesh_abaqus`: High-order, curved boundary information created by +- `p4est_mesh_from_hohqmesh_abaqus`: High-order, polygonial boundary information created by [`HOHQMesh.jl`](https://github.com/trixi-framework/HOHQMesh.jl) is available in the `meshfile`. The mesh polynomial degree `polydeg` of the boundaries is provided from the `meshfile`. The computation of @@ -1770,7 +1770,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, nnodes = length(nodes) # Setup the starting file index to read in element indices and the additional - # curved boundary information provided by HOHQMesh. + # higher-order, polygonial boundary information provided by HOHQMesh. file_idx = findfirst(contains("** mesh polynomial degree"), file_lines) + 1 # Create a work set of Gamma curves to create the node coordinates @@ -1803,7 +1803,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, for i in 1:4 quad_vertices[i, :] .= vertices[1:2, element_node_ids[i]] # 2D => 1:2 end - # Pull the information to check if boundary is curved in order to read in additional data + # Pull the information to check if boundary is curved/higher-order in order to read in additional data file_idx += 1 current_line = split(file_lines[file_idx]) # Note: This strategy is HOHQMesh-Abaqus.inp specific! @@ -1815,7 +1815,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, # Create the node coordinates on this particular element calc_node_coordinates!(node_coordinates, tree, nodes, quad_vertices) else - # Quadrilateral element has at least one curved side + # Quadrilateral element has at least one curved/higher-order side # Flip node ordering to make sure the element is right-handed for the interpolations m1 = 1 m2 = 2 @@ -1840,7 +1840,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, 2]) end else - # When curved_check[i] is 1 this curved boundary information is supplied by the mesh + # When curved_check[i] is 1 this curved/higher-order 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 @@ -2003,7 +2003,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, nnodes = length(nodes) # Setup the starting file index to read in element indices and the additional - # curved boundary information provided by HOHQMesh. + # curved/higher-order boundary information provided by HOHQMesh. file_idx = findfirst(contains("** mesh polynomial degree"), file_lines) + 1 # Create a work set of Gamma curves to create the node coordinates @@ -2040,7 +2040,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, for i in 1:8 hex_vertices[:, i] .= vertices[:, element_node_ids[i]] end - # Pull the information to check if boundary is curved in order to read in additional data + # Pull the information to check if boundary is curved/higher-order in order to read in additional data file_idx += 1 current_line = split(file_lines[file_idx]) curved_check[1] = parse(Int, current_line[2]) @@ -2053,7 +2053,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, # Create the node coordinates on this element calc_node_coordinates!(node_coordinates, tree, nodes, hex_vertices) else - # Hexahedral element has at least one curved side + # Hexahedral element has at least one curved/higher-order side for face in 1:6 if curved_check[face] == 0 # Face is a flat plane. @@ -2067,7 +2067,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, nodes[q]) end else # curved_check[face] == 1 - # Curved face boundary information is supplied by + # Curved/higher-order face boundary information is supplied by # the mesh file. Just read it into a work array for q in 1:nnodes, p in 1:nnodes file_idx += 1 From 393424aee63e1e9bd82ab29781179d67bb514818 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 13:33:21 +0100 Subject: [PATCH 21/32] rename --- ...irfoil.jl => elixir_navierstokes_SD7003airfoil.jl} | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) rename examples/p4est_2d_dgsem/{elixir_navierstokes_SD70003airfoil.jl => elixir_navierstokes_SD7003airfoil.jl} (93%) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl similarity index 93% rename from examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl rename to examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl index 554fd4fc4b7..2dd36ffc228 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD70003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl @@ -77,8 +77,10 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 0.01 * t_c) -#tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +# Drag/Lift coefficient measurements should then be done over the 30 to 35 t_c interval +# by restarting the simulation. + ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -104,7 +106,7 @@ lift_coefficient = AnalysisSurfaceIntegral((:Airfoil,), # For long simulation run, use a large interval. # For measurements once the simulation has settled in, one should use a -# significantly smaller interval, e.g. 50 to record the drag/lift coefficients. +# significantly smaller interval, e.g. TODO 50 to record the drag/lift coefficients. analysis_interval = 10_000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", @@ -114,7 +116,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, drag_coefficient_shear_force, lift_coefficient)) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 2.5) alive_callback = AliveCallback(alive_interval = 50) @@ -138,7 +140,6 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -# TODO: PERK ? sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()); From 41e5c3132bbb3ac9b8ffd5c417731ddf359d007b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 13:44:22 +0100 Subject: [PATCH 22/32] test + example --- .../elixir_navierstokes_SD7003airfoil.jl | 10 ++++--- test/test_parabolic_2d.jl | 26 +++++++++++++++++++ 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl index 2dd36ffc228..e8b22d8bfed 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl @@ -77,7 +77,10 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 30 * t_c) # Try to get into a state where initial pressure wave is gone +# Run simulation until initial pressure wave is gone. +# Note: This is a very long simulation! +tspan = (0.0, 30 * t_c) + # Drag/Lift coefficient measurements should then be done over the 30 to 35 t_c interval # by restarting the simulation. @@ -106,7 +109,7 @@ lift_coefficient = AnalysisSurfaceIntegral((:Airfoil,), # For long simulation run, use a large interval. # For measurements once the simulation has settled in, one should use a -# significantly smaller interval, e.g. TODO 50 to record the drag/lift coefficients. +# significantly smaller interval, e.g. 500 to record the drag/lift coefficients. analysis_interval = 10_000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", @@ -116,11 +119,10 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, drag_coefficient_shear_force, lift_coefficient)) -stepsize_callback = StepsizeCallback(cfl = 2.5) +stepsize_callback = StepsizeCallback(cfl = 2.2) alive_callback = AliveCallback(alive_interval = 50) -# For plots etc save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 62142e1fc3b..5d2012d5e37 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -785,6 +785,32 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_navierstokes_SD7003airfoil.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_SD7003airfoil.jl"), + l2=[ + 9.292899618740586e-5, + 0.0001350510200255721, + 7.964907891113045e-5, + 0.0002336568736996096 + ], + linf=[ + 0.2845637352223691, + 0.295808392241858, + 0.19309201225626166, + 0.7188927326929244 + ], + tspan=(0.0, 5e-3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory From fc2d1f1c0bd9ae22545def77ab862f5a48606d1d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 13:53:36 +0100 Subject: [PATCH 23/32] docs --- docs/literate/src/files/p4est_from_gmsh.jl | 2 +- docs/src/meshes/p4est_mesh.md | 1 + src/meshes/p4est_mesh.jl | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/p4est_from_gmsh.jl b/docs/literate/src/files/p4est_from_gmsh.jl index caaea4faa9a..7995902404c 100644 --- a/docs/literate/src/files/p4est_from_gmsh.jl +++ b/docs/literate/src/files/p4est_from_gmsh.jl @@ -315,7 +315,7 @@ end #hide #md # ``` # which forces `gmsh` to generate quadrilateral elements instead of the default triangles. # This is strictly required to be able to use the mesh later with `p4est`, which supports only straight-sided quads, -# i.e., `C2D4, CPS4, S4` in 2D and `C3D` in 3D. +# i.e., `CPS4, S4, ...` in 2D and `C3D8` in 3D. # See for more details the (short) [documentation](https://p4est.github.io/p4est-howto.pdf) on the interaction of `p4est` with `.inp` files. # In principle, it should also be possible to use the `recombine` function of `gmsh` to convert the triangles to quads, # but this is observed to be less robust than enforcing quads from the beginning. diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index a14551b3f46..915d9547573 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -386,6 +386,7 @@ transfinite map of the straight sided hexahedral element to find Also for a mesh in standard Abaqus format there are no qualitative changes when going from 2D to 3D. The most notable difference is that boundaries are formed in 3D by faces defined by four nodes while in 2D boundaries are edges consisting of two elements. +Note that standard Abaqus also defines quadratic element types. In Trixi.jl, these higher-order elements are currently only supported in 2D, i.e., 8-node quadrilaterals. A simple mesh file, which is used also in `examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl`, is given below: ``` *Heading diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 6fd0b6a0b0c..3e2e0877c48 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -659,7 +659,7 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, new_line = join(parts[1:5], ',') elseif occursin(quadratic_hexes, current_element_type) # Print the first (element), second to ninth (vertices 1-8) indices to file. - # The node order is fortunately the same for hexes/bricks of type "C3D20", "C3D27", see + # The node order is fortunately the same for hexes/bricks of type "C3D27", see # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael03.html new_line = join(parts[1:9], ',') end From dcab80cb80e9296c05672e2771c3bd082120eab3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 14:02:40 +0100 Subject: [PATCH 24/32] Comment --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 3e2e0877c48..b8e3d9d387e 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -597,7 +597,7 @@ function preprocess_standard_abaqus(meshfile, return meshfile_preproc, elements_begin_idx, sets_begin_idx end -# p4est can handle only linear elements. This function checks the `meshfile` +# p4est can handle only linear elements. This function checks the `meshfile_pre_proc` # for quadratic elements (highest order supported by standard Abaqus) and # replaces them with linear elements. The higher-order (quadratic) boundaries are handled # "internally" by Trixi as for the HOHQMesh-Abaqus case. From 09e92c6da7ceb5b4ee08ea8d76a48ab58cf988cc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 14:03:19 +0100 Subject: [PATCH 25/32] reset project toml --- Project.toml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Project.toml b/Project.toml index e200188003e..296883f5420 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -25,13 +24,11 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -54,7 +51,6 @@ TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [weakdeps] Convex = "f65535da-76fb-5f13-bab9-19810c17039a" From 8b4173f7aa34a0dc2a8909022f303c1bc9ca3741 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 14:09:48 +0100 Subject: [PATCH 26/32] news --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index c187f229519..6ed5a9716af 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,7 @@ for human readability. and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008]) - `LobattoLegendreBasis` and related datastructures made fully floating-type general, enabling calculations with higher than double (`Float64`) precision ([#2128]) +- In 2D, quadratic elements, i.e., 8-node (quadratic) quadrilaterals and 6-node are now supported in standard Abaqus `inp` format ([#2217]) ## Changes when updating to v0.9 from v0.8.x From b2bd8cb6ff52fcb5baba171905b9e2d4e2c2f244 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 20 Dec 2024 14:11:18 +0100 Subject: [PATCH 27/32] typo --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index b8e3d9d387e..4f5ede41517 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -654,7 +654,7 @@ function preprocess_standard_abaqus_for_p4est(meshfile_pre_proc, parts = split(line, ',') if occursin(quadratic_quads, current_element_type) # Print the first (element), second to fifth (vertices 1-4) indices to file. - # For node order of quadratic (secod-order) quads, see + # For node order of quadratic (second-order) quads, see # http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt06ch28s01ael02.html new_line = join(parts[1:5], ',') elseif occursin(quadratic_hexes, current_element_type) From b00ea5431ce3b8370f6c4572b547b53f8e1f0afb Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Dec 2024 09:01:35 +0100 Subject: [PATCH 28/32] Update src/meshes/p4est_mesh.jl --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 4f5ede41517..240ed8c1c89 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -707,7 +707,7 @@ end # Create the mesh connectivity, mapped node coordinates within each tree, reference nodes in [-1,1] # and a list of boundary names for the `P4estMesh`. For linear meshes, the tree node coordinates are # computed according to the `mapping` passed to this function using polynomial interpolants of degree `polydeg`. -# For quadratic (second-order) meshes, the tree node coordinates are read from the file, similar as for +# For quadratic (second-order) meshes, the tree node coordinates are read from the meshfile, similar as for # `p4est_connectivity_from_hohqmesh_abaqus`. function p4est_connectivity_from_standard_abaqus(meshfile, mapping, polydeg, n_dimensions, From 523de2b21bae724a5f2d322d4cc3ba27260e5526 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Dec 2024 09:01:46 +0100 Subject: [PATCH 29/32] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 6ed5a9716af..c830ec384d3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,7 +13,7 @@ for human readability. and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008]) - `LobattoLegendreBasis` and related datastructures made fully floating-type general, enabling calculations with higher than double (`Float64`) precision ([#2128]) -- In 2D, quadratic elements, i.e., 8-node (quadratic) quadrilaterals and 6-node are now supported in standard Abaqus `inp` format ([#2217]) +- In 2D, quadratic elements, i.e., 8-node (quadratic) quadrilaterals are now supported in standard Abaqus `inp` format ([#2217]) ## Changes when updating to v0.9 from v0.8.x From 1588a43657fc7cd60dedc521bd5c7b1c475b3cc6 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Dec 2024 09:02:41 +0100 Subject: [PATCH 30/32] Update src/meshes/p4est_mesh.jl Co-authored-by: Andrew Winters --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 240ed8c1c89..f7e66a96335 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -556,7 +556,7 @@ function preprocess_standard_abaqus(meshfile, current_element_type = match(r"\*ELEMENT, type=([^,]+)", line).captures[1] # Keep only quads (2D) or hexes (3D), i.e., eliminate all other elements - # like trusses (2D/3D) or quads in 3D. + # like trusses (2D/3D) or quads (3D). if n_dimensions == 2 && (occursin(linear_quads, current_element_type) || occursin(quadratic_quads, current_element_type)) From 13f29942fecd9ea20a05ddad7e62e34add0e5cd4 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Dec 2024 09:09:03 +0100 Subject: [PATCH 31/32] Update examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl Co-authored-by: Andrew Winters --- examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl index e8b22d8bfed..e52d84f0c0e 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl @@ -27,6 +27,7 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), # set the freestream flow parameters such that c_inf = 1.0 => Mach 0.2 rho_freestream = 1.4 + # Values correspond to `aoa = 4 * pi / 180` v1 = 0.19951281005196486 # 0.2 * cos(aoa) v2 = 0.01395129474882506 # 0.2 * sin(aoa) From 35100a4ccc98b219ec2fa8e32a0f7f733917bb64 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Dec 2024 09:09:19 +0100 Subject: [PATCH 32/32] Update src/meshes/p4est_mesh.jl --- src/meshes/p4est_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index f7e66a96335..1fcac40275c 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -454,7 +454,7 @@ function p4est_connectivity_from_hohqmesh_abaqus(meshfile, n_vertices::Int = connectivity_pw.num_vertices[] # Extract a copy of the element vertices to compute the tree node coordinates - # `vertices` store coordinates of all three dimensions (even for the 2D case) + # `vertices` stores coordinates of all three dimensions (even for the 2D case) # since the Abaqus `.inp` format always stores 3D coordinates. vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices))