Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Matrix-valued interpolations #958

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Matrix-valued interpolations #958

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ee30493
Matrix-valued interpolation from Lagrange.
termi-official Oct 19, 2023
4fa7df4
Merge master and make sure it works in principle.
termi-official May 28, 2024
c8a759a
[no ci] ws
termi-official May 28, 2024
4e1d008
Test works up to transpose. :)
termi-official May 28, 2024
9e95261
Fix matrix shape indices.
termi-official May 28, 2024
7ff5882
Rename.
termi-official May 28, 2024
6312f56
Revert docstring oopsie.
termi-official May 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions src/FEValues/FunctionValues.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,16 @@
typeof_dNdx(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, sdim, T}
typeof_dNdξ(::Type{T}, ::VectorInterpolation{vdim}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim, sdim, rdim} = SMatrix{vdim, rdim, T}

# Matrix vdim == sdim == rdim
typeof_N( ::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
typeof_dNdx(::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
typeof_dNdξ(::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}

# Matrix vdim != sdim != rdim
typeof_N( ::Type{T}, ::MatrixizedInterpolation{vdim1,vdim2}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim1, vdim2, sdim, rdim} = SMatrix{vdim1, vdim2, T}
typeof_dNdx(::Type{T}, ::MatrixizedInterpolation{vdim1,vdim2}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim1, vdim2, sdim, rdim} = SArray{Tuple{vdim1, vdim2, sdim}, T, 3, vdim1*vdim2*sdim}
typeof_dNdξ(::Type{T}, ::MatrixizedInterpolation{vdim1,vdim2}, ::VectorizedInterpolation{sdim, <: AbstractRefShape{rdim}}) where {T, vdim1, vdim2, sdim, rdim} = SArray{Tuple{vdim1, vdim2, rdim}, T, 3, vdim1*vdim2*sdim}

Check warning on line 36 in src/FEValues/FunctionValues.jl

View check run for this annotation

Codecov / codecov/patch

src/FEValues/FunctionValues.jl#L34-L36 

Added lines #L34 - L36 were not covered by tests

"""
FunctionValues{DiffOrder}(::Type{T}, ip_fun, qr::QuadratureRule, ip_geo::VectorizedInterpolation)

Expand Down Expand Up @@ -98,6 +108,7 @@
sdim_from_gradtype(::Type{<:AbstractTensor{<:Any,sdim}}) where sdim = sdim
sdim_from_gradtype(::Type{<:SVector{sdim}}) where sdim = sdim
sdim_from_gradtype(::Type{<:SMatrix{<:Any,sdim}}) where sdim = sdim
sdim_from_gradtype(::Type{<:SArray{<:Tuple{<:Any,<:Any,sdim}}}) where sdim = sdim

Check warning on line 111 in src/FEValues/FunctionValues.jl

View check run for this annotation

Codecov / codecov/patch

src/FEValues/FunctionValues.jl#L111 

Added line #L111 was not covered by tests

# For performance, these must be fully inferable for the compiler.
# args: valname (:CellValues or :FacetValues), shape_gradient_type, eltype(x)
Expand Down
2 changes: 1 addition & 1 deletion src/Ferrite.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ using OrderedCollections:
using SparseArrays:
SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse
using StaticArrays:
StaticArrays, SMatrix, SVector
StaticArrays, SArray, SMatrix, SVector
using WriteVTK:
WriteVTK, VTKCellTypes
using Tensors:
Expand Down
125 changes: 108 additions & 17 deletions src/interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@

const InterpolationByDim{dim} = Interpolation{<:AbstractRefShape{dim}}

abstract type ScalarInterpolation{ refshape, order} <: Interpolation{refshape, order, Nothing} end
abstract type VectorInterpolation{vdim, refshape, order} <: Interpolation{refshape, order, Nothing} end
abstract type ScalarInterpolation{ refshape, order} <: Interpolation{refshape, order, Nothing} end
abstract type VectorInterpolation{vdim, refshape, order} <: Interpolation{refshape, order, Nothing} end
abstract type MatrixInterpolation{vdim1, vdim2, refshape, order} <: Interpolation{refshape, order, Nothing} end

# Number of components for the interpolation.
n_components(::ScalarInterpolation) = 1
Expand Down Expand Up @@ -225,8 +226,8 @@
"""
shape_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int)

Optimized version combining the evaluation [`Ferrite.shape_value(::Interpolation)`](@ref)
and [`Ferrite.shape_gradient(::Interpolation)`](@ref).
Optimized version combining the evaluation [`shape_value(::Interpolation)`](@ref)
and [`shape_gradient(::Interpolation)`](@ref).
"""
function shape_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int)
return gradient(x -> shape_value(ip, x, i), ξ, :all)
Expand All @@ -236,8 +237,8 @@
"""
shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int)

