From 769bc019df5538de00991b84b445543d08cb710c Mon Sep 17 00:00:00 2001 From: Keita Nakamura Date: Mon, 5 Aug 2024 16:23:54 +0900 Subject: [PATCH] Use mesh to construct MPValue --- docs/literate/examples/dam_break.jl | 4 +- docs/literate/examples/elastic_impact.jl | 2 +- .../examples/implicit_jacobian_based.jl | 2 +- .../examples/implicit_jacobian_free.jl | 2 +- docs/literate/examples/rigid_body_contact.jl | 2 +- docs/literate/examples/tlmpm_vortex.jl | 2 +- src/Interpolations/mpvalue.jl | 51 +++++---- src/Tesserae.jl | 2 + src/mesh.jl | 2 + test/interpolations.jl | 108 +++++++++--------- 10 files changed, 93 insertions(+), 84 deletions(-) diff --git a/docs/literate/examples/dam_break.jl b/docs/literate/examples/dam_break.jl index 461cd43e..47743509 100644 --- a/docs/literate/examples/dam_break.jl +++ b/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 diff --git a/docs/literate/examples/elastic_impact.jl b/docs/literate/examples/elastic_impact.jl index 71102af0..d6496470 100644 --- a/docs/literate/examples/elastic_impact.jl +++ b/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) diff --git a/docs/literate/examples/implicit_jacobian_based.jl b/docs/literate/examples/implicit_jacobian_based.jl index 470a3516..4d13a70c 100644 --- a/docs/literate/examples/implicit_jacobian_based.jl +++ b/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) diff --git a/docs/literate/examples/implicit_jacobian_free.jl b/docs/literate/examples/implicit_jacobian_free.jl index d06ad7cd..aef0d4ca 100644 --- a/docs/literate/examples/implicit_jacobian_free.jl +++ b/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) diff --git a/docs/literate/examples/rigid_body_contact.jl b/docs/literate/examples/rigid_body_contact.jl index 32e85569..7c434039 100644 --- a/docs/literate/examples/rigid_body_contact.jl +++ b/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")) diff --git a/docs/literate/examples/tlmpm_vortex.jl b/docs/literate/examples/tlmpm_vortex.jl index d2f3bff0..30f23777 100644 --- a/docs/literate/examples/tlmpm_vortex.jl +++ b/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 diff --git a/src/Interpolations/mpvalue.jl b/src/Interpolations/mpvalue.jl index 5722d246..050b229e 100644 --- a/src/Interpolations/mpvalue.jl +++ b/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) diff --git a/src/Tesserae.jl b/src/Tesserae.jl index 81b47fdf..4f1a3b0c 100644 --- a/src/Tesserae.jl +++ b/src/Tesserae.jl @@ -82,6 +82,8 @@ export closevtm, closepvd +abstract type Interpolation end + include("utils.jl") include("mesh.jl") include("blockspace.jl") diff --git a/src/mesh.jl b/src/mesh.jl index bc695439..668554de 100644 --- a/src/mesh.jl +++ b/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) diff --git a/test/interpolations.jl b/test/interpolations.jl index 61590a34..b2890e2d 100644 --- a/test/interpolations.jl +++ b/test/interpolations.jl @@ -1,40 +1,41 @@ @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 @@ -42,16 +43,17 @@ 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)