# Hexahedral sheet extraction

I’m new to Julia and need to implement the HEXAHEDRAL SHEET EXTRACTION algorithm. Here is the algorithm:

Step 1: Define the twist plane, p, and initial hex element, h0, of the sheet;
Step 2: For each chord segment, si, of h0 on p;
Find neighboring hex element, hi;
If hi is already in the sheet,
continue;
else,
Add hi to the sheet and perform Step 2 on hi;
Step 3: Return hex sheet.

# Identify twist plane
function identify_twist_plane(mesh::HexMesh.Mesh, edge_indices::Tuple{Int, Int})
node1, node2 = edge_indices
shared_hexes = intersect(mesh.nodes[node1].hexahedrons, mesh.nodes[node2].hexahedrons)
return Set(shared_hexes)
end

# Extract sheets
function extract_sheets(mesh::HexMesh.Mesh, twist_plane::Set{Int})
new_hexahedrons = [i in twist_plane ? nothing : hex for (i, hex) in enumerate(mesh.hexahedrons)]
new_hexahedrons = filter(x -> x !== nothing, new_hexahedrons)
return HexMesh.Mesh(mesh.nodes, new_hexahedrons)
end

# Merge nodes
function merge_nodes(mesh::HexMesh.Mesh, node_pairs::Vector{Tuple{Int, Int}})
for (n1, n2) in node_pairs
for hex in mesh.hexahedrons
hex.nodes = [n == n2 ? n1 : n for n in hex.nodes]
end
mesh.nodes[n2] = nothing
end
mesh.nodes = filter(x -> x !== nothing, mesh.nodes)
end

# Perform sheet extraction
function perform_sheet_extraction(mesh::HexMesh.Mesh)
twist_plane = identify_twist_plane(mesh, (1, 2))
extract_sheets(mesh, twist_plane)
end

# Expand twist plane
function expand_twist_plane(mesh::HexMesh.Mesh, twist_plane::Set{Int}, current_hex::Int)
for edge in mesh.hexahedrons[current_hex].nodes
for neighbor_hex in mesh.nodes[edge].hexahedrons
if neighbor_hex in twist_plane
continue
else
push!(twist_plane, neighbor_hex)
expand_twist_plane(mesh, twist_plane, neighbor_hex)
end
end
end
end

# Find hex sheet
function find_hex_sheet(mesh::HexMesh.Mesh, initial_hex::Int, edge_indices::Tuple{Int, Int})
twist_plane = identify_twist_plane(mesh, edge_indices)
expand_twist_plane(mesh, twist_plane, initial_hex)
return twist_plane
end

However, the algorithm is not working properly, and I do not know what to do next. I have made an example call by giving a 3D L-shape as a mesh, and the output was a cube. And this is incorrect because the surface must be preserved.
I hope to get a solution or idea here to solve the problem. Thank you in advance.