From 67d94231898d5649f82f1dc0b69b2780ace9b462 Mon Sep 17 00:00:00 2001 From: bantingl Date: Wed, 29 Nov 2023 20:59:32 -0600 Subject: [PATCH 1/4] Storing of unstructured grid data in VTKHDF file --- Project.toml | 7 +++ ext/WriteVTKHDFExt.jl | 13 +++++ src/WriteVTK.jl | 6 +++ src/gridtypes/vtk_hdf.jl | 102 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 128 insertions(+) create mode 100644 ext/WriteVTKHDFExt.jl create mode 100644 src/gridtypes/vtk_hdf.jl diff --git a/Project.toml b/Project.toml index 55758812..d2aaf0e6 100644 --- a/Project.toml +++ b/Project.toml @@ -7,10 +7,17 @@ version = "1.18.1" Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" TranscodingStreams = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" VTKBase = "4004b06d-e244-455f-a6ce-a5f9919cc534" +# [weakdeps] +# HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" + +# [extensions] +# WriteVTKHDFExt = "HDF5" + [compat] CodecZlib = "0.7" FillArrays = "0.13, 1.0" diff --git a/ext/WriteVTKHDFExt.jl b/ext/WriteVTKHDFExt.jl new file mode 100644 index 00000000..3b280532 --- /dev/null +++ b/ext/WriteVTKHDFExt.jl @@ -0,0 +1,13 @@ +module WriteVTKHDFExt + + using HDF5 + import WriteVTK: vtk_grid, num_points + using WriteVTK: VTKHDF5, VTKHDFUnstructuredGrid, VTKUnstructuredGrid + using WriteVTK: UnstructuredCoord, CellVector + + __init__() = println("Pkg extension loaded!") + + + + +end \ No newline at end of file diff --git a/src/WriteVTK.jl b/src/WriteVTK.jl index a41ca6cf..b519acf1 100644 --- a/src/WriteVTK.jl +++ b/src/WriteVTK.jl @@ -13,6 +13,9 @@ export VTKPointData, VTKCellData, VTKFieldData export PolyData export pvtk_grid +# vtkhdf needed +export VTKHDF5 + import CodecZlib: ZlibCompressorStream import TranscodingStreams @@ -163,6 +166,9 @@ include("gridtypes/unstructured/polydata.jl") # Parallel DataSet Formats include("gridtypes/pvtk_grid.jl") +# vtkhdf formats +include("gridtypes/vtk_hdf.jl") + # Convenient high-level tools. include("utils/array.jl") include("utils/surface.jl") diff --git a/src/gridtypes/vtk_hdf.jl b/src/gridtypes/vtk_hdf.jl new file mode 100644 index 00000000..136a35d7 --- /dev/null +++ b/src/gridtypes/vtk_hdf.jl @@ -0,0 +1,102 @@ +using HDF5 + +abstract type VTKFileType end +struct VTKHDF5 <: VTKFileType end +struct VTKXML <: VTKFileType end + +struct VTKHDFUnstructuredGrid <: UnstructuredVTKDataset end + +struct VTKHDF5File <: VTKFile + h5file # an HDF5.File + version::AbstractString + grid_type::AbstractString + Npts::Int + Ncls::Int +end + +function vtk_grid( + ::VTKHDF5, + filename, + points, + cells; + kws... +) + vtk_grid( + VTKHDFUnstructuredGrid(), + filename, points, cells; kws...) +end + +function vtk_grid( + dtype::VTKHDFUnstructuredGrid, + filename::AbstractString, + points::UnstructuredCoords, + cells::CellVector; + kwargs...) + + Npts = num_points(dtype, points) + Ncls = length(cells) + + # vtk cell information + IntType = connectivity_type(cells) :: Type{<:Integer} + # Create data arrays. + offsets = Array{IntType}(undef, Ncls+1) + offsets[1] = 0 + + types = Array{UInt8}(undef, Ncls) + + Nconn = 0 # length of the connectivity array + if Ncls >= 1 # it IS possible to have no cells + offsets[2] = length(cells[1].connectivity) + end + + for (n, c) in enumerate(cells) + Npts_cell = length(c.connectivity) + Nconn += Npts_cell + types[n] = cell_type(c).vtk_id + if n >= 2 + offsets[n+1] = offsets[n] + Npts_cell + end + end + display(offsets) + + # Create connectivity array. + conn = Array{IntType}(undef, Nconn) + n = 1 + for c in cells, i in c.connectivity + # We transform to zero-based indexing, required by VTK. + conn[n] = i - one(i) + n += 1 + end + + h5file = h5open("$filename.vtkhdf", "w") + + # write attributes + VTKHDF_group = create_group(h5file, "VTKHDF") + HDF5.attributes(VTKHDF_group)["Version"] = [2, 0] + + type = "UnstructuredGrid" + h5_dspace = HDF5.dataspace(type) + h5_dtype = HDF5.datatype(type) + HDF5.h5t_set_cset(h5_dtype, HDF5.H5T_CSET_ASCII) + attr = create_attribute(VTKHDF_group, "Type", h5_dtype, h5_dspace) + write_attribute(attr, h5_dtype, type) + + VTKHDF_group["NumberOfConnectivityIds"] = [Nconn] + VTKHDF_group["NumberOfPoints"] = [Npts] + VTKHDF_group["NumberOfCells"] = [Ncls] + VTKHDF_group["Points"] = points + VTKHDF_group["Types"] = types + VTKHDF_group["Connectivity"] = conn + VTKHDF_group["Offsets"] = offsets + + VTKHDF5File(h5file, "2.0", type, Npts, Ncls) +end + +# check if file is open +Base.isopen(file::VTKHDF5File) = isopen(file.h5file) + +function Base.close(file::VTKHDF5File) + if isopen(file) + close(file.h5file) + end +end From 133b2cda92450f7b17f18d611f0106df808630f4 Mon Sep 17 00:00:00 2001 From: bantingl Date: Thu, 30 Nov 2023 21:00:33 -0600 Subject: [PATCH 2/4] Time varying unstructured cell data with static mesh --- src/WriteVTK.jl | 2 +- src/gridtypes/vtk_hdf.jl | 166 ++++++++++++++++++++++++++++++++++++--- 2 files changed, 158 insertions(+), 10 deletions(-) diff --git a/src/WriteVTK.jl b/src/WriteVTK.jl index b519acf1..70b9b034 100644 --- a/src/WriteVTK.jl +++ b/src/WriteVTK.jl @@ -14,7 +14,7 @@ export PolyData export pvtk_grid # vtkhdf needed -export VTKHDF5 +export VTKHDF5, vtkhdf_open_timeseries, vtkhdf_append_timeseries_dataset import CodecZlib: ZlibCompressorStream import TranscodingStreams diff --git a/src/gridtypes/vtk_hdf.jl b/src/gridtypes/vtk_hdf.jl index 136a35d7..3a5ed08b 100644 --- a/src/gridtypes/vtk_hdf.jl +++ b/src/gridtypes/vtk_hdf.jl @@ -14,13 +14,27 @@ struct VTKHDF5File <: VTKFile Ncls::Int end +# for unstructured +struct VTKHDFTimeSeries{T <: AbstractFieldData} + vtkhdf::VTKHDF5File + name::AbstractString +end + +# check if file is open +Base.isopen(file::VTKHDF5File) = isopen(file.h5file) + +function Base.close(file::VTKHDF5File) + if isopen(file) + close(file.h5file) + end +end + function vtk_grid( ::VTKHDF5, filename, points, cells; - kws... -) + kws...) vtk_grid( VTKHDFUnstructuredGrid(), filename, points, cells; kws...) @@ -36,6 +50,10 @@ function vtk_grid( Npts = num_points(dtype, points) Ncls = length(cells) + # need types, offsets, connectivity for unstructured + # copied from unstructured.jl, but offset has a zero at the beginnining + # for VTKHDF case + # vtk cell information IntType = connectivity_type(cells) :: Type{<:Integer} # Create data arrays. @@ -57,8 +75,7 @@ function vtk_grid( offsets[n+1] = offsets[n] + Npts_cell end end - display(offsets) - + # Create connectivity array. conn = Array{IntType}(undef, Nconn) n = 1 @@ -92,11 +109,142 @@ function vtk_grid( VTKHDF5File(h5file, "2.0", type, Npts, Ncls) end -# check if file is open -Base.isopen(file::VTKHDF5File) = isopen(file.h5file) +""" +Check for Steps group, maybe create it, and add to Steps. + +Return VTKHDFTimeSeries + +""" +function vtkhdf_open_timeseries( + vtkhdf::VTKHDF5File, + name::AbstractString, + data_type::Union{VTKCellData, VTKPointData}, + vec_dim=1 + ) + root = vtkhdf.h5file["VTKHDF"] + + if !haskey(root, "Steps") + steps = create_group(root, "Steps") + # initialize num_steps to zero + num_steps = 0 + attrs(steps)["NSteps"] = num_steps + + create_dataset(steps, "Values", Float64, dataspace((0,), (-1,)), chunk=(100,)) + + # offsets that are just all zeros + create_dataset(steps, "PartOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) + create_dataset(steps, "PointOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) + create_dataset(steps, "CellOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) + create_dataset(steps, "ConnectivityIdOffsets", Int64, dataspace((1,0), (-1,-1)), chunk=(1,100)) + create_dataset(steps, "NumberOfParts", Int64, dataspace((0,), (-1,)), chunk=(100,)) + else + steps = root["Steps"] + end -function Base.close(file::VTKHDF5File) - if isopen(file) - close(file.h5file) + # check for and maybe create CellDataOffsets of PointDataOffsets respectively + if data_type isa VTKCellData + if vec_dim == 1 + data_size = dataspace((0,),(-1,)) + chunk_size = chunk=(vtkhdf.Ncls,) + else + data_size = dataspace((vec_dim, 0), (-1, -1)) + chunk_size = chunk=(vec_dim, vtkhdf.Ncls) + end + + # where data is stored + if !haskey(root, "CellData") + CellData = create_group(root, "CellData") + else + CellData = root["CellData"] + end + create_dataset(CellData, name, Float64, data_size, chunk=chunk_size) + + if !haskey(steps, "CellDataOffsets") + CellDataOffsets = create_group(steps, "CellDataOffsets") + else + CellDataOffsets = steps["CellDataOffsets"] + end + # where offsets are stored + create_dataset(CellDataOffsets, name, Int64, dataspace((0,), (-1,)), chunk=(100,)) + end + + if data_type isa VTKPointData + if !haskey(root, "PointData") + CellData = create_group(root, "PointData") + else + CellData = root["PointData"] + end + # where data is stored + create_dataset(CellData, name, Float64, (0,)) + + if !haskey(steps, "PointDataOffsets") + PointDataOffsets = create_group(steps, "PointDataOffsets") + else + PointDataOffsets = steps["PointDataOffsets"] + end + create_dataset(PointDataOffsets, name, Int64, (0,)) + + end + + VTKHDFTimeSeries{typeof(data_type)}(vtkhdf, name) +end + +# check for vector case of size(N, M), instead of size(N,) +# also handle chunking? +function append_and_resize(dset, array) + dset_size, dset_max_size = HDF5.get_extent_dims(dset) + + array_dims = size(array) + data_size = array_dims[end] + prev_size = dset_size[end] + + # handle vector case, e.g. (3, N) + if length(array_dims) == 2 + new_size = (array_dims[1], data_size + prev_size) + HDF5.set_extent_dims(dset, new_size) + dset[:,1+prev_size:end] = array + else + new_size = (data_size + prev_size,) + HDF5.set_extent_dims(dset, new_size) + dset[1+prev_size:end] = array end end + +function vtkhdf_append_timeseries_dataset( + series::VTKHDFTimeSeries{VTKCellData}, + time::Float64, + data +) + root = series.vtkhdf.h5file["VTKHDF"] + steps = root["Steps"] + + # increment NSteps + attrs(steps)["NSteps"] += 1 + + # append and resize the underlying datasets + + # Values, append time + append_and_resize(steps["Values"], [time]) + + # CellDataOffsets, append previous offset + length(CellData) + current_len = size(root["CellData"][series.name])[end] + append_and_resize(steps["CellDataOffsets"][series.name], [current_len]) + + # CellData, append 'data' + append_and_resize(root["CellData"][series.name], data) + + # PartOffsets, append 0 + append_and_resize(steps["PartOffsets"], [0]) + + # PointOffsets, append 0 + append_and_resize(steps["PointOffsets"], [0]) + + # CellOffsets, append 0 + append_and_resize(steps["CellOffsets"], [0]) + + # ConnectivityIdOffsets, append 0 + append_and_resize(steps["ConnectivityIdOffsets"], [0][:,:]) + + # NumberOfParts, append 1 + append_and_resize(steps["NumberOfParts"], [1]) +end \ No newline at end of file From 5fe150dc72feffb3ccd4fda2e0848d4bdb73497a Mon Sep 17 00:00:00 2001 From: bantingl Date: Fri, 1 Dec 2023 08:15:49 -0600 Subject: [PATCH 3/4] Selectively increment Steps datasets to allow for multiple datasets for each time step --- src/gridtypes/vtk_hdf.jl | 94 ++++++++++++++++++++++------------------ 1 file changed, 53 insertions(+), 41 deletions(-) diff --git a/src/gridtypes/vtk_hdf.jl b/src/gridtypes/vtk_hdf.jl index 3a5ed08b..7b87f8ae 100644 --- a/src/gridtypes/vtk_hdf.jl +++ b/src/gridtypes/vtk_hdf.jl @@ -1,12 +1,12 @@ using HDF5 -abstract type VTKFileType end +abstract type VTKFileType end struct VTKHDF5 <: VTKFileType end struct VTKXML <: VTKFileType end struct VTKHDFUnstructuredGrid <: UnstructuredVTKDataset end -struct VTKHDF5File <: VTKFile +struct VTKHDF5File <: VTKFile h5file # an HDF5.File version::AbstractString grid_type::AbstractString @@ -15,7 +15,7 @@ struct VTKHDF5File <: VTKFile end # for unstructured -struct VTKHDFTimeSeries{T <: AbstractFieldData} +struct VTKHDFTimeSeries{T<:AbstractFieldData} vtkhdf::VTKHDF5File name::AbstractString end @@ -32,7 +32,7 @@ end function vtk_grid( ::VTKHDF5, filename, - points, + points, cells; kws...) vtk_grid( @@ -41,9 +41,9 @@ function vtk_grid( end function vtk_grid( - dtype::VTKHDFUnstructuredGrid, + dtype::VTKHDFUnstructuredGrid, filename::AbstractString, - points::UnstructuredCoords, + points::UnstructuredCoords, cells::CellVector; kwargs...) @@ -55,9 +55,9 @@ function vtk_grid( # for VTKHDF case # vtk cell information - IntType = connectivity_type(cells) :: Type{<:Integer} + IntType = connectivity_type(cells)::Type{<:Integer} # Create data arrays. - offsets = Array{IntType}(undef, Ncls+1) + offsets = Array{IntType}(undef, Ncls + 1) offsets[1] = 0 types = Array{UInt8}(undef, Ncls) @@ -66,7 +66,7 @@ function vtk_grid( if Ncls >= 1 # it IS possible to have no cells offsets[2] = length(cells[1].connectivity) end - + for (n, c) in enumerate(cells) Npts_cell = length(c.connectivity) Nconn += Npts_cell @@ -75,7 +75,7 @@ function vtk_grid( offsets[n+1] = offsets[n] + Npts_cell end end - + # Create connectivity array. conn = Array{IntType}(undef, Nconn) n = 1 @@ -91,14 +91,14 @@ function vtk_grid( VTKHDF_group = create_group(h5file, "VTKHDF") HDF5.attributes(VTKHDF_group)["Version"] = [2, 0] - type = "UnstructuredGrid" + type = "UnstructuredGrid" h5_dspace = HDF5.dataspace(type) h5_dtype = HDF5.datatype(type) HDF5.h5t_set_cset(h5_dtype, HDF5.H5T_CSET_ASCII) attr = create_attribute(VTKHDF_group, "Type", h5_dtype, h5_dspace) write_attribute(attr, h5_dtype, type) - VTKHDF_group["NumberOfConnectivityIds"] = [Nconn] + VTKHDF_group["NumberOfConnectivityIds"] = [Nconn] VTKHDF_group["NumberOfPoints"] = [Npts] VTKHDF_group["NumberOfCells"] = [Ncls] VTKHDF_group["Points"] = points @@ -118,9 +118,9 @@ Return VTKHDFTimeSeries function vtkhdf_open_timeseries( vtkhdf::VTKHDF5File, name::AbstractString, - data_type::Union{VTKCellData, VTKPointData}, + data_type::Union{VTKCellData,VTKPointData}, vec_dim=1 - ) +) root = vtkhdf.h5file["VTKHDF"] if !haskey(root, "Steps") @@ -130,12 +130,12 @@ function vtkhdf_open_timeseries( attrs(steps)["NSteps"] = num_steps create_dataset(steps, "Values", Float64, dataspace((0,), (-1,)), chunk=(100,)) - + # offsets that are just all zeros - create_dataset(steps, "PartOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) + create_dataset(steps, "PartOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) create_dataset(steps, "PointOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) create_dataset(steps, "CellOffsets", Int64, dataspace((0,), (-1,)), chunk=(100,)) - create_dataset(steps, "ConnectivityIdOffsets", Int64, dataspace((1,0), (-1,-1)), chunk=(1,100)) + create_dataset(steps, "ConnectivityIdOffsets", Int64, dataspace((1, 0), (-1, -1)), chunk=(1, 100)) create_dataset(steps, "NumberOfParts", Int64, dataspace((0,), (-1,)), chunk=(100,)) else steps = root["Steps"] @@ -144,11 +144,11 @@ function vtkhdf_open_timeseries( # check for and maybe create CellDataOffsets of PointDataOffsets respectively if data_type isa VTKCellData if vec_dim == 1 - data_size = dataspace((0,),(-1,)) - chunk_size = chunk=(vtkhdf.Ncls,) + data_size = dataspace((0,), (-1,)) + chunk_size = (vtkhdf.Ncls,) else data_size = dataspace((vec_dim, 0), (-1, -1)) - chunk_size = chunk=(vec_dim, vtkhdf.Ncls) + chunk_size = (vec_dim, vtkhdf.Ncls) end # where data is stored @@ -158,7 +158,7 @@ function vtkhdf_open_timeseries( CellData = root["CellData"] end create_dataset(CellData, name, Float64, data_size, chunk=chunk_size) - + if !haskey(steps, "CellDataOffsets") CellDataOffsets = create_group(steps, "CellDataOffsets") else @@ -167,7 +167,7 @@ function vtkhdf_open_timeseries( # where offsets are stored create_dataset(CellDataOffsets, name, Int64, dataspace((0,), (-1,)), chunk=(100,)) end - + if data_type isa VTKPointData if !haskey(root, "PointData") CellData = create_group(root, "PointData") @@ -183,7 +183,7 @@ function vtkhdf_open_timeseries( PointDataOffsets = steps["PointDataOffsets"] end create_dataset(PointDataOffsets, name, Int64, (0,)) - + end VTKHDFTimeSeries{typeof(data_type)}(vtkhdf, name) @@ -202,7 +202,7 @@ function append_and_resize(dset, array) if length(array_dims) == 2 new_size = (array_dims[1], data_size + prev_size) HDF5.set_extent_dims(dset, new_size) - dset[:,1+prev_size:end] = array + dset[:, 1+prev_size:end] = array else new_size = (data_size + prev_size,) HDF5.set_extent_dims(dset, new_size) @@ -210,6 +210,12 @@ function append_and_resize(dset, array) end end +Base.setindex!( + series::VTKHDFTimeSeries{VTKCellData}, + data, + time::Float64 +) = vtkhdf_append_timeseries_dataset(series, time, data) + function vtkhdf_append_timeseries_dataset( series::VTKHDFTimeSeries{VTKCellData}, time::Float64, @@ -218,33 +224,39 @@ function vtkhdf_append_timeseries_dataset( root = series.vtkhdf.h5file["VTKHDF"] steps = root["Steps"] - # increment NSteps - attrs(steps)["NSteps"] += 1 - # append and resize the underlying datasets - - # Values, append time - append_and_resize(steps["Values"], [time]) # CellDataOffsets, append previous offset + length(CellData) current_len = size(root["CellData"][series.name])[end] append_and_resize(steps["CellDataOffsets"][series.name], [current_len]) - # CellData, append 'data' append_and_resize(root["CellData"][series.name], data) - # PartOffsets, append 0 - append_and_resize(steps["PartOffsets"], [0]) + # check the number of CellData datasets matches number of time steps + num_datasets = current_len รท size(data)[end] + + if num_datasets == length(steps["Values"]) + # increment NSteps + attrs(steps)["NSteps"] += 1 + + # Values, append time + append_and_resize(steps["Values"], [time]) + + + + # PartOffsets, append 0 + append_and_resize(steps["PartOffsets"], [0]) - # PointOffsets, append 0 - append_and_resize(steps["PointOffsets"], [0]) + # PointOffsets, append 0 + append_and_resize(steps["PointOffsets"], [0]) - # CellOffsets, append 0 - append_and_resize(steps["CellOffsets"], [0]) + # CellOffsets, append 0 + append_and_resize(steps["CellOffsets"], [0]) - # ConnectivityIdOffsets, append 0 - append_and_resize(steps["ConnectivityIdOffsets"], [0][:,:]) + # ConnectivityIdOffsets, append 0 + append_and_resize(steps["ConnectivityIdOffsets"], [0][:, :]) - # NumberOfParts, append 1 - append_and_resize(steps["NumberOfParts"], [1]) + # NumberOfParts, append 1 + append_and_resize(steps["NumberOfParts"], [1]) + end # need some sort of reasonable else condition to see if time steps are in sync or not end \ No newline at end of file From f07c268d983948a0fdf2c892100f3363444c9b54 Mon Sep 17 00:00:00 2001 From: Juan Ignacio Polanco Date: Fri, 1 Dec 2023 16:37:11 +0100 Subject: [PATCH 4/4] Fix tests --- src/gridtypes/vtk_hdf.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gridtypes/vtk_hdf.jl b/src/gridtypes/vtk_hdf.jl index 7b87f8ae..73ae2682 100644 --- a/src/gridtypes/vtk_hdf.jl +++ b/src/gridtypes/vtk_hdf.jl @@ -1,4 +1,4 @@ -using HDF5 +using HDF5: HDF5 abstract type VTKFileType end struct VTKHDF5 <: VTKFileType end