Add feature to read VTK files

.gitignore

@@ -1,2 +1,3 @@
 /docs/build/
 **/Manifest.toml
+**/*.vtu

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"

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

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...;
@@ -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`](@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)
+    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

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 = [