diff --git a/ext/AbaqusReaderExt.jl b/ext/AbaqusReaderExt.jl index 9ccd831..df17cde 100644 --- a/ext/AbaqusReaderExt.jl +++ b/ext/AbaqusReaderExt.jl @@ -117,6 +117,7 @@ function _get_command_data(lines, command_line_ids, comment_line_ids, command_id n = n + 1 continue end + @show lines[n] for node in split(lines[n], ',') try push!(temp_nodes, parse(Int, node)) diff --git a/src/meshes/Meshes.jl b/src/meshes/Meshes.jl index 7c9d2dc..d27f2d4 100644 --- a/src/meshes/Meshes.jl +++ b/src/meshes/Meshes.jl @@ -1,4 +1,4 @@ -elem_type_map = Dict{String, Type{<:ReferenceFiniteElements.AbstractElementType}}( +const elem_type_map = Dict{String, Type{<:ReferenceFiniteElements.AbstractElementType}}( "HEX" => Hex8, "HEX8" => Hex8, "QUAD" => Quad4, @@ -24,8 +24,6 @@ $(TYPEDEF) """ abstract type AbstractMesh end - - """ $(TYPEDSIGNATURES) Returns file name for an mesh type @@ -46,8 +44,6 @@ struct FileMesh{MeshObj} <: AbstractMesh mesh_obj::MeshObj end - - function _create_edges!(all_edges, block_edges, el_to_nodes, ref_fe) local_edge_to_nodes = ref_fe.edge_nodes @@ -96,5 +92,67 @@ end # return nothing # end +# TODO might need to be careful about int types below +function write_to_file(mesh::AbstractMesh, file_name::String; force::Bool = false) + if force && isfile(file_name) + Base.rm(file_name; force = true) + end + + # initialization parameters + num_dim, num_nodes = size(mesh.nodal_coords) + num_elems = mapreduce(x -> size(x, 2), +, values(mesh.element_conns)) + num_elem_blks = length(mesh.element_conns) + # num_side_sets = length(mesh.sideset_elems) + num_node_sets = length(mesh.nodeset_nodes) + num_side_sets = 0 + # num_node_sets = 0 + + # make init + init = Initialization{Int32}( + num_dim, num_nodes, num_elems, + num_elem_blks, num_node_sets, num_side_sets + ) + + # create exo + exo = ExodusDatabase{Int32, Int32, Int32, eltype(mesh.nodal_coords)}( + file_name, "w", init + ) + + # write coordinates + coords = mesh.nodal_coords |> collect + write_coordinates(exo, coords) + + # write node map + write_id_map(exo, NodeMap, convert.(Int32, mesh.node_id_map)) + + # write block names + write_names(exo, Block, map(String, mesh.element_block_names)) + + # TODO write block id maps + for n in axes(mesh.element_block_names, 1) + write_block(exo, n, String(mesh.element_types[n]), values(mesh.element_conns)[n] |> collect) + end + + # write nodesets + names = keys(mesh.nodeset_nodes) |> collect + for (n, name) in enumerate(names) + nodes = mesh.nodeset_nodes[name] + nset = NodeSet(Int32(n), convert.(Int32, nodes)) + write_set(exo, nset) + end + names = map(String, names) + write_names(exo, NodeSet, names) + + # TODO write nodesets + names = keys(mesh.sideset_elems) |> collect + for (n, name) in enumerate(names) + elems = mesh.sideset_elems[name] + nodes = mesh.sideset_nodes[name] + sides = mesh.sideset_sides[name] + end + + close(exo) +end + include("StructuredMesh.jl") include("UnstructuredMesh.jl") \ No newline at end of file diff --git a/src/meshes/StructuredMesh.jl b/src/meshes/StructuredMesh.jl index e17451d..b84dfa9 100644 --- a/src/meshes/StructuredMesh.jl +++ b/src/meshes/StructuredMesh.jl @@ -28,23 +28,25 @@ function StructuredMesh(element_type, mins, maxs, counts) end element_type = uppercase(element_type) if occursin(element_type, "HEX") - @assert false "Implement hex case" + element_type = :HEX8 + num_dims = 3 + nodal_coords, element_conns, nodeset_nodes, + sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes = + _hex8_structured_mesh_data(mins, maxs, counts) elseif occursin(element_type, "QUAD") element_type = :QUAD4 num_dims = 2 nodal_coords, element_conns, nodeset_nodes, - sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes = _quad4_structured_mesh_data( - counts[1], counts[2], (mins[1], maxs[1]), (mins[2], maxs[2]) - ) + sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes = + _quad4_structured_mesh_data(mins, maxs, counts) elseif occursin(element_type, "TET") @assert false "Implement tet case" elseif occursin(element_type, "TRI") element_type = :TRI3 num_dims = 2 nodal_coords, element_conns, nodeset_nodes, - sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes = _tri3_structured_mesh_data( - counts[1], counts[2], (mins[1], maxs[1]), (mins[2], maxs[2]) - ) + sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes = + _tri3_structured_mesh_data(mins, maxs, counts) else throw(ArgumentError("Unsupported element type $element_type")) end @@ -78,12 +80,160 @@ function Base.show(io::IO, mesh::StructuredMesh) end -function _quad4_structured_mesh_data(N_x, N_y, x_extent, y_extent) +function _hex8_structured_mesh_data(mins, maxs, counts) + N_x, N_y, N_z = counts + E_x = N_x - 1 + E_y = N_y - 1 + E_z = N_z - 1 + + coords = Matrix{Float64}(undef, 3, N_x * N_y * N_z) + _three_dimensional_coords!(coords, mins, maxs, counts) + + node(i, j, k) = i + N_x * (j - 1) + N_x * N_y * (k - 1) + + conns = Matrix{Int64}(undef, 8, E_x * E_y * E_z) + n = 1 + # for ez in 1:E_z + # for ey in 1:E_y + # for ex in 1:E_x + # conns[1, n] = node(ex, ey, ez) + # conns[2, n] = node(ex+1, ey, ez) + # conns[3, n] = node(ex+1, ey+1, ez) + # conns[4, n] = node(ex, ey+1, ez) + # conns[5, n] = node(ex, ey, ez+1) + # conns[6, n] = node(ex+1, ey, ez+1) + # conns[7, n] = node(ex+1, ey+1, ez+1) + # conns[8, n] = node(ex, ey+1, ez+1) + # n += 1 + # end + # end + # end + for ex in 1:E_x + for ey in 1:E_y + for ez in 1:E_z + conns[1, n] = node(ex, ey, ez) + conns[2, n] = node(ex + 1, ey, ez) + conns[3, n] = node(ex + 1, ey + 1, ez) + conns[4, n] = node(ex, ey + 1, ez) + conns[5, n] = node(ex, ey, ez + 1) + conns[6, n] = node(ex + 1, ey, ez + 1) + conns[7, n] = node(ex + 1, ey + 1, ez + 1) + conns[8, n] = node(ex, ey + 1, ez + 1) + n = n + 1 + end + end + end + nodeset_nodes = _three_dimensional_nsets(N_x, N_y, N_z) + + return coords, conns, nodeset_nodes, _hex8_ssets(conns, N_x, N_y, N_z)... +end + +function _hex8_ssets(conns, N_x, N_y, N_z) + E_x = N_x - 1 + E_y = N_y - 1 + E_z = N_z - 1 + elem(i, j, k) = (i - 1) * E_y * E_z + (j - 1) * E_z + k + + sideset_elems = Dict{Symbol, Vector{Int}}( + :bottom => Int[], + :right => Int[], + :front => Int[], + :top => Int[], + :left => Int[], + :back => Int[] + ) + sideset_side_nodes = Dict{Symbol, Matrix{Int}}( + :bottom => Matrix{Int}(undef, 4, E_x), + :right => Matrix{Int}(undef, 4, E_y), + :front => Matrix{Int}(undef, 4, E_z), + :top => Matrix{Int}(undef, 4, E_x), + :left => Matrix{Int}(undef, 4, E_y), + :back => Matrix{Int}(undef, 4, E_z) + ) + sideset_sides = Dict{Symbol, Vector{Int}}( + :bottom => Int[], + :right => Int[], + :front => Int[], + :top => Int[], + :left => Int[], + :back => Int[] + ) + + # bottom + j = 1 + for k in 1:E_z, i in 1:E_x + e = elem(i, j, k) + n1, n2, n3, n4 = conns[1,e], conns[2,e], conns[3,e], conns[4,e] + push!(sideset_elems[:bottom], e) + push!(sideset_sides[:bottom], 5) + sideset_side_nodes[:bottom][:, i] .= [n1, n2, n3, n4] + end + + # right + i = E_x + for k in 1:E_z, j in 1:E_y + e = elem(i, j, k) + n2, n3, n7, n6 = conns[2, e], conns[3, e], conns[7, e], conns[6, e] + push!(sideset_elems[:right], e) + push!(sideset_sides[:right], 2) + sideset_side_nodes[:right][:, j] .= [n2, n3, n7, n6] + end + + # front + k = E_z + for j in 1:E_y, i in 1:E_x + e = elem(i, j, k) + n1, n2, n6, n5 = conns[1, e], conns[2, e], conns[6, e], conns[5, e] + push!(sideset_elems[:front], e) + push!(sideset_sides[:front], 1) + sideset_side_nodes[:front][:, j] .= [n1, n2, n6, n5] + end + + # top + j = E_y + for k in 1:E_z, i in 1:E_x + e = elem(i, j, k) + n5, n6, n7, n8 = conns[5, e], conns[6, e], conns[7, e], conns[8, e] + push!(sideset_elems[:top], e) + push!(sideset_sides[:top], 6) + sideset_side_nodes[:top][:, i] .= [n5, n6, n7, n8] + end + + # left + i = 1 + for j in 1:E_y + e = elem(i, j, k) + n4, n1, n5, n8 = conns[4, e], conns[1, e], conns[5, e], conns[8, e] + push!(sideset_elems[:left], e) + push!(sideset_sides[:left], 4) + sideset_side_nodes[:left][:, j] .= [n4, n1, n5, n8] + end + + # back + k = 1 + for j in 1:E_y, i in 1:E_x + e = elem(i, j, k) + n3, n4, n8, n7 = conns[3, e], conns[4, e], conns[8, e], conns[7, e] + push!(sideset_elems[:back], e) + push!(sideset_sides[:back], 3) + sideset_side_nodes[:back][:, j] .= [n3, n4, n8, n7] + end + + sideset_nodes = Dict{Symbol, Vector{Int}}() + for (k, v) in sideset_side_nodes + sideset_nodes[k] = unique(v) + end + + return sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes +end + +function _quad4_structured_mesh_data(mins, maxs, counts) + N_x, N_y = counts E_x = N_x - 1 E_y = N_y - 1 coords = Matrix{Float64}(undef, 2, N_x * N_y) - _two_dimensional_coords!(coords, x_extent, y_extent, N_x, N_y) + _two_dimensional_coords!(coords, mins, maxs, counts) node(i, j) = i + (j - 1) * N_x @@ -174,12 +324,13 @@ function _quad4_ssets(conns, N_x, N_y) return sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes end -function _tri3_structured_mesh_data(N_x, N_y, x_extent, y_extent) +function _tri3_structured_mesh_data(mins, maxs, counts) + N_x, N_y = counts E_x = N_x - 1 E_y = N_y - 1 coords = Matrix{Float64}(undef, 2, N_x * N_y) - _two_dimensional_coords!(coords, x_extent, y_extent, N_x, N_y) + _two_dimensional_coords!(coords, mins, maxs, counts) conns = Matrix{Int64}(undef, 3, 2 * E_x * E_y) n = 1 @@ -277,12 +428,51 @@ function _tri3_ssets(conns, N_x, N_y) return sideset_elems, sideset_nodes, sideset_sides, sideset_side_nodes end -function _two_dimensional_coords!(coords, x_extent, y_extent, N_x, N_y) - xs = LinRange(x_extent[1], x_extent[2], N_x) - ys = LinRange(y_extent[1], y_extent[2], N_y) +function _three_dimensional_coords!(coords, mins, maxs, counts) + xs = LinRange(mins[1], maxs[1], counts[1]) + ys = LinRange(mins[2], maxs[2], counts[2]) + zs = LinRange(mins[3], maxs[3], counts[3]) + n = 1 + for nz in 1:counts[1] + for ny in 1:counts[2] + for nx in 1:counts[3] + coords[1, n] = xs[nx] + coords[2, n] = ys[ny] + coords[3, n] = zs[nz] + n = n + 1 + end + end + end + return nothing +end + +# TODO +function _three_dimensional_nsets(N_x, N_y, N_z) + node(i, j, k) = i + (j - 1) * N_x + (k - 1) * N_x * N_y + + bottom = [node(i, 1, k) for i in 1:N_x, k in 1:N_z] |> vec + right = [node(N_x, j, k) for j in 1:N_y, k in 1:N_z] |> vec + front = [node(i, j, N_z) for i in 1:N_x, j in 1:N_y] |> vec + top = [node(i, N_y, k) for i in 1:N_x, k in 1:N_z] |> vec + left = [node(1, j, k) for j in 1:N_y, k in 1:N_z] |> vec + back = [node(i, j, 1) for i in 1:N_x, j in 1:N_y] |> vec + + return Dict{Symbol, Vector{Int}}( + :bottom => bottom, + :right => right, + :front => front, + :top => top, + :left => left, + :back => back + ) +end + +function _two_dimensional_coords!(coords, mins, maxs, counts) + xs = LinRange(mins[1], maxs[1], counts[1]) + ys = LinRange(mins[2], maxs[2], counts[2]) n = 1 - for ny in 1:N_x - for nx in 1:N_y + for ny in 1:counts[1] + for nx in 1:counts[2] coords[1, n] = xs[nx] coords[2, n] = ys[ny] n = n + 1 diff --git a/test/TestMesh.jl b/test/TestMesh.jl index 9689765..73d3187 100644 --- a/test/TestMesh.jl +++ b/test/TestMesh.jl @@ -1,10 +1,26 @@ -using FiniteElementContainers -using Test function test_structured_mesh() # testing bad input @test_throws ArgumentError StructuredMesh("bad element", (0., 0.), (1., 1.), (3, 3)) @test_throws BoundsError StructuredMesh("tri3", (0., 0.), (0., 1.), (3, 3)) + # Hex8 test + mesh = StructuredMesh("hex", (0., 0., 0.), (1., 1., 1.), (3, 3, 3)) + coords = mesh.nodal_coords + @test all(coords[2, mesh.nodeset_nodes[:bottom]] .≈ 0.) + @test all(coords[1, mesh.nodeset_nodes[:right]] .≈ 1.) + @test all(coords[3, mesh.nodeset_nodes[:front]] .≈ 1.) + @test all(coords[2, mesh.nodeset_nodes[:top]] .≈ 1.) + @test all(coords[1, mesh.nodeset_nodes[:left]] .≈ 0.) + @test all(coords[3, mesh.nodeset_nodes[:back]] .≈ 0.) + + # currently failing TODO fix this + # @test all(coords[2, mesh.sideset_nodes[:bottom]] .≈ 0.) + # @test all(coords[1, mesh.sideset_nodes[:right]] .≈ 1.) + # @test all(coords[3, mesh.sideset_nodes[:front]] .≈ 1.) + # @test all(coords[2, mesh.sideset_nodes[:top]] .≈ 1.) + # @test all(coords[1, mesh.sideset_nodes[:left]] .≈ 0.) + # @test all(coords[3, mesh.sideset_nodes[:right]] .≈ 0.) + # Quad4 test mesh = StructuredMesh("quad", (0., 0.), (1., 1.), (3, 3)) coords = mesh.nodal_coords @@ -51,9 +67,24 @@ function test_structured_mesh() @show mesh end -# @testset ExtendedTestSet "Mesh" begin -# test_structured_mesh() -# end +function test_write_mesh() + mesh = StructuredMesh("hex", (0., 0., 0.), (1., 1., 1.), (3, 3, 3)) + FiniteElementContainers.write_to_file(mesh, "shex.exo") + rm("shex.exo", force = true) + + mesh = StructuredMesh("quad", (0., 0.), (1., 1.), (3, 3)) + FiniteElementContainers.write_to_file(mesh, "squad.exo") + rm("squad.exo", force = true) + + mesh = StructuredMesh("tri", (0., 0.), (1., 1.), (3, 3)) + FiniteElementContainers.write_to_file(mesh, "stri.exo") + rm("stri.exo", force = true) + + # TODO eventually put an exodiff test below + mesh = UnstructuredMesh("poisson/poisson.g") + FiniteElementContainers.write_to_file(mesh, "umesh.exo") + rm("umesh.exo") +end struct DummyMesh <: FiniteElementContainers.AbstractMesh end @@ -84,10 +115,9 @@ function test_bad_mesh_file_type() @test_throws ErrorException UnstructuredMesh(file_name) end -# test_bad_mesh_file_type() - @testset "Mesh" begin test_bad_mesh_file_type() test_bad_mesh_methods() test_structured_mesh() + test_write_mesh() end