From fdf0b284e0f87c5e5c101cd31fdc61cadf833ab6 Mon Sep 17 00:00:00 2001 From: Antoine Marteau Date: Mon, 5 Aug 2024 11:56:29 +1000 Subject: [PATCH] Added type for traceless symmetric tensors - Added SymTracelessTensorValue (alias QTensorValue) It is a symmetric tensor Multivalue with null trace, the last diagonal value should not be provided in the constructor, the latter computes it and stores it in the data field - Added AbstractSymTensorValue{D} <: Multivalue{Tuple{D,D}} as abstract supertype for SymTensorValue and SymTracelessTensorValue --- src/TensorValues/Indexing.jl | 6 +- src/TensorValues/Operations.jl | 252 ++++++++++-------- src/TensorValues/SymTensorValueTypes.jl | 8 +- .../SymTracelessTensorValueTypes.jl | 176 ++++++++++++ src/TensorValues/TensorValues.jl | 5 + test/TensorValuesTests/IndexingTests.jl | 23 +- test/TensorValuesTests/OperationsTests.jl | 162 ++++++++++- test/TensorValuesTests/ReinterpretTests.jl | 6 + test/TensorValuesTests/TypesTests.jl | 84 +++++- 9 files changed, 599 insertions(+), 123 deletions(-) create mode 100644 src/TensorValues/SymTracelessTensorValueTypes.jl diff --git a/src/TensorValues/Indexing.jl b/src/TensorValues/Indexing.jl index 8f9f996d8..ee47e7f9d 100644 --- a/src/TensorValues/Indexing.jl +++ b/src/TensorValues/Indexing.jl @@ -14,7 +14,7 @@ function getindex(arg::TensorValue{D},i::Integer,j::Integer) where D arg.data[index] end -function getindex(arg::SymTensorValue{D},i::Integer,j::Integer) where D +function getindex(arg::AbstractSymTensorValue{D},i::Integer,j::Integer) where D index = _2d_sym_tensor_linear_index(D,i,j) arg.data[index] end @@ -31,7 +31,7 @@ end getindex(arg::VectorValue, ci::CartesianIndex{1}) = getindex(arg,ci[1]) getindex(arg::TensorValue,ci::CartesianIndex{2}) = getindex(arg,ci[1],ci[2]) -getindex(arg::SymTensorValue,ci::CartesianIndex{2}) = getindex(arg,ci[1],ci[2]) +getindex(arg::AbstractSymTensorValue,ci::CartesianIndex{2}) = getindex(arg,ci[1],ci[2]) getindex(arg::ThirdOrderTensorValue,ci::CartesianIndex{3}) = getindex(arg,ci[1],ci[2],ci[3]) getindex(arg::SymFourthOrderTensorValue,ci::CartesianIndex{4}) = getindex(arg,ci[1],ci[2],ci[3],ci[4]) @@ -44,7 +44,7 @@ getindex(arg::ThirdOrderTensorValue, i::Integer) = arg.data[i] data_index(::Type{<:VectorValue},i) = i data_index(::Type{<:TensorValue{D}},i,j) where D = _2d_tensor_linear_index(D,i,j) -data_index(::Type{<:SymTensorValue{D}},i,j) where D = _2d_sym_tensor_linear_index(D,i,j) +data_index(::Type{<:AbstractSymTensorValue{D}},i,j) where D = _2d_sym_tensor_linear_index(D,i,j) data_index(::Type{<:ThirdOrderTensorValue{D1,D2}},i,j,k) where {D1,D2} = _3d_tensor_linear_index(D1,D2,i,j,k) data_index(::Type{<:SymFourthOrderTensorValue{D}},i,j,k,l) where D = _4d_sym_tensor_linear_index(D,i,j,k,l) diff --git a/src/TensorValues/Operations.jl b/src/TensorValues/Operations.jl index e81d90b55..4f831de65 100644 --- a/src/TensorValues/Operations.jl +++ b/src/TensorValues/Operations.jl @@ -45,6 +45,10 @@ for op in (:+,:-) T(r) end + function ($op)(a::MultiValue,b::MultiValue) + @notimplemented "Not implemented or undefined operation \"$($op)\" on MultiValues of these shapes" + end + function ($op)(a::MultiValue{S},b::MultiValue{S}) where S r = map(($op), a.data, b.data) T = _eltype($op,r,a,b) @@ -52,17 +56,52 @@ for op in (:+,:-) M(r) end - function ($op)(a::TensorValue,b::SymTensorValue) + function ($op)(a::TensorValue{D,D},b::SymTensorValue{D}) where D + map(($op), a, TensorValue(get_array(b))) + end + + function ($op)(a::SymTensorValue{D},b::TensorValue{D,D}) where D + map(($op), TensorValue(get_array(a)), b) + end + + function ($op)(a::TensorValue{D,D},b::SymTracelessTensorValue{D}) where D map(($op), a, TensorValue(get_array(b))) end - function ($op)(a::SymTensorValue,b::TensorValue) + function ($op)(a::SymTracelessTensorValue{D},b::TensorValue{D,D}) where D map(($op), TensorValue(get_array(a)), b) end + function ($op)(a::SymTracelessTensorValue{D},b::SymTensorValue{D}) where D + r = map(($op), a.data, b.data) + T = _eltype($op,r,a,b) + M = change_eltype(b,T) + M(r) + end + + function ($op)(a::SymTensorValue{D},b::SymTracelessTensorValue{D}) where D + r = map(($op), a.data, b.data) + T = _eltype($op,r,a,b) + M = change_eltype(a,T) + M(r) + end + + function ($op)(a::SymTracelessTensorValue) + r = map($op, a.data[1:end-1]) + typeof(a)(r) + end + + function ($op)(a::SymTracelessTensorValue{D},b::SymTracelessTensorValue{D}) where D + r = map(($op), a.data[1:end-1], b.data[1:end-1]) + T = _eltype($op,r,a,b) + M = change_eltype(a,T) + M(r) + end + end end + ############################################################### # Matrix Division ############################################################### @@ -98,28 +137,45 @@ end for op in (:+,:-,:*) @eval begin function ($op)(a::MultiValue,b::Number) - r = _bc($op,a.data,b) - T = _eltype($op,r,a,b) - M = change_eltype(a,T) - M(r) + r = _bc($op,a.data,b) + T = _eltype($op,r,a,b) + M = change_eltype(a,T) + M(r) end function ($op)(a::Number,b::MultiValue) - r = _bc($op,a,b.data) - T = _eltype($op,r,a,b) - M = change_eltype(b,T) - M(r) + r = _bc($op,a,b.data) + T = _eltype($op,r,a,b) + M = change_eltype(b,T) + M(r) end end end +function (*)(a::Number,b::SymTracelessTensorValue) + r = _bc(*,a,b.data[1:end-1]) + T = _eltype(*,r,a,b) + M = change_eltype(b,T) + M(r) +end + +function (*)(a::SymTracelessTensorValue,b::Number) + b*a +end + function (/)(a::MultiValue,b::Number) - r = _bc(/,a.data,b) - T = _eltype(/,r,a,b) - P = change_eltype(a,T) - P(r) + r = _bc(/,a.data,b) + T = _eltype(/,r,a,b) + P = change_eltype(a,T) + P(r) end +const _err = " with number is undefined for traceless tensors" +function -(::SymTracelessTensorValue,::Number) error("Addition" *_err) end +function +(::SymTracelessTensorValue,::Number) error("Subtraction"*_err) end +function -(::Number,::SymTracelessTensorValue) error("Addition" *_err) end +function +(::Number,::SymTracelessTensorValue) error("Subtraction"*_err) end + @inline function _eltype(op,r,a,b) eltype(r) end @@ -146,18 +202,18 @@ dot(a::MultiValue{Tuple{D}}, b::MultiValue{Tuple{D}}) where D = inner(a,b) dot(a::MultiValue,b::MultiValue) = @notimplemented @generated function dot(a::A,b::B) where {A<:MultiValue{Tuple{D1}},B<:MultiValue{Tuple{D1,D2}}} where {D1,D2} - ss = String[] - for j in 1:D2 - s = "" - for i in 1:D1 - ak = data_index(A,i) - bk = data_index(B,i,j) - s *= "a.data[$ak]*b.data[$bk]+" - end - push!(ss,s[1:(end-1)]*", ") + ss = String[] + for j in 1:D2 + s = "" + for i in 1:D1 + ak = data_index(A,i) + bk = data_index(B,i,j) + s *= "a.data[$ak]*b.data[$bk]+" end - str = join(ss) - Meta.parse("VectorValue{$D2}($str)") + push!(ss,s[1:(end-1)]*", ") + end + str = join(ss) + Meta.parse("VectorValue{$D2}($str)") end function dot(a::A,b::B) where {A<:MultiValue{Tuple{0}},B<:MultiValue{Tuple{0,D2}}} where D2 @@ -166,30 +222,30 @@ function dot(a::A,b::B) where {A<:MultiValue{Tuple{0}},B<:MultiValue{Tuple{0,D2} end @generated function dot(a::A,b::B) where {A<:MultiValue{Tuple{D1,D2}},B<:MultiValue{Tuple{D2}}} where {D1,D2} - ss = String[] - for i in 1:D1 - s = "" - for j in 1:D2 - ak = data_index(A,i,j) - bk = data_index(B,j) - s *= "a.data[$ak]*b.data[$bk]+" - end - push!(ss,s[1:(end-1)]*", ") + ss = String[] + for i in 1:D1 + s = "" + for j in 1:D2 + ak = data_index(A,i,j) + bk = data_index(B,j) + s *= "a.data[$ak]*b.data[$bk]+" end - str = join(ss) - Meta.parse("VectorValue{$D1}($str)") + push!(ss,s[1:(end-1)]*", ") + end + str = join(ss) + Meta.parse("VectorValue{$D1}($str)") end @generated function dot(a::MultiValue{Tuple{D1,D3}}, b::MultiValue{Tuple{D3,D2}}) where {D1,D2,D3} - ss = String[] - for j in 1:D2 - for i in 1:D1 - s = join([ "a[$i,$k]*b[$k,$j]+" for k in 1:D3]) - push!(ss,s[1:(end-1)]*", ") - end + ss = String[] + for j in 1:D2 + for i in 1:D1 + s = join([ "a[$i,$k]*b[$k,$j]+" for k in 1:D3]) + push!(ss,s[1:(end-1)]*", ") end - str = join(ss) - Meta.parse("TensorValue{$D1,$D2}(($str))") + end + str = join(ss) + Meta.parse("TensorValue{$D1,$D2}(($str))") end # a_ij = b_ijk*c_k @@ -261,34 +317,35 @@ function inner(a::MultiValue, b::MultiValue) end @generated function inner(a::MultiValue{S}, b::MultiValue{S}) where S - str = join([" a[$i]*b[$i] +" for i in 1:length(a) ]) - Meta.parse(str[1:(end-1)]) + str = join([" a[$i]*b[$i] +" for i in 1:length(a) ]) + Meta.parse(str[1:(end-1)]) end -@generated function inner(a::SymTensorValue{D}, b::SymTensorValue{D}) where D +@generated function inner(a::AbstractSymTensorValue{D}, b::AbstractSymTensorValue{D}) where D str = "" for i in 1:D - for j in 1:D - k = data_index(a,i,j) - str *= " a.data[$k]*b.data[$k] +" + str *= "+ a[$i,$i]*b[$i,$i]" + end + str *= " + 2*(" + for i in 1:D + for j in i+1:D + str *= "+ a[$i,$j]*b[$i,$j]" end end - Meta.parse(str[1:(end-1)]) + str *= ")" + Meta.parse(str) end -@generated function inner(a::SymFourthOrderTensorValue{D}, b::SymTensorValue{D}) where D +@generated function inner(a::SymFourthOrderTensorValue{D}, b::AbstractSymTensorValue{D}) where D str = "" for i in 1:D for j in i:D - s = "" for k in 1:D for l in 1:D - ak = data_index(a,i,j,k,l) - bk = data_index(b,k,l) - s *= " a.data[$ak]*b.data[$bk] +" + str *= "+ a[$i,$j,$k,$l]*b[$k,$l]" end end - str *= s[1:(end-1)]*", " + str *= ", " end end Meta.parse("SymTensorValue{D}($str)") @@ -388,9 +445,9 @@ const ⋅² = double_contraction ############################################################### for op in (:sum,:maximum,:minimum) - @eval begin - $op(a::MultiValue) = $op(a.data) - end + @eval begin + $op(a::MultiValue) = $op(a.data) + end end # Outer product (aka dyadic product) @@ -426,15 +483,13 @@ end Meta.parse("ThirdOrderTensorValue{D,D1,D2}($str)") end -@generated function outer(a::SymTensorValue{D},b::SymTensorValue{D}) where D +@generated function outer(a::AbstractSymTensorValue{D},b::AbstractSymTensorValue{D}) where D str = "" for i in 1:D for j in i:D - ak = data_index(a,i,j) for k in 1:D for l in k:D - bk = data_index(b,k,l) - str *= "a.data[$ak]*b.data[$bk], " + str *= "a[$i,$j]*b[$k,$l], " end end end @@ -481,6 +536,9 @@ function det(a::MultiValue{Tuple{3,3}}) end inv(a::MultiValue{Tuple{D1,D2}}) where {D1,D2} = TensorValue(inv(get_array(a))) +# those still have better perf than the D=2,3 specialization below +inv(a::AbstractSymTensorValue{D}) where D = SymTensorValue(inv(get_array(a))) +inv(a::SymTracelessTensorValue{2}) = SymTracelessTensorValue(inv(get_array(a))) function inv(a::MultiValue{Tuple{1,1}}) r = 1/a[1] @@ -520,7 +578,8 @@ end """ meas(a::MultiValue{Tuple{D}}) where D = sqrt(inner(a,a)) meas(a::MultiValue{Tuple{D,D}}) where D = abs(det(a)) -meas(a::TensorValue{0,D,T}) where {T,D} = one(T) +#meas( ::TensorValue{0,D,T}) where {T,D} = one(T) +#meas( ::MultiValue{Tuple{0,0},T}) where {T} = one(T) function meas(v::MultiValue{Tuple{1,D}}) where D t = VectorValue(v.data) @@ -553,14 +612,20 @@ function conj(a::T) where {T<:MultiValue} T(r) end +function conj(a::SymTracelessTensorValue) + r = map(conj, a.data) + SymTracelessTensorValue(r[1:end-1]) +end + ############################################################### # Trace ############################################################### @generated function tr(v::MultiValue{Tuple{D,D}}) where D - str = join([" v[$i,$i] +" for i in 1:D ]) - Meta.parse(str[1:(end-1)]) + str = join([" v[$i,$i] +" for i in 1:D ]) + Meta.parse(str[1:(end-1)]) end +tr(::SymTracelessTensorValue{D,T}) where {D,T} = zero(T) @generated function tr(v::MultiValue{Tuple{A,A,B}}) where {A,B} lis = LinearIndices((A,A,B)) @@ -596,7 +661,7 @@ transpose(a::MultiValue{Tuple{D,D}}) where D = @notimplemented Meta.parse("TensorValue{D2,D1}($str)") end -@generated function transpose(a::TensorValue{D1,D2}) where {D1,D2} +@generated function transpose(a::TensorValue{D1,D2,T}) where {D1,D2,T} str = "" for i in 1:D1 for j in 1:D2 @@ -604,18 +669,18 @@ end str *= "a.data[$k], " end end - Meta.parse("TensorValue{D2,D1}($str)") + Meta.parse("TensorValue{D2,D1,T}($str)") end @inline function adjoint(a::TensorValue{D1,D2,T}) where {D1,D2,T<:Real} transpose(a) end -adjoint(a::SymTensorValue) = conj(a) +adjoint(a::AbstractSymTensorValue) = conj(a) -@inline adjoint(a::SymTensorValue{D,T} where {D,T<:Real}) = transpose(a) +@inline adjoint(a::AbstractSymTensorValue{D,T} where {D,T<:Real}) = transpose(a) -transpose(a::SymTensorValue) = a +transpose(a::AbstractSymTensorValue) = a ############################################################### # Symmetric part @@ -634,24 +699,14 @@ transpose(a::SymTensorValue) = a Meta.parse("SymTensorValue{D}($str)") end +symmetric_part(v::AbstractSymTensorValue) = v + ############################################################### # diag ############################################################### -function LinearAlgebra.diag(a::TensorValue{1,1}) - VectorValue(a.data[1]) -end - -function LinearAlgebra.diag(a::TensorValue{2,2}) - VectorValue(a.data[1],a.data[4]) -end - -function LinearAlgebra.diag(a::TensorValue{3,3}) - VectorValue(a.data[1],a.data[5],a.data[9]) -end - -function LinearAlgebra.diag(a::TensorValue) - @notimplemented +function LinearAlgebra.diag(a::MultiValue{Tuple{D,D},T}) where {D,T} + VectorValue((a[i,i] for i in 1:D)...) end ############################################################### @@ -667,25 +722,6 @@ function Base.broadcasted(f,a::TensorValue,b::TensorValue) TensorValue(map(f,a.data,b.data)) end -############################################################### -# Define new operations for Gridap types -############################################################### - -#for op in (:symmetric_part,) -# @eval begin -# ($op)(a::GridapType) = operate($op,a) -# end -#end -# -#for op in (:inner,:outer,:double_contraction)#,:(:)) -# @eval begin -# ($op)(a::GridapType,b::GridapType) = operate($op,a,b) -# ($op)(a::GridapType,b::Number) = operate($op,a,b) -# ($op)(a::Number, b::GridapType) = operate($op,a,b) -# ($op)(a::GridapType,b::Function) = operate($op,a,b) -# ($op)(a::Function, b::GridapType) = operate($op,a,b) -# end -#end - - - +function Base.broadcasted(f,a::AbstractSymTensorValue,b::AbstractSymTensorValue) + SymTensorValue(map(f,a.data,b.data)) +end diff --git a/src/TensorValues/SymTensorValueTypes.jl b/src/TensorValues/SymTensorValueTypes.jl index 1f9412aef..3382f23cb 100644 --- a/src/TensorValues/SymTensorValueTypes.jl +++ b/src/TensorValues/SymTensorValueTypes.jl @@ -1,11 +1,15 @@ ############################################################### # SymTensorValue Type ############################################################### +""" +Abstract type representing any symmetric second-order tensor +""" +abstract type AbstractSymTensorValue{D,T,L} <: MultiValue{Tuple{D,D},T,2,L} end """ -Type representing a symmetric second-order tensor +Type representing a symmetric second-order tensor (with D(D-1)/2 independant components) """ -struct SymTensorValue{D,T,L} <: MultiValue{Tuple{D,D},T,2,L} +struct SymTensorValue{D,T,L} <: AbstractSymTensorValue{D,T,L} data::NTuple{L,T} function SymTensorValue{D,T}(data::NTuple{L,T}) where {D,T,L} @check L == D*(D+1)÷2 diff --git a/src/TensorValues/SymTracelessTensorValueTypes.jl b/src/TensorValues/SymTracelessTensorValueTypes.jl new file mode 100644 index 000000000..4bd847e3a --- /dev/null +++ b/src/TensorValues/SymTracelessTensorValueTypes.jl @@ -0,0 +1,176 @@ +############################################################### +# SymTracelessTensorValue Type +############################################################### + +""" +Type representing a traceless symmetric second-order tensor, +used to model the Q tensor in nematic liquid cristals + +The last diagonal value is determined by minus the sum of the other and musn't be provided +""" +struct SymTracelessTensorValue{D,T,L} <: AbstractSymTensorValue{D,T,L} + data::NTuple{L,T} + function SymTracelessTensorValue{D,T}(data::NTuple{L,T}) where {D,T,L} + @check L == D*(D+1)÷2-1 + new{D,T,L+1}( (data..., _minus_trace(data,Val(D))) ) + end + function SymTracelessTensorValue{0,T}(data::NTuple{0,T}) where {T} + new{0,T,0}(data) + end + function SymTracelessTensorValue{1,T}(::NTuple{0,T}) where {T} + new{1,T,1}( (zero(T),) ) + end +end + +@generated function _minus_trace(data::NTuple{L,T},::Val{D}) where {D,T,L} + str = "" + for i in 1:D-1 + k = _2d_sym_tensor_linear_index(D,i,i) + str *= "- data[$k]" + end + Meta.parse("($str)") +end + +const QTensorValue = SymTracelessTensorValue + +############################################################### +# Constructors (SymTracelessTensorValue) +############################################################### + +# Empty SymTracelessTensorValue constructor + +SymTracelessTensorValue() = SymTracelessTensorValue{0,Int}(NTuple{0,Int}()) +SymTracelessTensorValue{0}() = SymTracelessTensorValue{0,Int}(NTuple{0,Int}()) +SymTracelessTensorValue{0,T}() where {T} = SymTracelessTensorValue{0,T}(NTuple{0,T}()) +SymTracelessTensorValue(data::NTuple{0}) = SymTracelessTensorValue{0,Int}(data) +SymTracelessTensorValue{0}(data::NTuple{0}) = SymTracelessTensorValue{0,Int}(data) + +# 1D SymTracelessTensorValue missing constructor + +SymTracelessTensorValue{1}() = SymTracelessTensorValue{1,Int}(NTuple{0,Int}()) +SymTracelessTensorValue{1}(data::NTuple{0}) = SymTracelessTensorValue{1,Int}(data) + +# SymTracelessTensorValue single NTuple argument constructor + +@generated function SymTracelessTensorValue(data::NTuple{L,T}) where {L,T} + msg = "Invalid number of scalar arguments in SymTracelessTensorValue constructor" + V = (sqrt(9+8*L)-1)/2 + @check floor(Int,V) == ceil(Int,V) msg + D = Int(V) + quote + SymTracelessTensorValue{$D,T}(data) + end +end +SymTracelessTensorValue{D}(data::NTuple{L,T}) where {D,L,T} = SymTracelessTensorValue{D,T}(data) +SymTracelessTensorValue{D,T1}(data::NTuple{L,T2}) where {D,L,T1,T2} = SymTracelessTensorValue{D,T1}(NTuple{L,T1}(data)) +SymTracelessTensorValue{D,T1,L}(data::NTuple{Lm,T2}) where {D,L,Lm,T1,T2} = SymTracelessTensorValue{D,T1}(NTuple{Lm,T1}(data)) + +# SymTracelessTensorValue single Tuple argument constructor + +SymTracelessTensorValue(data::Tuple) = SymTracelessTensorValue(promote(data...)) +SymTracelessTensorValue{D}(data::Tuple) where {D} = SymTracelessTensorValue{D}(promote(data...)) +SymTracelessTensorValue{D,T1}(data::Tuple) where {D,T1} = SymTracelessTensorValue{D,T1}(NTuple{length(data),T1}(data)) +SymTracelessTensorValue{D,T1,L}(data::Tuple) where {D,T1,L} = SymTracelessTensorValue{D,T1}(NTuple{L-1,T1}(data)) + +# SymTracelessTensorValue Vararg constructor + +SymTracelessTensorValue(data::Number...) = SymTracelessTensorValue(data) +SymTracelessTensorValue{D}(data::Number...) where {D} = SymTracelessTensorValue{D}(data) +SymTracelessTensorValue{D,T1}(data::Number...) where {D,T1} = SymTracelessTensorValue{D,T1}(data) +SymTracelessTensorValue{D,T1,L}(data::Number...) where {D,T1,L} = SymTracelessTensorValue{D,T1}(data) + +# SymTracelessTensorValue single AbstractMatrix argument constructor + +#From Square Matrices +@generated function _flatten_upper_triangle_traceless(data::AbstractArray,::Val{D}) where D + str = "" + for i in 1:D-1 + for j in i:D + str *= "data[$i,$j], " + end + end + Meta.parse("($str)") +end + +SymTracelessTensorValue(data::AbstractMatrix{T}) where {T} = ((D1,D2)=size(data); SymTracelessTensorValue{D1}(data)) +SymTracelessTensorValue{D}(data::AbstractMatrix{T}) where {D,T} = SymTracelessTensorValue{D,T}(_flatten_upper_triangle_traceless(data,Val{D}())) +SymTracelessTensorValue{D,T1}(data::AbstractMatrix{T2}) where {D,T1,T2} = SymTracelessTensorValue{D,T1}(_flatten_upper_triangle_traceless(data,Val{D}())) +SymTracelessTensorValue{D,T1,L}(data::AbstractMatrix{T2}) where {D,T1,T2,L} = SymTracelessTensorValue{D,T1,L}(_flatten_upper_triangle_traceless(data,Val{D}())) + +############################################################### +# Conversions (SymTracelessTensorValue) +############################################################### + +@generated function _SymTracelessTensorValue_to_array(arg::SymTracelessTensorValue{D,T,L}) where {D,T,L} + str = "" + for j in 1:D + for i in 1:D + str *= "arg[$i,$j], " + end + end + Meta.parse("SMatrix{D,D,T}(($str))") +end + +# Direct conversion +convert(::Type{<:SymTracelessTensorValue{D,T}}, arg::AbstractArray) where {D,T} = SymTracelessTensorValue{D,T}(arg) +convert(::Type{<:SymTracelessTensorValue{D,T}}, arg::Tuple) where {D,T} = SymTracelessTensorValue{D,T}(arg) + +# Inverse conversion +convert(::Type{<:MMatrix{D,D,T}}, arg::SymTracelessTensorValue) where {D,T} = MMatrix{D,D,T}(_SymTracelessTensorValue_to_array(arg)) +convert(::Type{<:SMatrix{D,D,T}}, arg::SymTracelessTensorValue) where {D,T} = _SymTracelessTensorValue_to_array(arg) +convert(::Type{<:NTuple{L,T}}, arg::SymTracelessTensorValue) where {L,T} = NTuple{L,T}(Tuple(arg)) + +# Internal conversion +convert(::Type{<:SymTracelessTensorValue{D,T}}, arg::SymTracelessTensorValue{D}) where {D,T} = SymTracelessTensorValue{D,T}(Tuple(arg)) +convert(::Type{<:SymTracelessTensorValue{D,T}}, arg::SymTracelessTensorValue{D,T}) where {D,T} = arg + +############################################################### +# Other constructors and conversions (SymTracelessTensorValue) +############################################################### + +zero(::Type{<:SymTracelessTensorValue{0,T}}) where {T} = SymTracelessTensorValue{0,T}() +@generated function zero(::Type{<:SymTracelessTensorValue{D,T}}) where {D,T} + L=D*(D+1)÷2-1 + quote + SymTracelessTensorValue{D,T}(tfill(zero(T),Val{$L}())) + end +end + +zero(::Type{<:SymTracelessTensorValue{D,T,L}}) where {D,T,L} = SymTracelessTensorValue{D,T}(tfill(zero(T),Val{L}())) +zero(::SymTracelessTensorValue{D,T,L}) where {D,T,L} = zero(SymTracelessTensorValue{D,T,L}) + +rand(::AbstractRNG, ::Random.SamplerType{<:SymTracelessTensorValue{0,T}}) where {T} = SymTracelessTensorValue{0,T}() +@generated function rand(rng::AbstractRNG, + ::Random.SamplerType{<:SymTracelessTensorValue{D,T}}) where {D,T} + L=D*(D+1)÷2 + quote + rand(rng, SymTracelessTensorValue{D,T,$L}) + end +end +rand(rng::AbstractRNG,::Random.SamplerType{<:SymTracelessTensorValue{D,T,L}}) where {D,T,L} = + SymTracelessTensorValue{D,T}(Tuple(rand(rng, SVector{L-1,T}))) + +Mutable(::Type{<:SymTracelessTensorValue{D,T}}) where {D,T} = MMatrix{D,D,T} +Mutable(::SymTracelessTensorValue{D,T}) where {D,T} = Mutable(SymTracelessTensorValue{D,T}) +mutable(a::SymTracelessTensorValue{D}) where D = MMatrix{D,D}(Tuple(get_array(a))) + +change_eltype(::Type{SymTracelessTensorValue{D,T1,L}},::Type{T2}) where {D,T1,T2,L} = SymTracelessTensorValue{D,T2,L} +change_eltype(::SymTracelessTensorValue{D,T1,L},::Type{T2}) where {D,T1,T2,L} = change_eltype(SymTracelessTensorValue{D,T1,L},T2) + +get_array(arg::SymTracelessTensorValue{D,T,L}) where {D,T,L} = convert(SMatrix{D,D,T}, arg) + +############################################################### +# Introspection (SymTracelessTensorValue) +############################################################### + +eltype(::Type{<:SymTracelessTensorValue{D,T}}) where {D,T} = T +eltype(::SymTracelessTensorValue{D,T}) where {D,T} = eltype(SymTracelessTensorValue{D,T}) + +size(::Type{<:SymTracelessTensorValue{D}}) where {D} = (D,D) +size(::SymTracelessTensorValue{D}) where {D} = size(SymTracelessTensorValue{D}) + +length(::Type{<:SymTracelessTensorValue{D}}) where {D} = D*D +length(::SymTracelessTensorValue{D}) where {D} = length(SymTracelessTensorValue{D}) + +num_components(::Type{<:SymTracelessTensorValue{D}}) where {D} = length(SymTracelessTensorValue{D}) +num_components(::SymTracelessTensorValue{D}) where {D} = num_components(SymTracelessTensorValue{D}) diff --git a/src/TensorValues/TensorValues.jl b/src/TensorValues/TensorValues.jl index 513d046a0..747e21d71 100644 --- a/src/TensorValues/TensorValues.jl +++ b/src/TensorValues/TensorValues.jl @@ -41,7 +41,10 @@ using Random export MultiValue export VectorValue export TensorValue +export AbstractSymTensorValue export SymTensorValue +export SymTracelessTensorValue +export QTensorValue export SymFourthOrderTensorValue export ThirdOrderTensorValue @@ -89,6 +92,8 @@ include("TensorValueTypes.jl") include("SymTensorValueTypes.jl") +include("SymTracelessTensorValueTypes.jl") + include("SymFourthOrderTensorValueTypes.jl") include("ThirdOrderTensorValueTypes.jl") diff --git a/test/TensorValuesTests/IndexingTests.jl b/test/TensorValuesTests/IndexingTests.jl index 5f6951076..378c23054 100644 --- a/test/TensorValuesTests/IndexingTests.jl +++ b/test/TensorValuesTests/IndexingTests.jl @@ -40,24 +40,33 @@ for (k,ti) in enumerate(t) end s = SymTensorValue{2}(11,21,22) +q = SymTracelessTensorValue{2}(11,21) t = TensorValue(convert(SMatrix{2,2,Int},s)) +p = TensorValue(convert(SMatrix{2,2,Int},q)) -@test size(s) == (2,2) -@test length(s) == 4 -@test lastindex(s) == length(s) -@test s[end] == 22 +@test size(s) == size(q) == (2,2) +@test length(s) == length(q) == 4 +@test lastindex(s) == lastindex(q) == length(s) +@test s[end] == 22 +@test q[end] == -11 for (k,i) in enumerate(eachindex(t)) @test s[i] == t[k] end +for (k,i) in enumerate(eachindex(p)) + @test q[i] == p[k] +end -@test s[2,1] == 21 - -@test s[2] == 21 +@test s[2,1] == q[2,1] == 21 +@test s[2] == q[2] == 21 +@test q[1] == -q[4] for (k,si) in enumerate(t) @test si == s[k] end +for (k,qi) in enumerate(p) + @test qi == q[k] +end v = @SMatrix zeros(2,3) w = TensorValue(v) diff --git a/test/TensorValuesTests/OperationsTests.jl b/test/TensorValuesTests/OperationsTests.jl index d7fd0b4c3..6c93d91a3 100644 --- a/test/TensorValuesTests/OperationsTests.jl +++ b/test/TensorValuesTests/OperationsTests.jl @@ -104,20 +104,93 @@ c = a + b r = SymTensorValue(6,8,10) @test c==r +c = b - a +r = SymTensorValue(4,4,4) +@test c==r + a = TensorValue(1,2,3,4) b = SymTensorValue(5,6,7) c = a + b +d = b + a r = TensorValue(6,8,9,11) @test c==r +@test d==r + +c = a - b +r = TensorValue(-4,-4,-3,-3) +@test c==r + +c = b - a +r = TensorValue(4,4,3,3) +@test c==r + +a = SymTracelessTensorValue(1,2) +b = SymTracelessTensorValue(5,6) + +c = -a +r = SymTracelessTensorValue(-1,-2) +@test c==r + +c = a + b +r = SymTracelessTensorValue(6,8) +@test c==r + +c = a - b +r = SymTracelessTensorValue(-4,-4) +@test c==r + +a = SymTensorValue(1,2,3) +b = SymTracelessTensorValue(5,6) + +c = a + b +d = b + a +r = SymTensorValue(6,8,-2) +@test c==r +@test d==r + +c = a - b +r = SymTensorValue(-4,-4,8) +@test c==r + +c = b - a +r = SymTensorValue(4,4,-8) +@test c==r a = SymTensorValue(5,6,7) b = TensorValue(1,2,3,4) +c = a + b +d = b + a +r = TensorValue(6,8,9,11) +@test c==r +@test d==r + c = a - b r = TensorValue(4,4,3,3) @test c==r +c = b - a +r = TensorValue(-4,-4,-3,-3) +@test c==r + +a = SymTracelessTensorValue(5,6) +b = TensorValue(1,2,3,4) + +c = a + b +d = b + a +r = TensorValue(6,8,9,-1) +@test c==r +@test d==r + +c = a - b +r = TensorValue(4,4,3,-9) +@test c==r + +c = b - a +r = TensorValue(-4,-4,-3,9) +@test c==r + # Matrix Division a = VectorValue(1,2,3) @@ -134,6 +207,7 @@ c = st\a t = TensorValue(1,2,3,4,5,6,7,8,9) st = SymTensorValue(1,2,3,5,6,9) +qt = SymTracelessTensorValue(1,2,3,5,6) s4ot = one(SymFourthOrderTensorValue{2,Int}) a = VectorValue(1,2,3) @@ -192,6 +266,18 @@ c = st + 2 r = SymTensorValue(3,4,5,7,8,11) @test c == r + +c = 2 * qt +@test isa(c,SymTracelessTensorValue{3}) +r = SymTracelessTensorValue(2,4,6,10,12) +@test c == r + +c = qt * 2 +@test isa(c,SymTracelessTensorValue{3}) +r = SymTracelessTensorValue(2,4,6,10,12) +@test c == r + + c = 2 * s4ot @test isa(c,SymFourthOrderTensorValue{2}) r = SymFourthOrderTensorValue(2,0,0, 0,1,0, 0,0,2) @@ -214,8 +300,10 @@ b = VectorValue(2,1,6) t = TensorValue(1,2,3,4,5,6,7,8,9) s = TensorValue(9,8,3,4,5,6,7,2,1) -st = SymTensorValue(1,2,3,5,6,9) +st = SymTensorValue(1,2,3,5,6,9) st2 = SymTensorValue(9,6,5,3,2,1) +qt = SymTracelessTensorValue(1,2,3,5,6) +qt2 = SymTracelessTensorValue(9,6,5,3,2) c = a ⋅ b @test isa(c,Int) @@ -231,6 +319,11 @@ c = st ⋅ a r = VectorValue(14,30,42) @test c == r +c = qt ⋅ a +@test isa(c,VectorValue{3,Int}) +r = VectorValue(14,30,-3) +@test c == r + c = s ⋅ t @test isa(c,TensorValue{3,3,Int}) r = TensorValue(38,24,18,98,69,48,158,114,78) @@ -241,11 +334,31 @@ c = st ⋅ st2 r = TensorValue(36, 78, 108, 18, 39, 54, 12, 26, 36) @test c == r +c = qt ⋅ qt2 +@test isa(c,TensorValue{3,3,Int}) +r = TensorValue(36, 78, 33, 18, 39, 24, -27, -52, 99) +@test c == r + +c = st ⋅ qt2 +@test isa(c,TensorValue{3,3,Int}) +r = TensorValue(36, 78, 108, 18, 39, 54, -27, -52, -81) +@test c == r + +c = qt2 ⋅ st +@test isa(c,TensorValue{3,3,Int}) +r = TensorValue(36, 18, -27, 78, 39, -52, 108, 54, -81) +@test c == r + c = a ⋅ st @test isa(c,VectorValue{3,Int}) r = VectorValue(14,30,42) @test c == r +c = a ⋅ qt +@test isa(c,VectorValue{3,Int}) +r = VectorValue(14,30,-3) +@test c == r + a1 = VectorValue(1,0) b1 = VectorValue(1,2) @@ -307,6 +420,11 @@ c = st ⊙ st2 @test isa(c,Int) @test c == inner(TensorValue(get_array(st)),TensorValue(get_array(st2))) +c = inner(qt,qt2) +c = qt ⊙ qt2 +@test isa(c,Int) +@test c == inner(TensorValue(get_array(qt)),TensorValue(get_array(qt2))) + # Reductions a = VectorValue(1,2,3) @@ -389,7 +507,11 @@ c = inv(t) st = SymTensorValue(9,8,7,5,4,1) @test det(st) == det(TensorValue(get_array(st))) -@test inv(st) == inv(TensorValue(get_array(st))) +@test inv(st) ≈ inv(TensorValue(get_array(st))) + +qt = SymTracelessTensorValue(9,8,7,5,4) +@test det(qt) == det(TensorValue(get_array(qt))) +@test inv(qt) ≈ inv(TensorValue(get_array(qt))) t = TensorValue(10) @test det(t) == 10 @@ -412,6 +534,9 @@ c = meas(t) st = SymTensorValue(1,2,3,5,6,9) @test meas(st) == meas(TensorValue(get_array(st))) +qt = SymTracelessTensorValue(1,2,3,5,6) +@test meas(qt) == meas(TensorValue(get_array(qt))) + v = TensorValue{1,2}(10,20) @test meas(v) == sqrt(500) @@ -465,8 +590,14 @@ t = TensorValue(1,2,3,4,5,6,7,8,9) st = SymTensorValue(1,2,3,5,6,9) @test tr(st) == tr(TensorValue(get_array(st))) +qt = SymTracelessTensorValue(1,2,3,5,6) +@test tr(qt) == tr(TensorValue(get_array(qt))) + @test get_array(symmetric_part(t)) == get_array(TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0)) @test symmetric_part(st) == symmetric_part(TensorValue(get_array(st))) +@test symmetric_part(st) === st +@test symmetric_part(qt) == symmetric_part(TensorValue(get_array(qt))) +@test symmetric_part(qt) === qt a = TensorValue(1,2,3,4) b = a' @@ -566,6 +697,15 @@ odot_contraction_array = 1*a[:,1,1] + 2*a[:,1,2] + 3*a[:,1,3] + 2*a[:,2,1] + 4*a[:,2,2] + 5*a[:,2,3] + 3*a[:,3,1] + 5*a[:,3,2] + 6*a[:,3,3] @test odot_contraction == odot_contraction_array +a = reshape(Vector(1:27),(3,3,3)) +a_tensor = ThirdOrderTensorValue(a...) +b_tensor = SymTracelessTensorValue((1:5)...) +b = Matrix(get_array(b_tensor)) +odot_contraction = Vector(get_array(a_tensor ⋅² b_tensor)) +odot_contraction_array = 1*a[:,1,1] + 2*a[:,1,2] + 3*a[:,1,3] + 2*a[:,2,1] + + 4*a[:,2,2] + 5*a[:,2,3] + 3*a[:,3,1] + 5*a[:,3,2] + (-5)*a[:,3,3] +@test odot_contraction == odot_contraction_array + # double Contractions w/ products Sym4TensorIndexing = [1111, 1121, 1131, 1122, 1132, 1133, 2111, 2121, 2131, 2122, 2132, 2133, 3111, 3121, 3131, 3122, 3132, 3133, 2211, 2221, 2231, 2222, 2232, 2233, @@ -738,4 +878,22 @@ c = a .* b @test diag(a) == VectorValue(1,4) +a = SymTensorValue(1,2,4) +b = SymTensorValue(1.,2.,4.) +c = a .* b +@test isa(c,SymTensorValue) +@test c.data == map(*,a.data,b.data) + +@test diag(a) == VectorValue(1,4) + + +# Componant wise operations on sym. traceless tensors yield sym. tensors +a = SymTracelessTensorValue(1,2) +b = SymTracelessTensorValue(1.,2.) +c = a .* b +@test isa(c,SymTensorValue) +@test c.data == map(*,a.data,b.data) + +@test diag(a) == VectorValue(1,-1) + end # module OperationsTests diff --git a/test/TensorValuesTests/ReinterpretTests.jl b/test/TensorValuesTests/ReinterpretTests.jl index bba4228af..23cf505a0 100644 --- a/test/TensorValuesTests/ReinterpretTests.jl +++ b/test/TensorValuesTests/ReinterpretTests.jl @@ -27,4 +27,10 @@ A = reinterpret(V) R = [1 1 1 1 1; 2 2 2 2 2; 3 3 3 3 3] @test A == R +v = SymTracelessTensorValue(1,2) +V = fill(v,5) +A = reinterpret(V) +R = [1 1 1 1 1; 2 2 2 2 2; -1 -1 -1 -1 -1] +@test A == R + end # module ReinterpretTests diff --git a/test/TensorValuesTests/TypesTests.jl b/test/TensorValuesTests/TypesTests.jl index c86d8391e..900d22231 100644 --- a/test/TensorValuesTests/TypesTests.jl +++ b/test/TensorValuesTests/TypesTests.jl @@ -103,6 +103,60 @@ s = SymTensorValue{2,Int}(11,21.0,22) @test isa(s,SymTensorValue{2,Int}) @test convert(SMatrix{2,2,Int},s) == [11.0 21.0;21.0 22.0] +# Constructors (SymTracelessTensorValue) + +q = SymTracelessTensorValue( (11,21) ) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Int},q) == [11 21;21 -11] + +q = SymTracelessTensorValue(11,21) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Float64},q) == [11.0 21.0;21.0 -11.0] + +q = SymTracelessTensorValue{2}( (11,21) ) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Int},q) == [11 21;21 -11] + +q = SymTracelessTensorValue{2}(11,21) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Float64},q) == [11.0 21.0;21.0 -11.0] + +q = SymTracelessTensorValue{2,Int}( (11,21) ) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Int},q) == [11 21;21 -11] + +q = SymTracelessTensorValue{2,Float64}(11,21) +@test isa(q,SymTracelessTensorValue{2,Float64}) +@test convert(SMatrix{2,2,Float64},q) == [11.0 21.0;21.0 -11.0] + +q = SymTracelessTensorValue{0,Int}( () ) +@test isa(q,SymTracelessTensorValue{0,Int}) +@test convert(SMatrix{0,0,Int},q) == Array{Any,2}(undef,0,0) + +q = SymTracelessTensorValue{0,Int}() +@test isa(q,SymTracelessTensorValue{0,Int}) +@test convert(SMatrix{0,0,Int},q) == Array{Any,2}(undef,0,0) + +q = SymTracelessTensorValue{1,Int}( () ) +@test isa(q,SymTracelessTensorValue{1,Int}) +@test convert(SMatrix{1,1,Int},q) == [0;;] + +q = SymTracelessTensorValue{1,Int}() +@test isa(q,SymTracelessTensorValue{1,Int}) +@test convert(SMatrix{1,1,Int},q) == [0;;] + +q = SymTracelessTensorValue(11,21.0) +@test isa(q,SymTracelessTensorValue{2,Float64}) +@test convert(SMatrix{2,2,Float64},q) == [11.0 21.0;21.0 -11.0] + +q = SymTracelessTensorValue{2}(11,21.0) +@test isa(q,SymTracelessTensorValue{2,Float64}) +@test convert(SMatrix{2,2,Float64},q) == [11.0 21.0;21.0 -11.0] + +q = SymTracelessTensorValue{2,Int}(11,21.0) +@test isa(q,SymTracelessTensorValue{2,Int}) +@test convert(SMatrix{2,2,Int},q) == [11.0 21.0;21.0 -11.0] + # Constructors (SymFourthOrderTensorValue) s = SymFourthOrderTensorValue( (1111,1121,1122, 2111,2121,2122, 2211,2221,2222) ) @@ -241,6 +295,10 @@ z = zero(SymTensorValue{3,Int}) @test isa(z,SymTensorValue{3,Int,6}) @test convert(SMatrix{3,3,Int},z) == zeros(Int,(3,3)) +z = zero(SymTracelessTensorValue{3,Int}) +@test isa(z,SymTracelessTensorValue{3,Int,6}) +@test convert(SMatrix{3,3,Int},z) == zeros(Int,(3,3)) + z = zero(ThirdOrderTensorValue{3,3,3,Int,27}) @test isa(z,ThirdOrderTensorValue{3,3,3,Int,27}) @test Tuple(z) == Tuple(zeros(Int,(27))) @@ -279,6 +337,10 @@ r = rand(SymTensorValue{3,Int}) @test isa(r,SymTensorValue{3,Int,6}) @test r ≠ rand(typeof(r)) +r = rand(SymTracelessTensorValue{3,Int}) +@test isa(r,SymTracelessTensorValue{3,Int,6}) +@test r ≠ rand(typeof(r)) + r = rand(SymFourthOrderTensorValue{3,Int}) @test isa(r,SymFourthOrderTensorValue{3,Int,36}) @test r ≠ rand(typeof(r)) @@ -311,6 +373,13 @@ b = convert(V,a) b = V[a,a,a,] @test isa(b,Vector{V}) +a = (11,21) +V = SymTracelessTensorValue{2,Int,3} +b = convert(V,a) +@test isa(b,V) +b = V[a,a,a] +@test isa(b,Vector{V}) + a = (1111,1121,1122, 2111,2121,2122, 2211,2221,2222) V = SymFourthOrderTensorValue{2,Int,9} b = convert(V,a) @@ -403,6 +472,7 @@ v = VectorValue(m) @test num_components(VectorValue(1,2,3)) == 3 @test num_components(TensorValue(1,2,3,4)) == 4 @test num_components(SymTensorValue(1,2,3)) == 4 +@test num_components(SymTracelessTensorValue(1,2)) == 4 @test num_components(SymFourthOrderTensorValue(1111,1121,1122, 2111,2121,2122, 2211,2221,2222)) == 16 a = VectorValue(1,2,3,4) @@ -434,9 +504,21 @@ b[1,1] = a[1,1] b[1,2] = a[1,2] b[2,1] = a[2,1] b[2,2] = a[2,2] +a = SymTensorValue(11,21,22) bt = SymTensorValue{2,Int64}(b) -@test bt .== a +@test all(bt .== a) +a = SymTracelessTensorValue(11,21) +@test change_eltype(a,Float64) == SymTracelessTensorValue{2,Float64,3} +@test isa(Tuple(a),Tuple) +@test Tuple(a) == a.data +b = Matrix{Int64}(undef,2,2) +b[1,1] = a[1,1] +b[1,2] = a[1,2] +b[2,1] = a[2,1] +b[2,2] = a[2,2] +bt = SymTracelessTensorValue{2,Int64}(b) +@test all(bt .== a) a = SymFourthOrderTensorValue(1111,1121,1122, 2111,2121,2122, 2211,2221,2222) @test change_eltype(a,Float64) == SymFourthOrderTensorValue{2,Float64,9}