Use mesh to construct MPValue

docs/literate/examples/dam_break.jl

@@ -92,9 +92,9 @@ function main(transfer::Transfer = FLIP(1.0))
 
     ## Interpolation
     it = KernelCorrection(QuadraticBSpline())
-    mpvalues = map(p -> MPValue(Vec{2}, it), eachindex(particles))
+    mpvalues = map(p -> MPValue(it, grid.X), eachindex(particles))
     mpvalues_cell = map(CartesianIndices(size(grid).-1)) do cell
-        mp = MPValue(Vec{2}, it)
+        mp = MPValue(it, grid.X)
         xc = cellcenter(cell, grid.X)
         update!(mp, xc, grid.X)
     end

docs/literate/examples/elastic_impact.jl

@@ -87,7 +87,7 @@ function main(transfer::Transfer = FLIP(1.0))
     @show length(particles)
 
     ## Interpolation
-    mpvalues = generate_mpvalues(Vec{2, Float64}, QuadraticBSpline(), length(particles))
+    mpvalues = generate_mpvalues(QuadraticBSpline(), grid.x, length(particles))
 
     ## Material model (neo-Hookean)
     function caucy_stress(F)

docs/literate/examples/implicit_jacobian_based.jl

@@ -67,7 +67,7 @@ function main()
     ## Interpolation
     ## Use the kernel correction to properly handle the boundary conditions
     it = KernelCorrection(QuadraticBSpline())
-    mpvalues = generate_mpvalues(Vec{3}, it, length(particles))
+    mpvalues = generate_mpvalues(it, grid.X, length(particles))
 
     ## Neo-Hookean model
     function kirchhoff_stress(F)

docs/literate/examples/implicit_jacobian_free.jl

@@ -69,7 +69,7 @@ function main()
 
     ## Interpolation
     ## Use the kernel correction to properly handle the boundary conditions
-    mpvalues = generate_mpvalues(Vec{3}, KernelCorrection(QuadraticBSpline()), length(particles))
+    mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.X, length(particles))
 
     ## Neo-Hookean model
     function kirchhoff_stress(F)

docs/literate/examples/rigid_body_contact.jl

@@ -97,7 +97,7 @@ function main()
     @show length(particles)
 
     ## Interpolation
-    mpvalues = generate_mpvalues(Vec{2}, KernelCorrection(QuadraticBSpline()), length(particles))
+    mpvalues = generate_mpvalues(KernelCorrection(QuadraticBSpline()), grid.x, length(particles))
 
     ## Output
     outdir = mkpath(joinpath("output", "rigid_body_contact"))

docs/literate/examples/tlmpm_vortex.jl

@@ -84,7 +84,7 @@ function main()
     @show length(particles)
 
     ## Precompute linear kernel values
-    mpvalues = generate_mpvalues(Vec{2, Float64}, LinearBSpline(), length(particles))
+    mpvalues = generate_mpvalues(LinearBSpline(), grid.X, length(particles))
     for p in eachindex(particles, mpvalues)
         update!(mpvalues[p], particles.x[p], grid.X)
     end

src/Interpolations/mpvalue.jl

@@ -1,7 +1,14 @@
-abstract type Interpolation end
 abstract type Kernel <: Interpolation end
 
-function create_property(::Type{Vec{dim, T}}, it::Interpolation; diff=gradient) where {dim, T}
+#=
+To create a new interpolation, following methods need to be implemented.
+* Tesserae.create_property(::Type{T}, it::Interpolation, mesh; kwargs...) -> NamedTuple
+* Tesserae.update_property!(mp::MPValue, it::Interpolation, pt, mesh)
+* Tesserae.neighboringnodes(it::Interpolation, pt, mesh)
+* Tesserae.NeighboringNodesType(it::Interpolation, mesh)
+=#
+
+function create_property(::Type{T}, it::Interpolation, mesh::CartesianMesh{dim}; diff=gradient) where {dim, T}
     dims = nfill(gridspan(it), Val(dim))
     (diff === nothing || diff === identity) && return (; w=zeros(T, dims))
     diff === gradient && return (; w=fill(zero(T), dims), ∇w=fill(zero(Vec{dim, T}), dims))
@@ -46,17 +53,16 @@ struct MPValue{It, Prop <: NamedTuple, Indices <: AbstractArray{<: Any, 0}}
     indices::Indices
 end
 
-function _MPValue(it::Union{Nothing, Interpolation}, prop::NamedTuple)
-    @assert hasproperty(prop, :w)
-    @assert prop.w isa AbstractArray{<: Real}
-    dim = ndims(prop.w)
-    indices = fill(EmptyCartesianIndices(Val(dim)))
-    MPValue(it, prop, indices)
+function MPValue(it::Interpolation, prop::NamedTuple, ::Type{Tinds}) where {Tinds}
+    MPValue(it, prop, Array{Tinds}(undef))
 end
 
