diff --git a/src/basic_operations.jl b/src/basic_operations.jl index 0b5824a2..49277faf 100644 --- a/src/basic_operations.jl +++ b/src/basic_operations.jl @@ -19,22 +19,12 @@ @inline Base.:/(S::AbstractTensor, n::Number) = map(x->(x/n), S) # map implementations -@generated function Base.map{T <: AbstractTensor}(f, S::T) - TensorType = get_base(S) - N = n_components(TensorType) - expr = Expr(:tuple, [:(f(get_data(S)[$i])) for i in 1:N]...) - quote - $(Expr(:meta, :inline)) - @inbounds return $TensorType($expr) - end +@inline function Base.map{T <: AbstractTensor}(f, S::T) + return apply_all(S, @inline function(i) @inboundsret f(S.data[i]); end) end -@generated function Base.map{T <: AllTensors}(f, S1::T, S2::T) - TensorType = get_base(S1) - N = n_components(TensorType) - expr = Expr(:tuple, [:(f(get_data(S1)[$i], get_data(S2)[$i])) for i in 1:N]...) - quote - $(Expr(:meta, :inline)) - @inbounds return $TensorType($expr) - end +@inline Base.map{order, dim, T1, T2}(f, S1::AbstractTensor{order, dim, T1}, S2::AbstractTensor{order, dim, T2}) = ((SS1, SS2) = promote(S1, S2); map(f, SS1, SS2)) + +@inline function Base.map{T <: AllTensors}(f, S1::T, S2::T) + return apply_all(S1, @inline function(i) @inboundsret f(S1.data[i], S2.data[i]); end) end diff --git a/src/constructors.jl b/src/constructors.jl index fdb5137e..3f923b2e 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -16,7 +16,6 @@ end end -# Tensor from AbstractArray @generated function (::Type{Tensor{order, dim}}){order, dim}(data::AbstractArray) N = n_components(Tensor{order,dim}) exp = Expr(:tuple, [:(data[$i]) for i in 1:N]...) @@ -24,11 +23,23 @@ end if length(data) != $N throw(ArgumentError("wrong number of elements, expected $($N), got $(length(data))")) end - # println("BAD DISPATCH") Tensor{order, dim}($exp) end end +@inline function apply_all{order, dim}(S::Union{Tensor{order, dim}, SymmetricTensor{order, dim}}, f::Function) + apply_all(get_base(typeof(S)), f) +end + +# Tensor from AbstractArray +function (T::Type{Tensor{order, dim}}){order, dim}(data::AbstractArray) + N = n_components(T) + length(data) != n_components(T) && throw(ArgumentError("wrong number of elements, expected $N, got $(length(data))")) + return apply_all(T, @inline function(i) @inboundsret data[i]; end) +end + + + # SymmetricTensor from AbstractArray @generated function (::Type{SymmetricTensor{order, dim}}){order, dim}(data::AbstractArray) N = n_components(Tensor{order,dim}) @@ -91,16 +102,8 @@ end @eval @inline Base.$op(t::AllTensors) = $op(typeof(t)) end -@generated function Base.fill{T <: AbstractTensor}(el::Union{Number, Function}, S::Type{T}) - TensorType = get_base(get_type(S)) - N = n_components(TensorType) - ele = el <: Number ? :(el) : :(el()) - expr = Expr(:tuple, [ele for i in 1:N]...) - quote - $(Expr(:meta, :inline)) - @inbounds return $TensorType($expr) - end -end +Base.fill{T <: AbstractTensor}(el::Number, S::Type{T}) = apply_all(get_base(T), i -> el) +Base.fill{T <: AbstractTensor}(f::Function, S::Type{T}) = apply_all(get_base(T), i -> f()) # Array with zero/ones @inline Base.zeros{T <: AbstractTensor}(::Type{T}, dims...) = fill(zero(T), dims...) diff --git a/src/math_ops.jl b/src/math_ops.jl index f6108763..ad5a0823 100644 --- a/src/math_ops.jl +++ b/src/math_ops.jl @@ -263,17 +263,15 @@ julia> trace(dev(A)) 0.0 ``` """ -@generated function dev{dim}(S::SecondOrderTensor{dim}) - Tt = get_base(S) - idx(i,j) = compute_index(Tt, i, j) - f = (i,j) -> i == j ? :((get_data(S)[$(idx(i,j))] - tr/3)) : - :(get_data(S)[$(idx(i,j))]) - exp = tensor_create(Tt, f) - return quote - $(Expr(:meta, :inline)) - tr = trace(S) - $Tt($exp) - end +@inline function dev(S::SecondOrderTensor) + Tt = get_base(typeof(S)) + tr = trace(S) / 3 + Tt( + @inline function(i, j) + @inbounds v = i == j ? S[i,j] - tr : S[i,j] + v + end + ) end # http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/ diff --git a/src/symmetric.jl b/src/symmetric.jl index 84705b63..86c5c63e 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -24,43 +24,28 @@ julia> symmetric(A) """ @inline symmetric(S1::SymmetricTensors) = S1 -@generated function symmetric{dim, T}(S::Tensor{2, dim, T}) - idx(i, j) = compute_index(Tensor{2, dim}, i, j) - expr = Expr(:tuple) - for j in 1:dim, i in j:dim - if i == j - push!(expr.args, :(get_data(S)[$(idx(i, j))])) - else - push!(expr.args, :((get_data(S)[$(idx(i, j))] + get_data(S)[$(idx(j, i))])/2)) - end - end - return quote - $(Expr(:meta, :inline)) - @inbounds return SymmetricTensor{2, dim}($expr) - end +@inline function symmetric{dim}(S::Tensor{2, dim}) + SymmetricTensor{2, dim}(@inline function(i, j) @inboundsret i == j ? S[i,j] : (S[i,j] + S[j,i]) / 2 end) end + + """ ```julia minorsymmetric(::FourthOrderTensor) ``` Computes the minor symmetric part of a fourth order tensor, returns a `SymmetricTensor{4}`. """ -@generated function minorsymmetric{dim}(S::Tensor{4, dim}) - idx(i, j, k, l) = compute_index(Tensor{4, dim}, i, j, k, l) - expr = Expr(:tuple) - for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim - if i == j && k == l - push!(expr.args, :(get_data(S)[$(idx(i, j, k, l))])) - else - push!(expr.args, :((get_data(S)[$(idx(i, j, k, l))] + get_data(S)[$(idx(j, i, k, l))] + - get_data(S)[$(idx(i, j, k, l))] + get_data(S)[$(idx(i, j, l, k))]) / 4)) +@inline function minorsymmetric{dim}(S::Tensor{4, dim}) + SymmetricTensor{4, dim}( + @inline function(i, j, k, l) + @inbounds if i == j && k == l + return S[i,j,k,l] + else + return (S[i,j,k,l] + S[j,i,k,l] + S[i,j,k,l] + S[i,j,l,k]) / 4 + end end - end - return quote - $(Expr(:meta, :inline)) - @inbounds return SymmetricTensor{4, dim}($expr) - end + ) end @inline minorsymmetric(S::SymmetricTensors) = S @@ -73,20 +58,16 @@ majorsymmetric(::FourthOrderTensor) ``` Computes the major symmetric part of a fourth order tensor, returns a `Tensor{4}`. """ -@generated function majorsymmetric{dim}(S::FourthOrderTensor{dim}) - idx(i, j, k, l) = compute_index(get_base(S), i, j, k, l) - expr = Expr(:tuple) - for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim - if i == j == k == l || i == k && j == l - push!(expr.args, :(get_data(S)[$(idx(i, j, k, l))])) - else - push!(expr.args, :((get_data(S)[$(idx(i, j, k, l))] + get_data(S)[$(idx(k, l, i, j))]) / 2)) +@inline function majorsymmetric{dim}(S::FourthOrderTensor{dim}) + Tensor{4, dim}( + @inline function(i, j, k, l) + @inbounds if i == j == k == l || i == k && j == l + return S[i,j,k,l] + else + return (S[i,j,k,l] + S[k,l,i,j]) / 2 + end end - end - return quote - $(Expr(:meta, :inline)) - @inbounds return Tensor{4, dim}($expr) - end + ) end """ @@ -100,17 +81,10 @@ Computes the skew-symmetric (anti-symmetric) part of a second order tensor, retu # Symmetry checks @inline Base.issymmetric(t::Tensor{2, 1}) = true -@inline function Base.issymmetric(t::Tensor{2, 2}) - data = get_data(t) - @inbounds return data[2] == data[3] -end +@inline Base.issymmetric(t::Tensor{2, 2}) = @inboundsret t[1,2] == t[2,1] + @inline function Base.issymmetric(t::Tensor{2, 3}) - data = get_data(t) - @inbounds begin - return (data[2] == data[4] && - data[3] == data[7] && - data[6] == data[8]) - end + return @inboundsret t[1,2] == t[2,1] && t[1,3] == t[3,1] && t[2,3] == t[3,2] end function isminorsymmetric{dim}(t::Tensor{4, dim}) diff --git a/src/tensor_ops_errors.jl b/src/tensor_ops_errors.jl index 7a272009..8fba80a1 100644 --- a/src/tensor_ops_errors.jl +++ b/src/tensor_ops_errors.jl @@ -2,22 +2,22 @@ # Remove `*` as infix operator between tensors function Base.:*(S1::AbstractTensor, S2::AbstractTensor) - error("Don't use `*` for multiplication between tensors. Use `⋅` (`\\cdot`) for single contraction and `⊡` (`\\boxdot`) for double contraction.") + error("use `⋅` (`\\cdot`) for single contraction and `⊡` (`\\boxdot`) for double contraction instead of `*`") end # Remove `'*` as infix operator between tensors function Base.Ac_mul_B(S1::AbstractTensor, S2::AbstractTensor) - error("Don't use `A'*B`, use `tdot(A,B)` (or `A'⋅B`) instead.") + error("use `tdot(A,B)` (or `A'⋅B`) instead of `A'*B`") end # Remove `.'*` as infix operator between tensors function Base.At_mul_B(S1::AbstractTensor, S2::AbstractTensor) - error("Don't use `A.'*B`, use `tdot(A,B)` (or `A.'⋅B`) instead.") + error("use `tdot(A,B)` (or `A.'⋅B`) instead if A.'*B") end # Remove `\` as infix operator between tensors function Base.:\(S1::AbstractTensor, S2::AbstractTensor) - error("Don't use `A\\B`, use `inv(A) ⋅ B` instead.") + error("use `inv(A) ⋅ B` instead of A\\B") end # Remove + and - between number and Tensor (issue #75) diff --git a/src/tensor_products.jl b/src/tensor_products.jl index 43925f7d..5bd3264b 100644 --- a/src/tensor_products.jl +++ b/src/tensor_products.jl @@ -123,27 +123,13 @@ julia> A ⊗ B 0.654957 0.48365 ``` """ -@generated function otimes{dim}(S1::Vec{dim}, S2::Vec{dim}) - exps = Expr(:tuple, [:(get_data(S1)[$i] * get_data(S2)[$j]) for i in 1:dim, j in 1:dim]...) - quote - $(Expr(:meta, :inline)) - @inbounds return Tensor{2, dim}($exps) - end +@inline function otimes{dim}(S1::Vec{dim}, S2::Vec{dim}) + return Tensor{2, dim}(@inline function(i,j) @inboundsret S1[i] * S2[j]; end) end -@generated function otimes{dim}(S1::SecondOrderTensor{dim}, S2::SecondOrderTensor{dim}) - TensorType = getreturntype(otimes, get_base(S1), get_base(S2)) - idxS1(i, j) = compute_index(get_base(S1), i, j) - idxS2(i, j) = compute_index(get_base(S2), i, j) - exps = Expr(:tuple) - for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim - push!(exps.args, :(get_data(S1)[$(idxS1(i, j))] * get_data(S2)[$(idxS2(k, l))])) - end - expr = remove_duplicates(TensorType, exps) - quote - $(Expr(:meta, :inline)) - @inbounds return $TensorType($expr) - end +@inline function otimes{dim}(S1::SecondOrderTensor{dim}, S2::SecondOrderTensor{dim}) + TensorType = getreturntype(otimes, get_base(typeof(S1)), get_base(typeof(S2))) + TensorType(@inline function(i,j,k,l) @inboundsret S1[i,j] * S2[k,l]; end) end const ⊗ = otimes diff --git a/src/transpose.jl b/src/transpose.jl index 0a12ae2d..8ca619bd 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -24,15 +24,8 @@ julia> A' """ @inline Base.transpose(S::Vec) = S -@generated function Base.transpose{dim}(S::Tensor{2, dim}) - expr = Expr(:tuple) - for j in 1:dim, i in 1:dim - push!(expr.args, :(get_data(S)[$(compute_index(Tensor{2, dim}, j, i))])) - end - quote - $(Expr(:meta, :inline)) - @inbounds return Tensor{2, dim}($expr) - end +@inline function Base.transpose{dim}(S::Tensor{2, dim}) + Tensor{2, dim}(@inline function(i, j) @inbounds v = S[j,i]; v; end) end Base.transpose(S::SymmetricTensor{2}) = S @@ -43,15 +36,8 @@ minortranspose(::FourthOrderTensor) ``` Computes the minor transpose of a fourth order tensor. """ -@generated function minortranspose{dim}(S::Tensor{4, dim}) - expr = Expr(:tuple) - for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim - push!(expr.args, :(get_data(S)[$(compute_index(Tensor{4, dim}, j, i, l, k))])) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return Tensor{4, dim}($expr) - end +@inline function minortranspose{dim}(S::Tensor{4, dim}) + Tensor{4, dim}(@inline function(i, j, k, l) @inboundsret S[j,i,l,k]; end) end minortranspose(S::SymmetricTensor{4}) = S @@ -63,15 +49,8 @@ majortranspose(::FourthOrderTensor) ``` Computes the major transpose of a fourth order tensor. """ -@generated function majortranspose{dim}(S::FourthOrderTensor{dim}) - expr = Expr(:tuple) - for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim - push!(expr.args, :(get_data(S)[$(compute_index(get_base(S), k, l, i, j))])) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return Tensor{4, dim}($expr) - end +@inline function majortranspose{dim}(S::FourthOrderTensor{dim}) + Tensor{4, dim}(@inline function(i, j, k, l) @inboundsret S[k,l,i,j]; end) end Base.ctranspose(S::AllTensors) = transpose(S) diff --git a/src/utilities.jl b/src/utilities.jl index 4242e45b..499db2d4 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,3 +1,14 @@ +macro inboundsret(ex) + return quote + @inbounds v = $(esc(ex)) + v + end +end + +function tensor_create_linear{order, dim}(T::Union{Type{Tensor{order, dim}}, Type{SymmetricTensor{order, dim}}}, f) + return Expr(:tuple, [f(i) for i=1:n_components(T)]...) +end + function tensor_create{order, dim}(::Type{Tensor{order, dim}}, f) if order == 1 ex = Expr(:tuple, [f(i) for i=1:dim]...) @@ -70,23 +81,23 @@ end # return types # double contraction function dcontract end -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{Tensor{4, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{SymmetricTensor{4, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{Tensor{4, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{4, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{Tensor{2, dim}}) = Tensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{SymmetricTensor{2, dim}}) = Tensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{Tensor{2, dim}}) = SymmetricTensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{SymmetricTensor{2, dim}}) = SymmetricTensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{2, dim}}, ::Type{Tensor{4, dim}}) = Tensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{2, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{2, dim}}, ::Type{Tensor{4, dim}}) = Tensor{2, dim} -getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{2, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{Tensor{4, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{SymmetricTensor{4, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{Tensor{4, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{4, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{Tensor{2, dim}}) = Tensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{4, dim}}, ::Type{SymmetricTensor{2, dim}}) = Tensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{Tensor{2, dim}}) = SymmetricTensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{4, dim}}, ::Type{SymmetricTensor{2, dim}}) = SymmetricTensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{2, dim}}, ::Type{Tensor{4, dim}}) = Tensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{Tensor{2, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{2, dim}}, ::Type{Tensor{4, dim}}) = Tensor{2, dim} +@pure getreturntype{dim}(::typeof(dcontract), ::Type{SymmetricTensor{2, dim}}, ::Type{SymmetricTensor{4, dim}}) = SymmetricTensor{2, dim} # otimes function otimes end -getreturntype{dim}(::typeof(otimes), ::Type{Tensor{1, dim}}, ::Type{Tensor{1, dim}}) = Tensor{2, dim} -getreturntype{dim}(::typeof(otimes), ::Type{Tensor{2, dim}}, ::Type{Tensor{2, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(otimes), ::Type{SymmetricTensor{2, dim}}, ::Type{Tensor{2, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(otimes), ::Type{Tensor{2, dim}}, ::Type{SymmetricTensor{2, dim}}) = Tensor{4, dim} -getreturntype{dim}(::typeof(otimes), ::Type{SymmetricTensor{2, dim}}, ::Type{SymmetricTensor{2, dim}}) = SymmetricTensor{4, dim} +@pure getreturntype{dim}(::typeof(otimes), ::Type{Tensor{1, dim}}, ::Type{Tensor{1, dim}}) = Tensor{2, dim} +@pure getreturntype{dim}(::typeof(otimes), ::Type{Tensor{2, dim}}, ::Type{Tensor{2, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(otimes), ::Type{SymmetricTensor{2, dim}}, ::Type{Tensor{2, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(otimes), ::Type{Tensor{2, dim}}, ::Type{SymmetricTensor{2, dim}}) = Tensor{4, dim} +@pure getreturntype{dim}(::typeof(otimes), ::Type{SymmetricTensor{2, dim}}, ::Type{SymmetricTensor{2, dim}}) = SymmetricTensor{4, dim} diff --git a/test/test_ops.jl b/test/test_ops.jl index 646b165c..353603ae 100644 --- a/test/test_ops.jl +++ b/test/test_ops.jl @@ -7,7 +7,7 @@ end @testsection "tensor ops" begin for T in (Float32, Float64, F64), dim in (1,2,3) - +println(T) AA = rand(Tensor{4, dim, T}) BB = rand(Tensor{4, dim, T}) A = rand(Tensor{2, dim, T}) @@ -250,6 +250,14 @@ end # of testsection @test rotate(a, b, 0) ≈ a @test rotate(a, b, π) ≈ rotate(a, b, -π) @test rotate(a, b, π/2) ≈ rotate(a, -b, -π/2) +<<<<<<< HEAD end end # of testsection end # of testsection +======= +end # of testset +end +println() +end +end # of testset +>>>>>>> kc/simplify