Convert to voigt and mandel as SVector/SMatrix

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
+
+

test/runtests.jl

@@ -4,6 +4,7 @@ using TimerOutputs
 using LinearAlgebra
 using Random
 using Statistics: mean
+using StaticArrays
 
 macro testsection(str, block)
     return quote

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