diff --git a/src/voigt.jl b/src/voigt.jl index c868d46a..04db296a 100644 --- a/src/voigt.jl +++ b/src/voigt.jl @@ -1,9 +1,13 @@ const DEFAULT_VOIGT_ORDER = ([1], [1 3; 4 2], [1 6 5; 9 2 4; 8 7 3]) """ - tovoigt(A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...) + tovoigt([type::Type{<:AbstractArray}, ]A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...) Converts a tensor to "Voigt"-format. +Optional argument: +- `type`: determines the returned Array type. Possible types are `Array` and `SArray` (see + [`StaticArrays`](https://juliaarrays.github.io/StaticArrays.jl/stable/)). + Keyword arguments: - `offdiagscale`: determines the scaling factor for the offdiagonal elements. This argument is only applicable for `SymmetricTensor`s. `tomandel` can also @@ -43,24 +47,39 @@ julia> tovoigt(Tensor{4,2}(1:16)) 4 16 12 8 3 15 11 7 2 14 10 6 + +julia> tovoigt(SMatrix, Tensor{4,2}(1:16)) +4×4 SMatrix{4, 4, Int64, 16} with indices SOneTo(4)×SOneTo(4): + 1 13 9 5 + 4 16 12 8 + 3 15 11 7 + 2 14 10 6 ``` + """ function tovoigt end -@inline function tovoigt(A::Tensor{2, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M} - @inbounds tovoigt!(Vector{T}(undef, M), A; order=order) +# default to regular Array +@inline function tovoigt(A::AbstractTensor; kwargs...) + return tovoigt(Array, A; kwargs...) end -@inline function tovoigt(A::Tensor{4, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M} - @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; order=order) +@inline tovoigt(::Type{Array}, A::SecondOrderTensor; kwargs...) = tovoigt(Vector, A; kwargs...) +@inline tovoigt(::Type{Array}, A::FourthOrderTensor; kwargs...) = tovoigt(Matrix, A; kwargs...) + +@inline function tovoigt(::Type{<:Vector}, A::Tensor{2, dim, T, M}; order=nothing) where {dim, T, M} + @inbounds _tovoigt!(Vector{T}(undef, M), A, order) end -@inline function tovoigt(A::SymmetricTensor{2, dim, T, M}; offdiagscale=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M} - @inbounds tovoigt!(Vector{T}(undef, M), A; offdiagscale=offdiagscale, order=order) +@inline function tovoigt(::Type{<:Matrix}, A::Tensor{4, dim, T, M}; order=nothing) where {dim, T, M} + @inbounds _tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A, order) end -@inline function tovoigt(A::SymmetricTensor{4, dim, T, M}; offdiagscale=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M} - @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; offdiagscale=offdiagscale, order=order) +@inline function tovoigt(::Type{<:Vector}, A::SymmetricTensor{2, dim, T, M}; offdiagscale=one(T), order=nothing) where {dim, T, M} + @inbounds _tovoigt!(Vector{T}(undef, M), A, order; offdiagscale=offdiagscale) +end +@inline function tovoigt(::Type{<:Matrix}, A::SymmetricTensor{4, dim, T, M}; offdiagscale=one(T), order=nothing) where {dim, T, M} + @inbounds _tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A, order; offdiagscale=offdiagscale) end """ - tovoigt!(v::Array, A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...) + tovoigt!(v::AbstractArray, A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...) Converts a tensor to "Voigt"-format using the following index order: `[11, 22, 33, 23, 13, 12, 32, 31, 21]`. @@ -105,25 +124,60 @@ julia> tovoigt!(x, T; offset=1) ``` """ function tovoigt! end -Base.@propagate_inbounds function tovoigt!(v::AbstractVector, A::Tensor{2, dim}; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim} +Base.@propagate_inbounds function tovoigt!(v::AbstractVector, A::Tensor{2}; offset::Int=0, order=nothing) + _tovoigt!(v, A, order; offset=offset) +end +Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix, A::Tensor{4}; offset_i::Int=0, offset_j::Int=0, order=nothing) + _tovoigt!(v, A, order; offset_i=offset_i, offset_j=offset_j) +end +Base.@propagate_inbounds function tovoigt!(v::AbstractVector{T}, A::SymmetricTensor{2}; offdiagscale=one(T), offset::Int=0, order=nothing) where T + _tovoigt!(v, A, order; offdiagscale=offdiagscale, offset=offset) +end +Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix{T}, A::SymmetricTensor{4}; offdiagscale=one(T), offset_i::Int=0, offset_j::Int=0, order=nothing) where T + _tovoigt!(v, A, order; offdiagscale=offdiagscale, offset_i=offset_i, offset_j=offset_j) +end + +# default voigt order (faster than custom voigt order) +Base.@propagate_inbounds function _tovoigt!(v::AbstractVecOrMat{T}, A::SecondOrderTensor{dim,T}, ::Nothing; offset=0, offdiagscale=one(T)) where {dim,T} + tuple_data, = _to_voigt_tuple(A, offdiagscale) + for i in eachindex(tuple_data) + v[offset+i] = tuple_data[i] + end + return v +end +Base.@propagate_inbounds function _tovoigt!(v::AbstractVecOrMat{T}, A::SymmetricTensor{4,dim,T}, ::Nothing; offdiagscale=one(T), offset_i=0, offset_j=0) where {dim,T} + tuple_data, N = _to_voigt_tuple(A, offdiagscale) + cartesian = CartesianIndices(((offset_i+1):(offset_i+N), (offset_j+1):(offset_j+N))) + for i in eachindex(tuple_data) + v[cartesian[i]] = tuple_data[i] + end + return v +end + +Base.@propagate_inbounds function _tovoigt!(v::AbstractVecOrMat{T}, A::Tensor{4,dim,T}, ::Nothing; kwargs...) where {dim,T} + return _tovoigt!(v, A, DEFAULT_VOIGT_ORDER[dim]; kwargs...) +end + +# custom voigt order (slower than default voigt order) +Base.@propagate_inbounds function _tovoigt!(v::AbstractVector, A::Tensor{2, dim}, order::AbstractVecOrMat; offset::Int=0) where {dim} for j in 1:dim, i in 1:dim v[offset + order[i, j]] = A[i, j] end return v end -Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix, A::Tensor{4, dim}; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim} +Base.@propagate_inbounds function _tovoigt!(v::AbstractMatrix, A::Tensor{4, dim}, order::AbstractVecOrMat; offset_i::Int=0, offset_j::Int=0) where {dim} for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim v[offset_i + order[i, j], offset_j + order[k, l]] = A[i, j, k, l] end return v end -Base.@propagate_inbounds function tovoigt!(v::AbstractVector{T}, A::SymmetricTensor{2, dim}; offdiagscale=one(T), offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim} +Base.@propagate_inbounds function _tovoigt!(v::AbstractVector{T}, A::SymmetricTensor{2, dim}, order::AbstractVecOrMat; offdiagscale=one(T), offset::Int=0) where {T, dim} for j in 1:dim, i in 1:j v[offset + order[i, j]] = i == j ? A[i, j] : A[i, j] * offdiagscale end return v end -Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix{T}, A::SymmetricTensor{4, dim}; offdiagscale=one(T), offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim} +Base.@propagate_inbounds function _tovoigt!(v::AbstractMatrix{T}, A::SymmetricTensor{4, dim}, order::AbstractVecOrMat; offdiagscale=one(T), offset_i::Int=0, offset_j::Int=0) where {T, dim} for l in 1:dim, k in 1:l, j in 1:dim, i in 1:j v[offset_i + order[i, j], offset_j + order[k, l]] = (i == j && k == l) ? A[i, j, k, l] : @@ -138,8 +192,9 @@ end Convert the tensor `A` to voigt-form using the Mandel convention, see [`tovoigt`](@ref). """ -@inline tomandel(A::SymmetricTensor{o, dim, T}; kwargs...) where{o,dim,T} = tovoigt(A; offdiagscale=√(2one(T)), kwargs...) -@inline tomandel(A::Tensor; kwargs...) = tovoigt(A; kwargs...) +@inline tomandel(A::AbstractTensor; kwargs...) = tomandel(Array, A; kwargs...) +@inline tomandel(::Type{AT}, A::SymmetricTensor{o, dim, T}; kwargs...) where{AT,o,dim,T} = tovoigt(AT, A; offdiagscale=√(2one(T)), kwargs...) +@inline tomandel(::Type{AT}, A::Tensor; kwargs...) where AT = tovoigt(AT, A; kwargs...) """ tomandel!(v, A; kwargs...) @@ -150,7 +205,7 @@ Base.@propagate_inbounds tomandel!(v::AbstractVecOrMat{T}, A::SymmetricTensor; k Base.@propagate_inbounds tomandel!(v::AbstractVecOrMat, A::Tensor; kwargs...) = tovoigt!(v, A; kwargs...) """ - fromvoigt(S::Type{<:AbstractTensor}, A::Array{T}; kwargs...) + fromvoigt(S::Type{<:AbstractTensor}, A::AbstractArray{T}; kwargs...) Converts an array `A` stored in Voigt format to a Tensor of type `S`. @@ -206,3 +261,109 @@ Convert the Array `v` in voigt-format, following the Mandel convention, to a ten """ Base.@propagate_inbounds frommandel(TT::Type{<: SymmetricTensor}, v::AbstractVecOrMat{T}; kwargs...) where{T} = fromvoigt(TT, v; offdiagscale=√(2one(T)), kwargs...) Base.@propagate_inbounds frommandel(TT::Type{<: Tensor}, v::AbstractVecOrMat; kwargs...) = fromvoigt(TT, v; kwargs...) + +############################################################## +# Reorder tensor data to tuple in default voigt format order # +############################################################## +function __to_voigt_tuple(A::Type{TT}, s=one(T)) where {TT<:SecondOrderTensor{dim,T}} where {dim,T} + # Define internals for generation + idx_fun(i, j) = compute_index(get_base(A), i, j) + maxind(j) = TT<:SymmetricTensor ? j : dim + N = n_components(get_base(A)) + + exps = Expr(:tuple) + append!(exps.args, [nothing for _ in 1:N]) # "Preallocate" to allow indexing directly + + for j in 1:dim, i in 1:maxind(j) + voigt_ind = DEFAULT_VOIGT_ORDER[dim][i,j] + if i==j + exps.args[voigt_ind] = :(get_data(A)[$(idx_fun(i, j))]) + else + exps.args[voigt_ind] = :(s*get_data(A)[$(idx_fun(i, j))]) + end + end + return exps, N +end + +function __to_voigt_tuple(A::Type{TT}, s=one(T)) where {TT<:FourthOrderTensor{dim,T}} where {dim,T} + # Define internals for generation + idx_fun(i, j, k, l) = compute_index(get_base(A), i, j, k, l) # why no change needed above? + maxind(j) = TT<:SymmetricTensor ? j : dim + voigt_lin_index(vi, vj) = (vj-1)*N + vi + N = Int(sqrt(n_components(get_base(A)))) + + exps = Expr(:tuple) + append!(exps.args, [nothing for _ in 1:N^2]) # "Preallocate" to allow indexing directly + + for l in 1:dim, k in 1:maxind(l), j in 1:dim, i in 1:maxind(j) + voigt_lin_ind = voigt_lin_index(DEFAULT_VOIGT_ORDER[dim][i,j], DEFAULT_VOIGT_ORDER[dim][k,l]) + if i==j && k==l + exps.args[voigt_lin_ind] = :(get_data(A)[$(idx_fun(i, j, k, l))]) + elseif i!=j && k!=l + exps.args[voigt_lin_ind] = :(s*s*get_data(A)[$(idx_fun(i, j, k, l))]) + else + exps.args[voigt_lin_ind] = :(s*get_data(A)[$(idx_fun(i, j, k, l))]) + end + end + return exps, N +end + +@generated function _to_voigt_tuple(A::AbstractTensor{order, dim, T}, s=one(T)) where {order,dim,T} + exps, N = __to_voigt_tuple(A, s) + quote + $(Expr(:meta, :inline)) + @inbounds return $exps, $N + end +end + +######################## +# StaticArrays support # +######################## +@inline tovoigt(::Type{SArray}, A::SecondOrderTensor; kwargs...) = tovoigt(SVector, A; kwargs...) +@inline tovoigt(::Type{SArray}, A::FourthOrderTensor; kwargs...) = tovoigt(SMatrix, A; kwargs...) +@inline function tovoigt(::Type{<:SVector}, A::Tensor{2, dim, T, M}; order=nothing) where {dim, T, M} + @inbounds _to_static_voigt(A, order) +end +@inline function tovoigt(::Type{<:SMatrix}, A::Tensor{4, dim, T, M}; order=nothing) where {dim, T, M} + @inbounds _to_static_voigt(A, order) +end +@inline function tovoigt(::Type{<:SVector}, A::SymmetricTensor{2, dim, T, M}; offdiagscale=one(T), order=nothing) where {dim, T, M} + @inbounds _to_static_voigt(A, order; offdiagscale=offdiagscale) +end +@inline function tovoigt(::Type{<:SMatrix}, A::SymmetricTensor{4, dim, T}; offdiagscale=one(T), order=nothing) where {dim, T} + @inbounds _to_static_voigt(A, order; offdiagscale=offdiagscale) +end + +# default voigt order +Base.@propagate_inbounds function _to_static_voigt(A::TT, ::Nothing; offdiagscale=one(T)) where {TT<:SecondOrderTensor{dim,T}} where {dim,T} + tuple_data, N = _to_voigt_tuple(A, offdiagscale) + return SVector{N, T}(tuple_data) +end +Base.@propagate_inbounds function _to_static_voigt(A::TT, ::Nothing; offdiagscale=one(T)) where {TT<:FourthOrderTensor{dim,T}} where {dim,T} + tuple_data, N = _to_voigt_tuple(A, offdiagscale) + return SMatrix{N, N, T}(tuple_data) +end + +# custom voigt order +@inline function _to_static_voigt(A::Tensor{2, dim, T, M}, order::AbstractVecOrMat) where {dim, T, M} + v = MVector{M, T}(undef) + @inbounds _tovoigt!(v, A, order) + return SVector(v) +end +@inline function _to_static_voigt(A::Tensor{4, dim, T, M}, order::AbstractVecOrMat) where {dim, T, M} + m = MMatrix{Int(√M), Int(√M), T}(undef) + @inbounds _tovoigt!(m, A, order) + return SMatrix(m) +end +@inline function _to_static_voigt(A::SymmetricTensor{2, dim, T, M}, order::AbstractVecOrMat; offdiagscale=one(T)) where {dim, T, M} + v = MVector{M, T}(undef) + @inbounds _tovoigt!(v, A, order; offdiagscale=offdiagscale) + return SVector(v) +end +@inline function _to_static_voigt(A::SymmetricTensor{4, dim, T, M}, order::AbstractVecOrMat; offdiagscale=one(T)) where {dim, T, M} + m = MMatrix{Int(√M), Int(√M), T}(undef) + @inbounds _tovoigt!(m, A, order; offdiagscale=offdiagscale) + return SMatrix(m) +end + + diff --git a/test/runtests.jl b/test/runtests.jl index 036373be..99487afe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using TimerOutputs using LinearAlgebra using Random using Statistics: mean +using StaticArrays macro testsection(str, block) return quote diff --git a/test/test_ops.jl b/test/test_ops.jl index 957626a1..e12e9f6f 100644 --- a/test/test_ops.jl +++ b/test/test_ops.jl @@ -321,6 +321,20 @@ end @test isa(fromvoigt(SymmetricTensor{2,dim}, rand(num_components) .> 0.5), SymmetricTensor{2,dim,Bool}) @test isa(fromvoigt(SymmetricTensor{4,dim}, rand(num_components,num_components) .> 0.5), SymmetricTensor{4,dim,Bool}) end + + # Static voigt/mandel that use default order + s2 = [rand(SymmetricTensor{2,dim}) for dim in 1:3] + s4 = [rand(SymmetricTensor{4,dim}) for dim in 1:3] + t2 = [rand(Tensor{2,dim}) for dim in 1:3] + t4 = [rand(Tensor{4,dim}) for dim in 1:3] + for a in (s2, s4, t2, t4) + for b in a + @test tovoigt(SArray, b) == tovoigt(b) + @test tomandel(SArray, b) ≈ tomandel(b) + @test fromvoigt(Tensors.get_base(typeof(b)), tovoigt(SArray, b)) ≈ b + @test frommandel(Tensors.get_base(typeof(b)), tomandel(SArray, b)) ≈ b + end + end end end # of testsection end # of testsection