revise

Manifest.toml

@@ -144,9 +144,9 @@ version = "1.15.0"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
+git-tree-sha1 = "ac67408d9ddf207de5cfa9a97e114352430f01ed"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.15"
+version = "0.18.16"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -164,9 +164,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
 deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "9242eec9b7e2e14f9952e8ea1c7e31a50501d587"
+git-tree-sha1 = "a4532d110ce91bd744b99280193a317310960c46"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.104"
+version = "0.25.106"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -378,9 +378,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.LoopVectorization]]
 deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "fc712dd664097440f19a91a704299cca02134ca0"
+git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
-version = "0.12.168"
+version = "0.12.166"
 
     [deps.LoopVectorization.extensions]
     ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]

README.md

@@ -1,7 +1,7 @@
 # ChebParticleMesh
 
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://Xuanzhao Gao.github.io/ChebParticleMesh.jl/stable/)
-[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://Xuanzhao Gao.github.io/ChebParticleMesh.jl/dev/)
-[![Build Status](https://github.com/Xuanzhao Gao/ChebParticleMesh.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Xuanzhao Gao/ChebParticleMesh.jl/actions/workflows/CI.yml?query=branch%3Amain)
-[![Build Status](https://travis-ci.com/Xuanzhao Gao/ChebParticleMesh.jl.svg?branch=main)](https://travis-ci.com/Xuanzhao Gao/ChebParticleMesh.jl)
-[![Coverage](https://codecov.io/gh/Xuanzhao Gao/ChebParticleMesh.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/Xuanzhao Gao/ChebParticleMesh.jl)
+<!-- [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://HPMolSim.github.io/ChebParticleMesh.jl/stable/) -->
+<!-- [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://HPMolSim.github.io/ChebParticleMesh.jl/dev/) -->
+[![Build Status](https://github.com/HPMolSim/ChebParticleMesh.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/HPMolSim/ChebParticleMesh.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Build Status](https://travis-ci.com/HPMolSim/ChebParticleMesh.jl.svg?branch=main)](https://travis-ci.com/HPMolSim/ChebParticleMesh.jl)
+[![Coverage](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl)

src/ChebParticleMesh.jl

@@ -1,10 +1,13 @@
 module ChebParticleMesh
 
-using LinearAlgebra, ExTinyMD, FastTransforms, SpecialFunctions, LoopVectorization
+using LinearAlgebra, ExTinyMD, FastTransforms, SpecialFunctions, LoopVectorization, FFTW
 
 export horner
-export ChebCoef, pwcheb_eval, pwcheb_eval!
-export GridInfo, PadIndex, ImageIndex, id_image2pad, id_pad2image, grid_revise_image!, grid_revise_pad!, grid_image2pad!, grid_pad2image!
+export ChebCoef, pwcheb_eval, pwcheb_eval!, f_eval, f_eval!
+export GridInfo, PadIndex, ImageIndex
+export id_image2pad, id_pad2image, grid_revise_image!, grid_revise_pad!, grid_image2pad!, grid_pad2image!
+export nearest_grid_id, pad_grid_id, image_grid_id, image_grid_pos, pad_grid_pos
+export interpolate_single!, interpolate!
 
 include("types.jl")
 include("utils.jl")
@@ -13,7 +16,7 @@ include("chebpoly.jl")
 include("grid.jl")
 include("index.jl")
 
-include("interpolation.jl")
+include("interpolate.jl")
 include("scaling.jl")
 include("gather.jl")
 

src/chebpoly.jl

@@ -10,7 +10,7 @@ function ChebCoef(f::Function, h::T, w::Int, ν::Int) where{T<:Union{Float32, Fl
 """
 function ChebCoef(f::Function, h::T, w::Int, ν::Int) where{T<:Union{Float32, Float64}}
     coef = pwcheb_coef(f, h, w, ν)
-    return ChebCoef{T}(h, w, ν, coef)
+    return ChebCoef{T}(f, h, w, ν, coef)
 end
 
 function pwcheb_coef(f::Function, h::T, window_cut::Int, ν::Int) where{T<:Union{Float32, Float64}}
@@ -50,7 +50,7 @@ function pwcheb_eval!(x0::T, cheb_approx::Array{T, 1}, chebcoef::ChebCoef{T}) wh
 
     h = chebcoef.h
     C = chebcoef.coef
-    @assert abs(x0) ≤ h / 2
+    @assert (abs(x0) ≤ h / 2) | (abs(x0) ≈ h / 2)
     @assert size(cheb_approx, 1) == size(C, 1) - 1
 
     if x0 ≥ zero(T)
@@ -73,5 +73,29 @@ function pwcheb_eval(x0::T, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float
     cheb_approx = zeros(T, size(chebcoef.coef, 1) - 1)
     cheb_approx = pwcheb_eval!(x0, cheb_approx, chebcoef)
 
+    return cheb_approx
+end
+
+function f_eval!(x0::T, cheb_approx::Array{T, 1}, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float64}}
+
+    h = chebcoef.h
+    C = chebcoef.coef
+    w = chebcoef.w
+    @assert (abs(x0) ≤ h / 2) | (abs(x0) ≈ h / 2)
+    @assert size(cheb_approx, 1) == size(C, 1) - 1
+
+    for i in - w : w
+        x = - x0 + i * h
+        cheb_approx[i + w + 1] = chebcoef.f(x)
+    end
+
+    return cheb_approx
+end
+
+function f_eval(x0::T, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float64}}
+
+    cheb_approx = zeros(T, 2 * chebcoef.w + 1)
+    cheb_approx = f_eval!(x0, cheb_approx, chebcoef)
+
     return cheb_approx
 end

src/grid.jl

@@ -21,13 +21,28 @@ function GridInfo(N_real::NTuple{3, Int}, w::NTuple{3, Int}, periodicity::NTuple
 
     trans_info = TransInfo(N_real, periodicity, image, pad)
 
-    return GridInfo{T}(N_real, w, periodicity, image, pad, L, h, N_image, N_pad, trans_info)
+    k = Vector{Array{T, 1}}()
+
+    for i in 1:3
+        ki = [2π * j / (L[i] + h[i] * 2 * pad[i]) for j in 0:(N_pad[i] - 1)]
+
+        for j in 1:N_pad[i]
+            if j - 1 > ceil(N_pad[i] / 2)
+                ki[j] -= 2π * N_pad[i] / (L[i] + h[i] * 2 * pad[i])
+            end
+        end
+
+        push!(k, ki)
+    end
+
+    return GridInfo{T}(N_real, w, periodicity, image, pad, L, h, N_image, N_pad, trans_info, k)
 end
 
 function GridBox(grid_info::GridInfo{T}) where{T<:Union{Float32, Float64}}
     pad_grid = zeros(T, grid_info.N_pad...)
     image_grid = zeros(T, grid_info.N_image...)
-    return GridBox{T}(pad_grid, image_grid)
+    cheb_value = [zeros(T, 2 * grid_info.w[i] + 1) for i in 1:3]
+    return GridBox{T}(pad_grid, image_grid, cheb_value)
 end
 
 """

src/index.jl

@@ -21,17 +21,35 @@ function pad_grid_id(x::Point{3, T}, gridinfo::GridInfo{T}) where{T}
     return PadIndex(idx, idy, idz)
 end
 
+function pad_grid_pos(id::PadIndex, gridinfo::GridInfo{T}) where{T}
+
+    x = (id.idx - gridinfo.pad[1] - T(0.5)) * gridinfo.h[1]
+    y = (id.idy - gridinfo.pad[2] - T(0.5)) * gridinfo.h[2]
+    z = (id.idz - gridinfo.pad[3] - T(0.5)) * gridinfo.h[3]
+
+    return Point(x, y, z)
+end
+
 function image_grid_id(x::Point{3, T}, gridinfo::GridInfo{T}) where{T}
 
     idx, idy, idz = nearest_grid_id(x, gridinfo.h)
 
-    idx = idx + gridinfo.w[1] + 1
-    idy = idy + gridinfo.w[2] + 1
-    idz = idz + gridinfo.w[3] + 1
+    idx = idx + gridinfo.image[1]
+    idy = idy + gridinfo.image[2]
+    idz = idz + gridinfo.image[3]
 
     return ImageIndex(idx, idy, idz)
 end
 
+function image_grid_pos(id::ImageIndex, gridinfo::GridInfo{T}) where{T}
+
+    x = (id.idx - gridinfo.image[1] - T(0.5)) * gridinfo.h[1]
+    y = (id.idy - gridinfo.image[2] - T(0.5)) * gridinfo.h[2]
+    z = (id.idz - gridinfo.image[3] - T(0.5)) * gridinfo.h[3]
+
+    return Point(x, y, z)
+end
+
 function id_image2pad_single(id_image_i::Int, N_i::T, image_i::Int, pad_i::Int, periodicity::Bool) where{T}
 
     id_pad_i = 0

src/interpolate.jl

@@ -0,0 +1,79 @@
+function interpolate_single!(q::T, pos::Point{3, T}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T}
+
+    cheb_value = gridbox.cheb_value
+    near_id_image = image_grid_id(pos, gridinfo)
+    near_pos_image = image_grid_pos(near_id_image, gridinfo)
+    for i in 1:3
+        dx = pos[i] - near_pos_image[i]
+        pwcheb_eval!(dx, cheb_value[i], chebcoefs[i])
+    end
+
+    @turbo for i in - gridinfo.w[1] : gridinfo.w[1]
+        for j in - gridinfo.w[2] : gridinfo.w[2]
+            for k in - gridinfo.w[3] : gridinfo.w[3]
+                gridbox.image_grid[
+                    near_id_image.idx + i, 
+                    near_id_image.idy + j, 
+                    near_id_image.idz + k] += 
+                    q * 
+                    cheb_value[1][i + gridinfo.w[1] + 1] * 
+                    cheb_value[2][j + gridinfo.w[2] + 1] * 
+                    cheb_value[3][k + gridinfo.w[3] + 1]
+            end
+        end
+    end
+
+    return nothing
+end
+
+function interpolate!(qs::Vector{T}, poses::Vector{Point{3, T}}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T}
+
+    @assert length(qs) == length(poses)
+    grid_revise_image!(gridbox)
+
+    for i in 1:length(qs)
+        interpolate_single!(qs[i], poses[i], gridinfo, gridbox, chebcoefs)
+    end
+
+    return nothing
+end
+
+function interpolate_single_direct!(q::T, pos::Point{3, T}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T}
+
+    cheb_value = gridbox.cheb_value
+    near_id_image = image_grid_id(pos, gridinfo)
+    near_pos_image = image_grid_pos(near_id_image, gridinfo)
+    for i in 1:3
+        dx = pos[i] - near_pos_image[i]
+        f_eval!(dx, cheb_value[i], chebcoefs[i])
+    end
+
+    @turbo for i in - gridinfo.w[1] : gridinfo.w[1]
+        for j in - gridinfo.w[2] : gridinfo.w[2]
+            for k in - gridinfo.w[3] : gridinfo.w[3]
+                gridbox.image_grid[
+                    near_id_image.idx + i, 
+                    near_id_image.idy + j, 
+                    near_id_image.idz + k] += 
+                    q * 
+                    cheb_value[1][i + gridinfo.w[1] + 1] * 
+                    cheb_value[2][j + gridinfo.w[2] + 1] * 
+                    cheb_value[3][k + gridinfo.w[3] + 1]
+            end
+        end
+    end
+
+    return nothing
+end
+
+function interpolate_direct!(qs::Vector{T}, poses::Vector{Point{3, T}}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T}
+
+    @assert length(qs) == length(poses)
+    grid_revise_image!(gridbox)
+
+    for i in 1:length(qs)
+        interpolate_single_direct!(qs[i], poses[i], gridinfo, gridbox, chebcoefs)
+    end
+
+    return nothing
+end

src/types.jl

@@ -1,4 +1,5 @@
 struct ChebCoef{T}
+    f::Function
     h::T
     w::Int
     ν::Int
@@ -18,11 +19,15 @@ struct GridInfo{T}
     N_pad::NTuple{3, Int}
 
     trans_info::Vector{Tuple{NTuple{3, Int}, NTuple{3, Int}, NTuple{3, Int}}}
+
+    k::Vector{Array{T, 1}}
 end
 
 mutable struct GridBox{T}
     pad_grid::Array{T, 3}
     image_grid::Array{T, 3}
+
+    cheb_value::Vector{Array{T, 1}}
 end
 
 abstract type AbstractIndex end

test/chebpoly.jl

@@ -6,6 +6,6 @@
     c = ChebCoef(f, h, w, ν)
     for x0 in [- h / 2 : 0.01 : h / 2...]
         p = [i * h for i in -w:w]
-        @test pwcheb_eval(x0, c) ≈ f.(p .- x0)
+        @test pwcheb_eval(x0, c) ≈ f_eval(x0, c) ≈ f.(p .- x0)
     end
 end

test/index.jl

@@ -1,4 +1,37 @@
-@testset "grid transform: image to pad" begin
+@testset "pos/id transform" begin
+    N_real = (20, 30, 40)
+    w = (12, 23, 3)
+    extra_pad = (24, 15, 6)
+    L = (100.0, 100.0, 100.0)
+
+    for periodicity_x in [true, false]
+        for periodicity_y in [true, false]
+            for periodicity_z in [true, false]
+                periodicity = (periodicity_x, periodicity_y, periodicity_z)
+                @testset "peridoicity is $(periodicity)" begin
+                    gridinfo = GridInfo(N_real, w, periodicity, extra_pad, L)
+                    for _ = 1:10
+                        pos = Point(L[1] * rand(), L[2] * rand(), L[3] * rand())
+
+                        id_image = image_grid_id(pos, gridinfo)
+                        pos_trans_image = image_grid_pos(id_image, gridinfo)
+                        for i in 1:3
+                            @test abs(pos[i] - pos_trans_image[i]) ≤ gridinfo.h[i] / 2
+                        end
+
+                        id_pad = pad_grid_id(pos, gridinfo)
+                        pos_trans_pad = pad_grid_pos(id_pad, gridinfo)
+                        for i in 1:3
+                            @test abs(pos[i] - pos_trans_pad[i]) ≤ gridinfo.h[i] / 2
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+
+@testset "grid transform" begin
     N_real = (5, 6, 7)
     w = (1, 3, 4)
     extra_pad = (4, 5, 6)

test/interpolate.jl

@@ -0,0 +1,5 @@
+@testset "Chebyshev interpolation vs direct interpolation" begin
+    f = x -> exp(-x^2 * 10.0)
+
+
+end

test/runtests.jl

@@ -1,8 +1,10 @@
 using ChebParticleMesh
 using Test
+using ExTinyMD
 
 @testset "ChebParticleMesh.jl" begin
     include("utils.jl")
     include("chebpoly.jl")
     include("index.jl")
+    include("interpolate.jl")
 end