Skip to content

Commit

Permalink
Skip Array data when constructing VTKData (#65)
Browse files Browse the repository at this point in the history
* Skip Array data when constructing VTKData

* Formatting

* Add contribution info
  • Loading branch information
williamjsdavis authored Sep 6, 2024
1 parent c914447 commit 0aa9453
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 13 deletions.
2 changes: 2 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,8 @@ Further contributions to ReadVTK have been made by the following people:
(Johannes-Gutenberg University Mainz, Germany)
* [Matthew Whisenant](https://volweb2.utk.edu/~mwhisena/)
(University of Tennessee, Knoxville)
* [William Davis](https://posgeo.wordpress.com/)
(Terra AI, USA)

## License and contributing
ReadVTK is licensed under the MIT license (see [LICENSE.md](LICENSE.md)).
Expand Down
9 changes: 8 additions & 1 deletion src/ReadVTK.jl
Original file line number Diff line number Diff line change
Expand Up @@ -415,14 +415,15 @@ end


# Retrieve a data section (should be `CellData` or `PointData`) from the VTK file

function get_data_section(vtk_file::VTKFile, section)
names = String[]
data_arrays = XMLElement[]

# Iterate over XML elemens in the section
for xml_element in child_elements(piece(vtk_file)[section][1])
# We do not know how to handle anything other than `DataArray`s
@assert LightXML.name(xml_element) == "DataArray"
LightXML.name(xml_element) != "DataArray" && continue

# Store the name and the XML element for each found data array
push!(names, attribute(xml_element, "Name", required = true))
Expand All @@ -437,6 +438,7 @@ end
get_cell_data(vtk_file::VTKFile)
Retrieve a lightweight `VTKData` object with the cell data of the given VTK file.
Only numeric data (i.e., `DataArray`) elements will be read.
See also: [`VTKData`](@ref), [`get_point_data`](@ref)
"""
Expand All @@ -446,6 +448,7 @@ get_cell_data(vtk_file::VTKFile) = get_data_section(vtk_file, "CellData")
get_cell_data(pvtk_file::PVTKFile)
Retrieve a lightweight vector with `PVTKData` objects with the cell data of the given PVTK files.
Only numeric data (i.e., `DataArray`) elements will be read.
See also: [`PVTKData`](@ref), [`get_cell_data`](@ref)
"""
Expand All @@ -464,6 +467,7 @@ end
get_point_data(vtk_file::VTKFile)
Retrieve a lightweight `VTKData` object with the point data of the given VTK file.
Only numeric data (i.e., `DataArray`) elements will be read.
See also: [`VTKData`](@ref), [`get_cell_data`](@ref)
"""
Expand All @@ -473,6 +477,7 @@ get_point_data(vtk_file::VTKFile) = get_data_section(vtk_file, "PointData")
get_point_data(pvtk_file::PVTKFile)
Retrieve a lightweight vector with `PVTKData` objects with the point data of the given PVTK files.
Only numeric data (i.e., `DataArray`) elements will be read.
See also: [`PVTKData`](@ref), [`get_cell_data`](@ref)
"""
Expand All @@ -490,6 +495,7 @@ end
get_coordinate_data(vtk_file::VTKFile)
Retrieve a lightweight `VTKData` object with the coordinate data of the given VTK file.
Only coordinates of numeric data (i.e., `DataArray`) elements will be read.
See also: [`VTKData`](@ref), [`get_point_data`](@ref), [`get_cell_data`](@ref)
"""
Expand All @@ -499,6 +505,7 @@ get_coordinate_data(vtk_file::VTKFile) = get_data_section(vtk_file, "Coordinates
get_coordinate_data(pvtk_file::PVTKFile)
Retrieve a lightweight `{VTKData` object with the coordinate data of the given VTK file.
Only coordinates of numeric data (i.e., `DataArray`) elements will be read.
See also: [`PVTKData`](@ref), [`get_point_data`](@ref), [`get_cell_data`](@ref)
"""
Expand Down
10 changes: 10 additions & 0 deletions test/pvtk_files.jl
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ for part in 1:4
pvtk["Velocity"] = [x + 2y + 3z + v for v in v_global, x in xs, y in ys, z in zs]
pvtk["Pressure"] = [x + 2y + 3z for x in xs_c, y in ys_c, z in zs_c]
pvtk["Phase"] = [trunc(Int64, x * 15) for x in xs_c, y in ys_c, z in zs_c]
pvtk["Temperature string"] = [x + 2y + 3z > 0.5 ? "high" : "low"
for x in xs, y in ys, z in zs]
return nothing
end
end
Expand Down Expand Up @@ -218,6 +220,14 @@ end
Phase_read = get_data_reshaped(cell_data["Phase"], cell_data = true)
@test P_global == P_read
@test Phase_global == Phase_read

@testset "error on get string cell data" begin
@test_throws KeyError cell_data["Temperature string"]
end

@testset "error on get string point data" begin
@test_throws KeyError point_data["Temperature string"]
end
end

# c) PVTU files
Expand Down
41 changes: 29 additions & 12 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -200,16 +200,25 @@ clean_directory(TEST_EXAMPLES_DIR) = @test_nowarn rm(TEST_EXAMPLES_DIR, recursiv

# generate random data
point_scalar_field = rand(Nx, Ny, Nz)
point_data_name = "Point scalar data"
point_scalar_name = "Point scalar data"

cell_scalar_field = rand(Nx - 1, Ny - 1, Nz - 1)
cell_data_name = "Cell scalar data"
cell_scalar_name = "Cell scalar data"

# generate random string data
point_string_field = map(x -> x > 0.5 ? "high" : "low", point_scalar_field)
point_string_name = "Point string data"

cell_string_field = map(x -> x > 0.5 ? "high" : "low", cell_scalar_field)
cell_string_name = "Cell string data"

# write vti file using WriteVTK
path = joinpath(TEST_EXAMPLES_DIR, "grid")
vtk_grid(path, x, y, z) do vtk
vtk[point_data_name, VTKPointData()] = point_scalar_field # scalar field attached to points
vtk[cell_data_name, VTKCellData()] = cell_scalar_field # scalar field attached to cells
vtk[point_scalar_name, VTKPointData()] = point_scalar_field # scalar field attached to points
vtk[cell_scalar_name, VTKCellData()] = cell_scalar_field # scalar field attached to cells
vtk[point_string_name, VTKPointData()] = point_string_field # string field attached to points
vtk[cell_string_name, VTKCellData()] = cell_string_field # string field attached to cells
return nothing
end

Expand All @@ -230,33 +239,41 @@ clean_directory(TEST_EXAMPLES_DIR) = @test_nowarn rm(TEST_EXAMPLES_DIR, recursiv
end

@testset "get scalar cell data" begin
cell_data_raw = get_data(get_cell_data(vtk)[cell_data_name])
cell_data_raw = get_data(get_cell_data(vtk)[cell_scalar_name])
cell_data_reshaped = reshape(cell_data_raw, ((Nx - 1), (Ny - 1), (Nz - 1)))
cell_data_reshaped1 = get_data_reshaped(get_cell_data(vtk)[cell_data_name],
cell_data_reshaped1 = get_data_reshaped(get_cell_data(vtk)[cell_scalar_name],
cell_data = true)

@test cell_data_reshaped == cell_scalar_field
@test cell_data_reshaped1 == cell_scalar_field
end

@testset "get scalar point data" begin
point_data_raw = get_data(get_point_data(vtk)[point_data_name])
point_data_raw = get_data(get_point_data(vtk)[point_scalar_name])
point_data_reshaped = reshape(point_data_raw, (Nx, Ny, Nz))
point_data_reshaped1 = get_data_reshaped(get_point_data(vtk)[point_data_name])
point_data_reshaped1 = get_data_reshaped(get_point_data(vtk)[point_scalar_name])

@test point_data_reshaped == point_scalar_field
@test point_data_reshaped1 == point_scalar_field
end

@testset "error on get string cell data" begin
@test_throws KeyError get_cell_data(vtk)[cell_string_name]
end

@testset "error on get string point data" begin
@test_throws KeyError get_point_data(vtk)[point_string_name]
end

# generate random 2D data
point_scalar_field = rand(Nx, Ny)
cell_scalar_field = rand(Nx - 1, Ny - 1)

# write 2D vti file using WriteVTK
path = joinpath(TEST_EXAMPLES_DIR, "grid_2D")
vtk_grid(path, x, y) do vtk
vtk[point_data_name, VTKPointData()] = point_scalar_field # scalar field attached to points
vtk[cell_data_name, VTKCellData()] = cell_scalar_field # scalar field attached to cells
vtk[point_scalar_name, VTKPointData()] = point_scalar_field # scalar field attached to points
vtk[cell_scalar_name, VTKCellData()] = cell_scalar_field # scalar field attached to cells
return nothing
end

Expand All @@ -275,13 +292,13 @@ clean_directory(TEST_EXAMPLES_DIR) = @test_nowarn rm(TEST_EXAMPLES_DIR, recursiv
end

@testset "get 2D scalar cell data" begin
cell_data_raw = get_data(get_cell_data(vtk)[cell_data_name])
cell_data_raw = get_data(get_cell_data(vtk)[cell_scalar_name])
cell_data_reshaped = reshape(cell_data_raw, ((Nx - 1), (Ny - 1)))
@test cell_data_reshaped == cell_scalar_field
end

@testset "get 2D scalar point data" begin
point_data_raw = get_data(get_point_data(vtk)[point_data_name])
point_data_raw = get_data(get_point_data(vtk)[point_scalar_name])
point_data_reshaped = reshape(point_data_raw, (Nx, Ny))
@test point_data_reshaped == point_scalar_field
end
Expand Down

0 comments on commit 0aa9453

Please sign in to comment.