diff --git a/Project.toml b/Project.toml index bbe4b9fe..b26afcb9 100644 --- a/Project.toml +++ b/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" diff --git a/src/KernelInterpolation.jl b/src/KernelInterpolation.jl index 029b5114..29f60aec 100644 --- a/src/KernelInterpolation.jl +++ b/src/KernelInterpolation.jl @@ -1,11 +1,12 @@ 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") @@ -13,6 +14,7 @@ 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 diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 00000000..5b979cab --- /dev/null +++ b/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 diff --git a/src/kernels/radialsymmetric_kernel.jl b/src/kernels/radialsymmetric_kernel.jl index c6d6ad9e..6c166a70 100644 --- a/src/kernels/radialsymmetric_kernel.jl +++ b/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) diff --git a/src/nodes.jl b/src/nodes.jl index 5d5e24be..778d0a27 100644 --- a/src/nodes.jl +++ b/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) diff --git a/test/test_unit.jl b/test/test_unit.jl index c17ef043..ee253ed3 100644 --- a/test/test_unit.jl +++ b/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,