Optimized version combining the evaluation [`Ferrite.shape_value(::Interpolation)`](@ref),
[`Ferrite.shape_gradient(::Interpolation)`](@ref), and the gradient of the latter.
Optimized version combining the evaluation [`shape_value(::Interpolation)`](@ref),
[`shape_gradient(::Interpolation)`](@ref), and the gradient of the latter.
"""
function shape_hessian_gradient_and_value(ip::Interpolation, ξ::Vec, i::Int)
return hessian(x -> shape_value(ip, x, i), ξ, :all)
Expand Down Expand Up @@ -380,15 +381,15 @@
"""
boundarydof_indices(::Type{<:BoundaryIndex})

boundarydof_indices(::Type{FaceIndex}) = Ferrite.facedof_indices
boundarydof_indices(::Type{EdgeIndex}) = Ferrite.edgedof_indices
boundarydof_indices(::Type{VertexIndex}) = Ferrite.vertexdof_indices
boundarydof_indices(::Type{FaceIndex}) = facedof_indices
boundarydof_indices(::Type{EdgeIndex}) = edgedof_indices
boundarydof_indices(::Type{VertexIndex}) = vertexdof_indices

Check warning on line 386 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L384-L386 

Added lines #L384 - L386 were not covered by tests

facetdof_indices(ip::InterpolationByDim{3}) = Ferrite.facedof_indices(ip)
facetdof_indices(ip::InterpolationByDim{2}) = Ferrite.edgedof_indices(ip)
facetdof_indices(ip::InterpolationByDim{1}) = Ferrite.vertexdof_indices(ip)
facetdof_interior_indices(ip::InterpolationByDim{3}) = Ferrite.facedof_interior_indices(ip)
facetdof_interior_indices(ip::InterpolationByDim{2}) = Ferrite.edgedof_interior_indices(ip)
facetdof_indices(ip::InterpolationByDim{3}) = facedof_indices(ip)
facetdof_indices(ip::InterpolationByDim{2}) = edgedof_indices(ip)
facetdof_indices(ip::InterpolationByDim{1}) = vertexdof_indices(ip)
facetdof_interior_indices(ip::InterpolationByDim{3}) = facedof_interior_indices(ip)
facetdof_interior_indices(ip::InterpolationByDim{2}) = edgedof_interior_indices(ip)
facetdof_interior_indices(ip::InterpolationByDim{1}) = ntuple(_ -> (), nvertices(ip))
dirichlet_facetdof_indices(ip::InterpolationByDim{3}) = dirichlet_facedof_indices(ip)
dirichlet_facetdof_indices(ip::InterpolationByDim{2}) = dirichlet_edgedof_indices(ip)
Expand All @@ -415,9 +416,9 @@
"""
dirichlet_boundarydof_indices(::Type{<:BoundaryIndex})

dirichlet_boundarydof_indices(::Type{FaceIndex}) = Ferrite.dirichlet_facedof_indices
dirichlet_boundarydof_indices(::Type{EdgeIndex}) = Ferrite.dirichlet_edgedof_indices
dirichlet_boundarydof_indices(::Type{VertexIndex}) = Ferrite.dirichlet_vertexdof_indices
dirichlet_boundarydof_indices(::Type{FaceIndex}) = dirichlet_facedof_indices

Check warning on line 419 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L419 

Added line #L419 was not covered by tests
dirichlet_boundarydof_indices(::Type{EdgeIndex}) = dirichlet_edgedof_indices
dirichlet_boundarydof_indices(::Type{VertexIndex}) = dirichlet_vertexdof_indices

Check warning on line 421 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L421 

Added line #L421 was not covered by tests
dirichlet_boundarydof_indices(::Type{FacetIndex}) = dirichlet_facetdof_indices

#########################
Expand Down Expand Up @@ -1535,6 +1536,95 @@

is_discontinuous(::Type{<:VectorizedInterpolation{<:Any, <:Any, <:Any, ip}}) where {ip} = is_discontinuous(ip)


# Foundation block for https://github.com/Ferrite-FEM/Ferrite.jl/issues/398 - Remove this below after the issue is resolved.
##################################################
# MatrixizedInterpolation{<:ScalarInterpolation} #
##################################################
struct MatrixizedInterpolation{vdim1, vdim2, refshape, order, SI <: ScalarInterpolation{refshape, order}} <: MatrixInterpolation{vdim1, vdim2, refshape,order}
ip::SI
function MatrixizedInterpolation{vdim1, vdim2}(ip::SI) where {vdim1, vdim2, refshape, order, SI <: ScalarInterpolation{refshape, order}}
return new{vdim1, vdim2, refshape, order, SI}(ip)
end
end

n_components(::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2} = vdim1*vdim2
adjust_dofs_during_distribution(ipm::MatrixizedInterpolation) = adjust_dofs_during_distribution(ipm.ip)

Check warning on line 1552 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1551-L1552 

Added lines #L1551 - L1552 were not covered by tests

# Vectorize to reference dimension by default
function MatrixizedInterpolation(ip::ScalarInterpolation{shape}) where {refdim, shape <: AbstractRefShape{refdim}}
return MatrixizedInterpolation{refdim,refdim}(ip)
end

Base.:(^)(ipv::VectorizedInterpolation{vdim1}, vdim2::Int) where {vdim1} = MatrixizedInterpolation{vdim1, vdim2}(ipv)
function Base.literal_pow(::typeof(^), ipv::VectorizedInterpolation{vdim1}, ::Val{vdim2}) where {vdim1,vdim2}
return MatrixizedInterpolation{vdim1, vdim2}(ipv.ip)

Check warning on line 1561 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1559-L1561 

Added lines #L1559 - L1561 were not covered by tests
end

function Base.show(io::IO, mime::MIME"text/plain", ipm::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2}
show(io, mime, ipm.ip)
print(io, "^", vdim1 , "×", vdim2)

Check warning on line 1566 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1564-L1566 

Added lines #L1564 - L1566 were not covered by tests
end

# Helper to get number of copies for DoF distribution
get_n_copies(::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2} = vdim1*vdim2

Check warning on line 1570 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1570 

Added line #L1570 was not covered by tests

function getnbasefunctions(ipm::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2}
return vdim1 * vdim2 * getnbasefunctions(ipm.ip)
end
function shape_value(ipm::MatrixizedInterpolation{vdim, vdim, shape}, ξ::Tensors.Vec{refdim, T}, I::Int) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T}
# First flatten to vector
i0, c0 = divrem(I - 1, vdim^2)
i = i0 + 1
v = shape_value(ipm.ip, ξ, i)

# Then compute matrix index
ci0, cj0 = divrem(c0, vdim)
ci = ci0 + 1
cj = cj0 + 1
return Tensor{2, vdim, T}((k, l) -> k == cj && l == ci ? v : zero(v))
end

function shape_value(ipm::MatrixizedInterpolation{vdim1, vdim2, shape}, ξ::Tensors.Vec{refdim, T}, I::Int) where {vdim1, vdim2, refdim, shape <: AbstractRefShape{refdim}, T}

Check warning on line 1588 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1588 

Added line #L1588 was not covered by tests
# First flatten to vector
i0, c0 = divrem(I - 1, vdim1*vdim2)
i = i0 + 1
v = shape_value(ipm.ip, ξ, i)

Check warning on line 1592 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1590-L1592 

Added lines #L1590 - L1592 were not covered by tests

# Then compute matrix index
ci0, cj0 = divrem(c0, vdim1)
ci = ci0 + 1
cj = cj0 + 1

Check warning on line 1597 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1595-L1597 

Added lines #L1595 - L1597 were not covered by tests

return SMatrix{vdim1, vdim2, T}(ntuple(i_ -> i_ == (cj0*vdim2 + ci) ? v : zero(v), vdim1*vdim2))

Check warning on line 1599 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1599 

Added line #L1599 was not covered by tests
end

# vdim1 == vdim2 == refdim
function shape_gradient_and_value(ipm::MatrixizedInterpolation{dim, dim, shape}, ξ::Vec{dim}, I::Int) where {dim, shape <: AbstractRefShape{dim}}
return invoke(shape_gradient_and_value, Tuple{Interpolation, Vec, Int}, ipm, ξ, I)
end

# vdim1 != vdim2 != refdim
function shape_gradient_and_value(ipm::MatrixizedInterpolation{vdim1, vdim2, shape}, ξ::V, I::Int) where {vdim1, vdim2, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}}
tosmat(v::Tensor{2}) = SMatrix{vdim1,vdim2}((v...,))
tosmat(v::SMatrix) = v
tosvec(v::Vec) = SVector((v...,))
tovec(sv::SVector) = Vec((sv...))
val = shape_value(ipm, ξ, I)
grad = ForwardDiff.jacobian(sv -> tosmat(shape_value(ipm, tovec(sv), I)), tosvec(ξ))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this will do the right thing unfortunately.
But for these matrix interpolations, it probably is better to calculate the gradients wrt. the scalar base-interpolation, and then just insert the values in the right places (and zeros elsewhere).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have any suggestion how we can cover this case in testing (in addition to the vectorial case)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think @lijas solution here was quite nice, just having the option of calling that function explicitly as setup here.

Maybe combined with calling something like

for Is in 1:getnbasefunctions(ipm.ip)
    g, v = shape_gradient_and_value(ipm.ip, ξ, Is)
    Im = (Is - 1)*getnbasefunctions(ipm.ip)
    for i in 1:vdim1, j in 1:vdim2
        Im += 1
        G, V = shape_gradient_and_value(ipm.ip, ξ, Im)
        @test G[i,j,:] \approx g
        @test V[i,j] \approx v
    end
end

return grad, val

Check warning on line 1615 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1608-L1615 

Added lines #L1608 - L1615 were not covered by tests
end

reference_coordinates(ipm::MatrixizedInterpolation) = reference_coordinates(ipm.ip)

is_discontinuous(::Type{<:MatrixizedInterpolation{<:Any, <:Any, <:Any, <:Any, ip}}) where {ip} = is_discontinuous(ip)

Check warning on line 1620 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1620 

Added line #L1620 was not covered by tests

get_gradient_interpolation(::VectorizedInterpolation{vdim, shape, order, <:Lagrange{shape, order}}) where {sdim,vdim,shape<:AbstractRefShape{sdim},order} = MatrixizedInterpolation{vdim, sdim}(DiscontinuousLagrange{shape, order-1}())

Check warning on line 1622 in src/interpolations.jl

View check run for this annotation

Codecov / codecov/patch

src/interpolations.jl#L1622 

Added line #L1622 was not covered by tests
get_gradient_interpolation_type(::Type{VectorizedInterpolation{vdim, shape, order, <:Lagrange{shape, order}}}) where {sdim,vdim,shape<:AbstractRefShape{sdim},order} = MatrixizedInterpolation{vdim, sdim, shape, order-1, DiscontinuousLagrange{shape, order-1}}

InterpolationInfo(ipm::MatrixizedInterpolation) = InterpolationInfo(ipm.ip, get_n_copies(ipm))


"""
mapping_type(ip::Interpolation)

