Add feature to save NodeSets to VTK file (#38)

* add feature to save NodeSets to VTK file

* format

* bump compat of WriteVTK to 1.7

Project.toml

@@ -10,6 +10,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 ForwardDiff = "0.10.36"
@@ -18,4 +19,5 @@ RecipesBase = "1.1"
 SpecialFunctions = "2"
 StaticArrays = "1"
 TypedPolynomials = "0.4.1"
+WriteVTK = "1.7"
 julia = "1.8"

src/KernelInterpolation.jl

@@ -1,18 +1,20 @@
 module KernelInterpolation
 
+using ForwardDiff: ForwardDiff
 using LinearAlgebra: norm, Symmetric, tr
 using RecipesBase: RecipesBase, @recipe, @series
 using SpecialFunctions: besselk, loggamma
 using StaticArrays: StaticArrays, MVector
 using TypedPolynomials: Variable, monomials, degree
-using ForwardDiff: ForwardDiff
+using WriteVTK: vtk_grid, MeshCell, VTKCellTypes
 
 include("kernels/kernels.jl")
 include("nodes.jl")
 include("differential_operators.jl")
 include("pdes.jl")
 include("interpolation.jl")
 include("visualization.jl")
+include("io.jl")
 include("util.jl")
 
 export get_name
@@ -24,12 +26,13 @@ export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
 export phi, Phi, order
 export Laplacian
 export PoissonEquation
-export NodeSet, separation_distance, dim, values_along_dim, random_hypercube,
+export NodeSet, separation_distance, dim, eachdim, values_along_dim, random_hypercube,
        random_hypercube_boundary, homogeneous_hypercube, homogeneous_hypercube_boundary,
        random_hypersphere, random_hypersphere_boundary
 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 examples_dir, get_examples, default_example, include_example
 
 end

src/io.jl

@@ -0,0 +1,19 @@
+"""
+    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`.
+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...;
+                  keys = "value_" .* string.(eachindex(functions)))
+    @assert dim(nodeset)<=3 "Only 1D, 2D, and 3D data can be saved to VTK files."
+    cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in eachindex(nodeset)]
+    points = values_along_dim.(Ref(nodeset), eachdim(nodeset))
+    vtk_grid(filename, points..., cells, append = false) do vtk
+        for (i, func) in enumerate(functions)
+            vtk["$(keys[i])"] = func.(nodeset)
+        end
+    end
+end

src/kernels/radialsymmetric_kernel.jl

@@ -412,7 +412,7 @@ function phi(kernel::Matern12Kernel, r::Real)
     y = kernel.shape_parameter * r
     return exp(-y)
 end
-order(kernel::Matern12Kernel) = 0
+order(::Matern12Kernel) = 0
 
 @doc raw"""
     Matern32Kernel{Dim}(; shape_parameter = 1.0)

src/nodes.jl

@@ -82,15 +82,16 @@ function update_separation_distance!(nodeset::NodeSet)
         nodeset.q = q
     end
 end
-dim(nodeset::NodeSet{Dim, RealT}) where {Dim, RealT} = Dim
+dim(::NodeSet{Dim, RealT}) where {Dim, RealT} = Dim
 # Functions to treat NodeSet as array
-Base.eltype(nodeset::NodeSet{Dim, RealT}) where {Dim, RealT} = RealT
+Base.eltype(::NodeSet{Dim, RealT}) where {Dim, RealT} = RealT
 Base.length(nodeset::NodeSet) = length(nodeset.nodes)
 Base.size(nodeset::NodeSet) = (length(nodeset), dim(nodeset))
 Base.iterate(nodeset::NodeSet, state = 1) = iterate(nodeset.nodes, state)
 Base.collect(nodeset::NodeSet) = collect(nodeset.nodes)
 Base.axes(nodeset::NodeSet) = axes(nodeset.nodes)
 Base.eachindex(nodeset::NodeSet) = eachindex(nodeset.nodes)
+eachdim(nodeset::NodeSet) = Base.OneTo(dim(nodeset))
 Base.isassigned(nodeset::NodeSet, i::Int) = isassigned(nodeset.nodes, i)
 function Base.similar(nodeset::NodeSet{Dim, RealT}) where {Dim, RealT}
     NodeSet{Dim, RealT}(similar(nodeset.nodes), Inf)

test/test_unit.jl

@@ -162,6 +162,7 @@ using Plots
         for node in nodeset1
             @test node isa MVector{2, Float64}
         end
+
         f(x) = x[1] + x[2]
         ff = @test_nowarn f.(nodeset1)
         @test ff == [0.0, 1.0, 1.0, 2.0]
@@ -170,6 +171,11 @@ using Plots
         dim2 = @test_nowarn values_along_dim(nodeset1, 2)
         @test dim2 == [0.0, 0.0, 1.0, 1.0]
 
+        # Saving the nodeset to a VTK file
+        @test_nowarn vtk_save("nodeset1", nodeset1)
+        # TODO: Check the VTK file after implementing a vtk_read function
+        @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]])
         @test dim(nodeset2) == 2
         @test length(nodeset2) == 4
@@ -536,6 +542,11 @@ using Plots
         @test nodeset(itp) == nodes
         @test dim(itp) == dim(kernel)
         @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
+        @test_nowarn rm("itp.vtu", force = true)
+
         expected_coefficients = [
             -2.225451664388596,
             0.31604241814819756,