From 9727bc14f2f7ec4054398c0069b04c02c61fd9cb Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 23 Apr 2024 14:12:13 +0200 Subject: [PATCH 1/2] add feature to read VTK files --- .gitignore | 1 + Project.toml | 2 ++ src/KernelInterpolation.jl | 3 ++- src/io.jl | 22 ++++++++++++++++++++++ test/test_unit.jl | 15 +++++++++++++-- 5 files changed, 40 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 2a2f96d5..688e35a1 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /docs/build/ **/Manifest.toml +**/*.vtu diff --git a/Project.toml b/Project.toml index b26afcb9..a033a589 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "1.0.0-DEV" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -15,6 +16,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] ForwardDiff = "0.10.36" LinearAlgebra = "1" +ReadVTK = "0.2" RecipesBase = "1.1" SpecialFunctions = "2" StaticArrays = "1" diff --git a/src/KernelInterpolation.jl b/src/KernelInterpolation.jl index 29f60aec..64fd38b6 100644 --- a/src/KernelInterpolation.jl +++ b/src/KernelInterpolation.jl @@ -2,6 +2,7 @@ module KernelInterpolation using ForwardDiff: ForwardDiff using LinearAlgebra: norm, Symmetric, tr +using ReadVTK: VTKFile, get_points, get_point_data, get_data using RecipesBase: RecipesBase, @recipe, @series using SpecialFunctions: besselk, loggamma using StaticArrays: StaticArrays, MVector @@ -32,7 +33,7 @@ export NodeSet, separation_distance, dim, eachdim, values_along_dim, random_hype export interpolation_kernel, nodeset, coefficients, kernel_coefficients, polynomial_coefficients, polynomial_basis, polyvars, system_matrix, interpolate, solve, kernel_inner_product, kernel_norm -export vtk_save +export vtk_save, vtk_read export examples_dir, get_examples, default_example, include_example end diff --git a/src/io.jl b/src/io.jl index 5b979cab..05f3e3fc 100644 --- a/src/io.jl +++ b/src/io.jl @@ -17,3 +17,25 @@ function vtk_save(filename, nodeset::NodeSet, functions...; end end end + +""" + vtk_read(filename) + +Read a set of nodes from a VTK file and return them as a `NodeSet`. Note that the data +will always be returned as a 3D `NodeSet`, even if the data is 1D or 2D. The point data +will be returned as a dictionary with the keys being the names of the data arrays in the VTK file. +""" +function vtk_read(filename) + vtk = VTKFile(filename) + nodeset = NodeSet(transpose(get_points(vtk))) + point_data = Dict{String, AbstractArray}() + # If there are no point_data `get_point_data` will throw an error + # I don't see a way to check for point_data without catching the error + try + for (key, data) in get_point_data(vtk) + point_data[key] = get_data(data) + end + catch + end + return nodeset, point_data +end diff --git a/test/test_unit.jl b/test/test_unit.jl index ee253ed3..326928f3 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -173,7 +173,12 @@ using Plots # Saving the nodeset to a VTK file @test_nowarn vtk_save("nodeset1", nodeset1) - # TODO: Check the VTK file after implementing a vtk_read function + nodeset1_2, point_data = @test_nowarn vtk_read("nodeset1.vtu") + @test length(nodeset1_2) == length(nodeset1) + for i in eachindex(nodeset1) + @test [nodeset1[i]; 0.0] == nodeset1_2[i] + end + @test length(point_data) == 0 @test_nowarn rm("nodeset1.vtu", force = true) nodeset2 = @test_nowarn NodeSet([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) @@ -544,7 +549,13 @@ using Plots @test dim(itp) == dim(nodes) # Saving the interpolation and the function to a VTK file @test_nowarn vtk_save("itp", nodes, f, itp; keys = ["f", "itp"]) - # TODO: Check the VTK file after implementing a vtk_read function + nodes2, point_data = @test_nowarn vtk_read("itp.vtu") + @test length(nodes) == length(nodes2) + for i in eachindex(nodes) + @test [nodes[i]; 0.0] == nodes2[i] + end + @test point_data["f"] == ff + @test point_data["itp"] == itp.(nodes) @test_nowarn rm("itp.vtu", force = true) expected_coefficients = [ From 098181cb3b4d322efabb9fbfab7721f55336d87c Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 23 Apr 2024 14:16:36 +0200 Subject: [PATCH 2/2] add refs to docstrings --- src/io.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/io.jl b/src/io.jl index 05f3e3fc..fd76e4aa 100644 --- a/src/io.jl +++ b/src/io.jl @@ -2,8 +2,8 @@ vtk_save(filename, nodeset::NodeSet, functions...; keys = "value_" .* string.(eachindex(functions))) -Save a `NodeSet` to a VTK file. You can optionally pass a list of functions to save -the values of the functions at the nodes. The functions can also be passed as `Interpolation`. +Save a [`NodeSet`](@ref) to a VTK file. You can optionally pass a list of functions to save +the values of the functions at the nodes. The functions can also be passed as [`Interpolation`](@ref). The optional keyword argument `keys` is used to specify the names of the data arrays in the VTK file. """ function vtk_save(filename, nodeset::NodeSet, functions...; @@ -21,8 +21,8 @@ end """ vtk_read(filename) -Read a set of nodes from a VTK file and return them as a `NodeSet`. Note that the data -will always be returned as a 3D `NodeSet`, even if the data is 1D or 2D. The point data +Read a set of nodes from a VTK file and return them as a [`NodeSet`](@ref). Note that the data +will always be returned as a 3D [`NodeSet`](@ref), even if the data is 1D or 2D. The point data will be returned as a dictionary with the keys being the names of the data arrays in the VTK file. """ function vtk_read(filename)