adding method to write mesh to file.

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))

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")

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