-MPValue(prop::NamedTuple) = _MPValue(nothing, prop)
-MPValue(::Type{Vec{dim, T}}, it::Interpolation; kwargs...) where {dim, T} = _MPValue(it, create_property(Vec{dim, T}, it; kwargs...))
-MPValue(::Type{Vec{dim}}, it::Interpolation; kwargs...) where {dim} = MPValue(Vec{dim, Float64}, it; kwargs...)
+function MPValue(::Type{T}, it::Interpolation, mesh::AbstractMesh; kwargs...) where {T}
+    prop = create_property(T, it, mesh; kwargs...)
+    Tinds = NeighboringNodesType(it, mesh)
+    MPValue(it, prop, Tinds)
+end
+MPValue(it::Interpolation, mesh::AbstractMesh; kwargs...) = MPValue(Float64, it, mesh; kwargs...)
 
 Base.propertynames(mp::MPValue) = propertynames(getfield(mp, :prop))
 @inline function Base.getproperty(mp::MPValue, name::Symbol)
@@ -128,20 +134,18 @@ struct MPValueVector{It, Prop <: NamedTuple, Indices, ElType <: MPValue{It}} <:
     indices::Indices
 end
 
-function generate_mpvalues(::Type{Vec{dim, T}}, it::Interpolation, n::Int; kwargs...) where {dim, T}
-    prop = map(create_property(Vec{dim, T}, it; kwargs...)) do prop
+function generate_mpvalues(::Type{T}, it::Interpolation, mesh::AbstractMesh, n::Int; kwargs...) where {T}
+    prop = map(create_property(T, it, mesh; kwargs...)) do prop
         fill(zero(eltype(prop)), size(prop)..., n)
     end
-    indices = fill(EmptyCartesianIndices(Val(dim)), n)
+    indices = Array{NeighboringNodesType(it, mesh)}(undef, n)
     It = typeof(it)
     Prop = typeof(prop)
     Indices = typeof(indices)
     ElType = Base._return_type(_getindex, Tuple{It, Prop, Indices, Int})
     MPValueVector{It, Prop, Indices, ElType}(it, prop, indices)
 end
-generate_mpvalues(::Type{Vec{dim}}, it::Interpolation, n::Int; diff=gradient) where {dim} = generate_mpvalues(Vec{dim, Float64}, it, n; diff)
-
-generate_mpvalues(::Type{V}, it::Interpolation, n::Int) where {V} = generate_mpvalues(V, it, Val(1), n)
+generate_mpvalues(it::Interpolation, mesh::AbstractMesh, n::Int; kwargs...) = generate_mpvalues(Float64, it, mesh, n; kwargs...)
 
 Base.IndexStyle(::Type{<: MPValueVector}) = IndexLinear()
 Base.size(x::MPValueVector) = size(getfield(x, :indices))
@@ -198,17 +202,16 @@ end
     true
 end
 
-function update!(mp::MPValue, it::Interpolation, pt, mesh::AbstractMesh)
-    set_neighboringnodes!(mp, neighboringnodes(interpolation(mp), pt, mesh))
+function update!(mp::MPValue, pt, mesh::AbstractMesh)
+    it = interpolation(mp)
+    set_neighboringnodes!(mp, neighboringnodes(it, pt, mesh))
     update_property!(mp, it, pt, mesh)
     mp
 end
-function update!(mp::MPValue, it::Interpolation, pt, mesh::AbstractMesh, filter::AbstractArray{Bool})
+function update!(mp::MPValue, pt, mesh::AbstractMesh, filter::AbstractArray{Bool})
     @debug @assert size(mesh) == size(filter)
-    set_neighboringnodes!(mp, neighboringnodes(interpolation(mp), pt, mesh))
+    it = interpolation(mp)
+    set_neighboringnodes!(mp, neighboringnodes(it, pt, mesh))
     update_property!(mp, it, pt, mesh, filter)
     mp
 end
-
-update!(mp::MPValue{<: Interpolation}, pt, mesh::AbstractMesh) = update!(mp, interpolation(mp), pt, mesh)
-update!(mp::MPValue{<: Interpolation}, pt, mesh::AbstractMesh, filter::AbstractArray{Bool}) = update!(mp, interpolation(mp), pt, mesh, filter)

src/Tesserae.jl

@@ -82,6 +82,8 @@ export
     closevtm,
     closepvd
 
+abstract type Interpolation end
+
 include("utils.jl")
 include("mesh.jl")
 include("blockspace.jl")

src/mesh.jl