Expand All @@ -1547,3 +1637,4 @@

mapping_type(::ScalarInterpolation) = IdentityMapping()
mapping_type(::VectorizedInterpolation) = IdentityMapping()
mapping_type(::MatrixizedInterpolation) = IdentityMapping()
2 changes: 2 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@ using LinearAlgebra
using SparseArrays
using OrderedCollections

import Ferrite: MatrixizedInterpolation

const HAS_EXTENSIONS = isdefined(Base, :get_extension)

# https://github.com/JuliaLang/julia/pull/47749
Expand Down
14 changes: 11 additions & 3 deletions test/test_cellvalues.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
(Lagrange{RefPyramid, 2}(), QuadratureRule{RefPyramid}(2)),
)

for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol))
for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol), MatrixizedInterpolation(scalar_interpol))
geom_interpol = scalar_interpol # Tests below assume this
n_basefunc_base = getnbasefunctions(scalar_interpol)
cv = @inferred CellValues(quad_rule, func_interpol, geom_interpol)
Expand All @@ -26,22 +26,26 @@

@test getnbasefunctions(cv) == n_basefuncs

x, n = valid_coordinates_and_normals(func_interpol)
x, n = valid_coordinates_and_normals(scalar_interpol)
reinit!(cv, x)
@test_call reinit!(cv, x)

