merge

Ferrite-FEM · Mar 18, 2017 · 48b5acf · 48b5acf
1 parent 7728952
commit 48b5acf
Show file tree

Hide file tree

Showing 9 changed files with 107 additions and 158 deletions.
diff --git 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
@@ -16,19 +16,30 @@
     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]...)
     return quote
         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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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}