@@ -137,6 +137,8 @@ CartesianIndices((1:5,))
 end
 @inline EmptyCartesianIndices(::Val{dim}) where {dim} = CartesianIndices(nfill(1:0, Val(dim)))
 
+NeighboringNodesType(::Interpolation, ::CartesianMesh{dim}) where {dim} = CartesianIndices{dim, NTuple{dim, UnitRange{Int}}}
+
 """
     whichcell(x::Vec, mesh::CartesianMesh)
 

test/interpolations.jl

@@ -1,57 +1,59 @@
 @testset "MPValue" begin
-    for dim in (1,2,3)
-        @test eltype(MPValue(Vec{dim}, QuadraticBSpline())) ==
-              eltype(MPValue(Vec{dim, Float64}, QuadraticBSpline()))
-        for T in (Float32, Float64)
-            for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline(), uGIMP())
-                for extension in (identity, KernelCorrection)
-                    it = extension(kernel)
-                    mp = @inferred MPValue(Vec{dim,T}, it)
-                    @test Tesserae.interpolation(mp) === it
-                    @test mp.w isa Array{T}
-                    @test mp.∇w isa Array{Vec{dim,T}}
-                    @test ndims(mp.w) == dim
-                    @test ndims(mp.∇w) == dim
-                    @test size(mp.w) == size(mp.∇w)
-                    @test all(size(neighboringnodes(mp)) .< size(mp.w))
+    @testset "CartesianMesh" begin
+        for dim in (1,2,3)
+            mesh = CartesianMesh(1, ntuple(d->(0,10), dim)...)
+            for T in (Float32, Float64)
+                for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline(), uGIMP())
+                    for extension in (identity, KernelCorrection)
+                        it = extension(kernel)
+                        mp = @inferred MPValue(T, it, mesh)
+                        @test Tesserae.interpolation(mp) === it
+                        @test mp.w isa Array{T}
+                        @test mp.∇w isa Array{Vec{dim,T}}
+                        @test ndims(mp.w) == dim
+                        @test ndims(mp.∇w) == dim
+                        @test size(mp.w) == size(mp.∇w)
+                        @test typeof(neighboringnodes(mp)) === CartesianIndices{dim, NTuple{dim, UnitRange{Int}}}
 
-                    # diff = nothing
-                    mp = @inferred MPValue(Vec{dim,T}, it; diff=nothing)
-                    @test hasproperty(mp, :w) && mp.w isa Array{T}
-                    @test !hasproperty(mp, :∇w)
-                    @test !hasproperty(mp, :∇∇w)
-                    @test ndims(mp.w) == dim
-                    # diff = gradient
-                    mp = @inferred MPValue(Vec{dim,T}, it; diff=gradient)
-                    @test hasproperty(mp, :w)  && mp.w  isa Array{T}
-                    @test hasproperty(mp, :∇w) && mp.∇w isa Array{Vec{dim,T}}
-                    @test !hasproperty(mp, :∇∇w)
-                    @test ndims(mp.w) == ndims(mp.∇w) == dim
-                    @test size(mp.w) == size(mp.∇w)
-                    # diff = hessian
-                    mp = @inferred MPValue(Vec{dim,T}, it; diff=hessian)
-                    @test hasproperty(mp, :w)   && mp.w   isa Array{T}
-                    @test hasproperty(mp, :∇w)  && mp.∇w  isa Array{Vec{dim,T}}
-                    @test hasproperty(mp, :∇∇w) && mp.∇∇w isa Array{<: SymmetricSecondOrderTensor{dim,T}}
-                    @test ndims(mp.w) == ndims(mp.∇w) == ndims(mp.∇∇w) == dim
-                    @test size(mp.w) == size(mp.∇w) == size(mp.∇∇w)
+                        # diff = nothing
+                        mp = @inferred MPValue(T, it, mesh; diff=nothing)
+                        @test hasproperty(mp, :w) && mp.w isa Array{T}
+                        @test !hasproperty(mp, :∇w)
+                        @test !hasproperty(mp, :∇∇w)
+                        @test ndims(mp.w) == dim
+                        # diff = gradient
+                        mp = @inferred MPValue(T, it, mesh; diff=gradient)
+                        @test hasproperty(mp, :w)  && mp.w  isa Array{T}
+                        @test hasproperty(mp, :∇w) && mp.∇w isa Array{Vec{dim,T}}
+                        @test !hasproperty(mp, :∇∇w)
+                        @test ndims(mp.w) == ndims(mp.∇w) == dim
+                        @test size(mp.w) == size(mp.∇w)
+                        # diff = hessian
+                        mp = @inferred MPValue(T, it, mesh; diff=hessian)
+                        @test hasproperty(mp, :w)   && mp.w   isa Array{T}
+                        @test hasproperty(mp, :∇w)  && mp.∇w  isa Array{Vec{dim,T}}
+                        @test hasproperty(mp, :∇∇w) && mp.∇∇w isa Array{<: SymmetricSecondOrderTensor{dim,T}}
+                        @test ndims(mp.w) == ndims(mp.∇w) == ndims(mp.∇∇w) == dim
+                        @test size(mp.w) == size(mp.∇w) == size(mp.∇∇w)
+                    end
                 end
             end
         end
     end
 end
 
 @testset "MPValueVector" begin
-    for dim in (1,2,3)
-        @test eltype(generate_mpvalues(Vec{dim}, QuadraticBSpline(), 2)) ==
-              eltype(generate_mpvalues(Vec{dim, Float64}, QuadraticBSpline(), 2))
-        for T in (Float32, Float64)
-            n = 100
-            mpvalues = @inferred generate_mpvalues(Vec{dim,T}, QuadraticBSpline(), n)
-            @test size(mpvalues) === (n,)
-            @test Tesserae.interpolation(mpvalues) === QuadraticBSpline()
-            @test all(eachindex(mpvalues)) do i
-                typeof(mpvalues[1]) === eltype(mpvalues)
+    @testset "CartesianMesh" begin
+        for dim in (1,2,3)
+            mesh = CartesianMesh(1, ntuple(d->(0,10), dim)...)
+            for T in (Float32, Float64)
+                n = 100
+                mpvalues = @inferred generate_mpvalues(T, QuadraticBSpline(), mesh, n)
+                @test size(mpvalues) === (n,)
+                @test Tesserae.interpolation(mpvalues) === QuadraticBSpline()
+                @test all(eachindex(mpvalues)) do i
+                    typeof(mpvalues[i]) === eltype(mpvalues)
+                end
             end
         end
     end
@@ -73,10 +75,10 @@ end
     end
 
     @testset "$it" for it in (LinearBSpline(), QuadraticBSpline(), CubicBSpline())
-        for T in (Float32, Float64), dim in (1,2,3)
+        for dim in (1,2,3)
             Random.seed!(1234)
-            mp = MPValue(Vec{dim}, it)
             mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...)
+            mp = MPValue(it, mesh)
             @test all(1:100) do _
                 x = rand(Vec{dim})
                 update!(mp, x, mesh)
@@ -91,8 +93,8 @@ end
     @testset "$it" for it in (SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline())
         for dim in (1,2,3)
             Random.seed!(1234)
-            mp = MPValue(Vec{dim}, it)
             mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...)
+            mp = MPValue(it, mesh)
             @test all(1:100) do _
                 x = rand(Vec{dim})
                 update!(mp, x, mesh)
@@ -108,8 +110,8 @@ end
         it = uGIMP()
         for dim in (1,2,3)
             Random.seed!(1234)
-            mp = MPValue(Vec{dim}, it)
             mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...)
+            mp = MPValue(it, mesh)
             l = 0.5*spacing(mesh) / 2
             @test all(1:100) do _
                 x = rand(Vec{dim})
@@ -128,8 +130,8 @@ end
         it = Wrapper(kernel)
         for dim in (1,2,3)
             Random.seed!(1234)
-            mp = MPValue(Vec{dim}, it)
             mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...)
+            mp = MPValue(it, mesh)
             l = 0.5*spacing(mesh) / 2
             @test all(1:100) do _
                 x = rand(Vec{dim})
@@ -183,11 +185,11 @@ end
         j = findfirst(==(i), neighboringnodes(mp))
         j === nothing ? zero(eltype(mp.w)) : mp.w[j]
     end
-    function kernelvalues(mesh::CartesianMesh{dim, T}, kernel, poly, index::CartesianIndex{dim}) where {dim, T}
-        mp = MPValue(Vec{dim}, KernelCorrection(kernel, poly))
+    function kernelvalues(mesh::CartesianMesh{dim}, kernel, poly, index::CartesianIndex{dim}) where {dim}
+        mp = MPValue(KernelCorrection(kernel, poly), mesh)
         L = kernel isa QuadraticBSpline ? 1.5 :
             kernel isa CubicBSpline     ? 2.0 : error()
-        X = ntuple(i -> range(max(mesh[1][i],index[i]-L-1), min(mesh[end][i],index[i]+L-1)-sqrt(eps(T)), step=0.04), Val(dim))
+        X = ntuple(i -> range(max(mesh[1][i],index[i]-L-1), min(mesh[end][i],index[i]+L-1)-sqrt(eps(Float64)), step=0.04), Val(dim))
         Z = Array{Float64}(undef, length.(X))
         for i in CartesianIndices(Z)
             @inbounds Z[i] = kernelvalue(mp, Vec(map(getindex, X, Tuple(i))), mesh, index)