# We test this by applying a given deformation gradient on all the nodes.
# Since this is a linear deformation we should get back the exact values
# from the interpolation.
u = zeros(Vec{rdim, Float64}, n_basefunc_base)
u2 = zeros(Tensor{2, rdim, Float64}, n_basefunc_base)
u_scal = zeros(n_basefunc_base)
W = rand(Tensor{3, rdim})
H = rand(Tensor{2, rdim})
V = rand(Tensor{1, rdim})
for i in 1:n_basefunc_base
u2[i] = W ⋅ x[i]
u[i] = H ⋅ x[i]
u_scal[i] = V ⋅ x[i]
end
u_vector = reinterpret(Float64, u)
u_matrix = reinterpret(Float64, u2)

for i in 1:getnquadpoints(cv)
if func_interpol isa Ferrite.ScalarInterpolation
Expand All @@ -53,7 +57,7 @@
rdim == 3 && @test function_curl(cv, i, u) ≈ Ferrite.curl_from_gradient(H)
function_value(cv, i, u)
function_value(cv, i, u_scal)
else# func_interpol isa Ferrite.VectorInterpolation
elseif func_interpol isa Ferrite.VectorInterpolation
@test function_gradient(cv, i, u_vector) ≈ H
@test (@test_deprecated function_gradient(cv, i, u)) ≈ H
@test function_symmetric_gradient(cv, i, u_vector) ≈ 0.5(H + H')
Expand All @@ -65,6 +69,10 @@
@test (@test_deprecated function_curl(cv, i, u)) ≈ Ferrite.curl_from_gradient(H)
end
@test function_value(cv, i, u_vector) ≈ (@test_deprecated function_value(cv, i, u))
else func_interpol isa Ferrite.MatrixInterpolation
@test function_gradient(cv, i, u_matrix) ≈ W
@test_throws MethodError function_curl(cv, i, u_matrix)
@test_throws MethodError function_divergence(cv, i, u_matrix)
end
end

Expand Down
2 changes: 1 addition & 1 deletion test/test_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ end
# Volume of cells (with planar edges) #
#######################################

calculate_volume(ip::VectorizedInterpolation, x) = calculate_volume(ip.ip, x)
calculate_volume(ip::Union{VectorizedInterpolation, MatrixizedInterpolation}, x) = calculate_volume(ip.ip, x)

function calculate_volume(::Lagrange{RefLine, 1}, x::Vector{Vec{dim, T}}) where {T, dim}
vol = norm(x[2] - x[1])
Expand Down
Loading