Basic support for 3rd order tensors (#205)

CHANGELOG.md

@@ -0,0 +1,26 @@
+# Tensors changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## Unreleased
+
+## [1.16.0]
+
+### Added
+ - Partial support for 3rd order Tensors [#205][github-205]
+    * All construction methods, e.g. `zero(Tensor{3})`, `rand(Tensor{3})`, `Tensor{3}((i,j,k)->f(i,j,k))`
+    * Gradient of 2nd order tensor wrt. vector
+    * `rotate(::Tensor{3})`
+    * `dcontract(::Tensor{D1}, ::Tensor{D2})` for (D1,D2) in ((2,3), (3,2), (3,4), (4,3))
+    * `otimes(::Vec, ::SecondOrderTensor)` and `otimes(::SecondOrderTensor, ::Vec)`
+    * `dot(::Tensor{D1}, ::Tensor{D2})` for (D1,D2) in ((3,1), (1,3), (2,3), (3,2))
+
+<!-- Release links -->
+[Unreleased]: https://github.com/Ferrite-FEM/Tensors.jl/compare/v1.15.1...HEAD
+[1.15.1]: https://github.com/Ferrite-FEM/Tensors.jl/compare/v1.15.0...v1.15.1
+
+<!-- GitHub pull request/issue links -->
+[github-205]: https://github.com/Ferrite-FEM/Tensors.jl/pull/205

Project.toml

@@ -1,6 +1,6 @@
 name = "Tensors"
 uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4"
-version = "1.15.0"
+version = "1.16.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/Tensors.jl

@@ -73,7 +73,7 @@ const Vec{dim, T, M} = Tensor{1, dim, T, dim}
 
 const AllTensors{dim, T} = Union{SymmetricTensor{2, dim, T}, Tensor{2, dim, T},
                                  SymmetricTensor{4, dim, T}, Tensor{4, dim, T},
-                                 Vec{dim, T}}
+                                 Vec{dim, T}, Tensor{3, dim, T}}
 
 
 const SecondOrderTensor{dim, T}   = Union{SymmetricTensor{2, dim, T}, Tensor{2, dim, T}}
@@ -118,16 +118,17 @@ Base.IndexStyle(::Type{<:Tensor}) = IndexLinear()
 ########
 Base.size(::Vec{dim})               where {dim} = (dim,)
 Base.size(::SecondOrderTensor{dim}) where {dim} = (dim, dim)
+Base.size(::Tensor{3,dim})          where {dim} = (dim, dim, dim)
 Base.size(::FourthOrderTensor{dim}) where {dim} = (dim, dim, dim, dim)
 
-# Also define lnegth for the type itself
+# Also define length for the type itself
 Base.length(::Type{Tensor{order, dim, T, M}}) where {order, dim, T, M} = M
 
 #########################
 # Internal constructors #
 #########################
-for TensorType in (SymmetricTensor, Tensor)
-    for order in (2, 4), dim in (1, 2, 3)
+for (TensorType, orders) in ((SymmetricTensor, (2,4)), (Tensor, (2,3,4)))
+    for order in orders, dim in (1, 2, 3)
         N = n_components(TensorType{order, dim})
         @eval begin
             @inline $TensorType{$order, $dim}(t::NTuple{$N, T}) where {T} = $TensorType{$order, $dim, T, $N}(t)

src/automatic_differentiation.jl

@@ -121,6 +121,14 @@ end
     end
     return ∇f
 end
+# Tensor{2} output, Vec input -> Tensor{3} gradient
+@inline function _extract_gradient(v::Tensor{2, 1, <: Dual}, ::Vec{1})
+    @inbounds begin
+        p1 = partials(v[1,1])
+        ∇f = Tensor{3, 1}((p1[1],))
+    end
+    return ∇f
+end
 # Tensor{2} output, Tensor{2} input -> Tensor{4} gradient
 @inline function _extract_gradient(v::Tensor{2, 2, <: Dual}, ::Tensor{2, 2})
     @inbounds begin
@@ -142,6 +150,15 @@ end
     end
     return ∇f
 end
+# Tensor{2} output, Vec input -> Tensor{3} gradient
+@inline function _extract_gradient(v::Tensor{2, 2, <: Dual}, ::Vec{2})
+    @inbounds begin
+        p1, p2, p3, p4 = partials(v[1,1]), partials(v[2,1]), partials(v[1,2]), partials(v[2,2])
+        ∇f = Tensor{3, 2}((p1[1], p2[1], p3[1], p4[1],
+                           p1[2], p2[2], p3[2], p4[2]))
+    end
+    return ∇f
+end
 # Tensor{2} output, Tensor{2} input -> Tensor{4} gradient
 @inline function _extract_gradient(v::Tensor{2, 3, <: Dual}, ::Tensor{2, 3})
     @inbounds begin
@@ -175,6 +192,19 @@ end
     return ∇f
 end
 
+# Tensor{2} output, Vec input -> Tensor{3} gradient
+@inline function _extract_gradient(v::Tensor{2, 3, <: Dual}, ::Vec{3})
+    @inbounds begin
+        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1])
+        p4, p5, p6 = partials(v[1,2]), partials(v[2,2]), partials(v[3,2])
+        p7, p8, p9 = partials(v[1,3]), partials(v[2,3]), partials(v[3,3])
+        ∇f = Tensor{3, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], p7[1], p8[1], p9[1],
+                           p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], p7[2], p8[2], p9[2],  
+                           p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], p7[3], p8[3], p9[3]))
+    end
+    return ∇f
+end
+
 # for non dual variable
 @inline function _extract_value(v::Any)
     return v

src/constructors.jl

@@ -7,6 +7,8 @@
         exp = tensor_create(TensorType, (i) -> :(f($i)))
     elseif order == 2
         exp = tensor_create(TensorType, (i,j) -> :(f($i, $j)))
+    elseif order == 3
+        exp = tensor_create(TensorType, (i,j,k) -> :(f($i, $j, $k)))
     elseif order == 4
         exp = tensor_create(TensorType, (i,j,k,l) -> :(f($i, $j, $k, $l)))
     end

src/indexing.jl

@@ -15,6 +15,13 @@ end
     return dim*(j-1) + i
 end
 
+@inline function compute_index(::Type{Tensor{3, dim}}, i::Int, j::Int, k::Int) where {dim}
+    lower_order = Tensor{2, dim}
+    I = compute_index(lower_order, i, j)
+    n = n_components(lower_order)
+    return (k-1) * n + I
+end
+
 @inline function compute_index(::Type{Tensor{4, dim}}, i::Int, j::Int, k::Int, l::Int) where {dim}
     lower_order = Tensor{2, dim}
     I = compute_index(lower_order, i, j)
@@ -62,8 +69,8 @@ end
 # Slice
 @inline Base.getindex(v::Vec, ::Colon) = v
 
-function Base.getindex(S::Union{SecondOrderTensor, FourthOrderTensor}, ::Colon)
-    throw(ArgumentError("S[:] not defined for S of order 2 or 4, use Array(S) to convert to an Array"))
+function Base.getindex(S::Union{SecondOrderTensor, Tensor{3}, FourthOrderTensor}, ::Colon)
+    throw(ArgumentError("S[:] not defined for S of order 2, 3, or 4, use Array(S) to convert to an Array"))
 end
 
 @inline @generated function Base.getindex(S::SecondOrderTensor{dim}, ::Colon, j::Int) where {dim}

src/math_ops.jl

@@ -21,6 +21,19 @@ julia> norm(A)
 @inline LinearAlgebra.norm(v::Vec) = sqrt(dot(v, v))
 @inline LinearAlgebra.norm(S::SecondOrderTensor) = sqrt(dcontract(S, S))
 
+@generated function LinearAlgebra.norm(S::Tensor{3,dim}) where {dim}
+    idx(i,j,k) = compute_index(get_base(S), i, j, k)
+    ex = Expr[]
+    for k in 1:dim, j in 1:dim, i in 1:dim
+        push!(ex, :(get_data(S)[$(idx(i,j,k))]))
+    end
+    exp = reducer(ex, ex)
+    return quote
+      $(Expr(:meta, :inline))
+      @inbounds return sqrt($exp)
+    end
+end
+
 # special case for Tensor{4, 3} since it is faster than unrolling
 @inline LinearAlgebra.norm(S::Tensor{4, 3}) = sqrt(mapreduce(abs2, +, S))
 
@@ -402,6 +415,10 @@ function rotate(x::Tensor{2}, args...)
     R = rotation_tensor(args...)
     return R ⋅ x ⋅ R'
 end
+function rotate(x::Tensor{3}, args...)
+    R = rotation_tensor(args...)
+    return otimesu(R, R) ⊡ x ⋅ R'
+end
 function rotate(x::Tensor{4}, args...)
     R = rotation_tensor(args...)
     return otimesu(R, R) ⊡ x ⊡ otimesu(R', R')

src/tensor_products.jl

@@ -68,6 +68,74 @@ end
     end
 end
 
+@generated function dcontract(S1::SecondOrderTensor{dim}, S2::Tensor{3,dim}) where {dim}
+    TensorType = getreturntype(dcontract, get_base(S1), get_base(S2))
+    idxS1(i, j) = compute_index(get_base(S1), i, j)
+    idxS2(i, j, k) = compute_index(get_base(S2), i, j, k)
+    exps = Expr(:tuple)
+    for k in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, j))])    for i in 1:dim, j in 1:dim][:]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(i, j, k))]) for i in 1:dim, j in 1:dim][:]
+        push!(exps.args, reducer(ex1, ex2, true))
+    end
+    expr = remove_duplicates(TensorType, exps) # TODO: Required?
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return $TensorType($expr)
+    end
+end
+
+@generated function dcontract(S1::Tensor{3,dim}, S2::SecondOrderTensor{dim}) where {dim}
+    TensorType = getreturntype(dcontract, get_base(S1), get_base(S2))
+    idxS1(i, j, k) = compute_index(get_base(S1), i, j, k)
+    idxS2(i, j) = compute_index(get_base(S2), i, j)
+    exps = Expr(:tuple)
+    for i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, k, l))]) for k in 1:dim, l in 1:dim][:]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(k, l))])    for k in 1:dim, l in 1:dim][:]
+        push!(exps.args, reducer(ex1, ex2, true))
+    end
+    expr = remove_duplicates(TensorType, exps)
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return $TensorType($expr)
+    end
+end
+
+@generated function dcontract(S1::Tensor{3,dim}, S2::FourthOrderTensor{dim}) where {dim}
+    TensorType = getreturntype(dcontract, get_base(S1), get_base(S2))
+    idxS1(i, j, k)    = compute_index(get_base(S1), i, j, k)
+    idxS2(i, j, k, l) = compute_index(get_base(S2), i, j, k, l)
+    exps = Expr(:tuple)
+    for k in 1:dim, j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, m, n))])    for m in 1:dim, n in 1:dim][:]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(m, n, j, k))]) for m in 1:dim, n in 1:dim][:]
+        push!(exps.args, reducer(ex1, ex2, true))
+    end
+    expr = remove_duplicates(TensorType, exps)
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return $TensorType($expr)
+    end
+end
+
+@generated function dcontract(S1::FourthOrderTensor{dim}, S2::Tensor{3,dim}) where {dim}
+    TensorType = getreturntype(dcontract, get_base(S1), get_base(S2))
+    idxS1(i, j, k, l)    = compute_index(get_base(S1), i, j, k, l)
+    idxS2(i, j, k) = compute_index(get_base(S2), i, j, k)
+    exps = Expr(:tuple)
+    for k in 1:dim, j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, j, m, n))]) for m in 1:dim, n in 1:dim][:]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(m, n, k))])    for m in 1:dim, n in 1:dim][:]
+        push!(exps.args, reducer(ex1, ex2, true))
+    end
+    expr = remove_duplicates(TensorType, exps)
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return $TensorType($expr)
+    end
+end
+
 @generated function dcontract(S1::FourthOrderTensor{dim}, S2::FourthOrderTensor{dim}) where {dim}
     TensorType = getreturntype(dcontract, get_base(S1), get_base(S2))
     idxS1(i, j, k, l) = compute_index(get_base(S1), i, j, k, l)
@@ -123,11 +191,20 @@ julia> A ⊗ B
     return Tensor{2, dim}(@inline function(i,j) @inbounds S1[i] * S2[j]; end)
 end
 
+@inline function otimes(S1::Vec{dim}, S2::SecondOrderTensor{dim}) where {dim}
+    return Tensor{3, dim}(@inline function(i,j,k) @inbounds S1[i] * S2[j,k]; end)
+end
+@inline function otimes(S1::SecondOrderTensor{dim}, S2::Vec{dim}) where {dim}
+    return Tensor{3, dim}(@inline function(i,j,k) @inbounds S1[i,j] * S2[k]; end)
+end
+
 @inline function otimes(S1::SecondOrderTensor{dim}, S2::SecondOrderTensor{dim}) where {dim}
     TensorType = getreturntype(otimes, get_base(typeof(S1)), get_base(typeof(S2)))
     TensorType(@inline function(i,j,k,l) @inbounds S1[i,j] * S2[k,l]; end)
 end
 
+# Defining {3}⊗{1} and {1}⊗{3} = {4} would also be valid...
+
 @inline otimes(S1::Number, S2::Number) = S1*S2
 @inline otimes(S1::AbstractTensor, S2::Number) = S1*S2
 @inline otimes(S1::Number, S2::AbstractTensor) = S1*S2
@@ -344,6 +421,64 @@ end
     end
 end
 
+@generated function LinearAlgebra.dot(S1::Tensor{3,dim}, S2::Vec{dim}) where {dim}
+    idxS1(i, j, k) = compute_index(get_base(S1), i, j, k)
+    exps = Expr(:tuple)
+    for j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, j, m))]) for m in 1:dim]
+        ex2 = Expr[:(get_data(S2)[$m])                for m in 1:dim]
+        push!(exps.args, reducer(ex1, ex2))
+    end
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return Tensor{2, dim}($exps)
+    end
+end
+
+@generated function LinearAlgebra.dot(S1::Vec{dim}, S2::Tensor{3,dim}) where {dim}
+    idxS2(i, j, k) = compute_index(get_base(S2), i, j, k)
+    exps = Expr(:tuple)
+    for j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$m]) for m in 1:dim]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(m, i, j))])                for m in 1:dim]
+        push!(exps.args, reducer(ex1, ex2))
+    end
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return Tensor{2, dim}($exps)
+    end
+end
+
+@generated function LinearAlgebra.dot(S1::SecondOrderTensor{dim}, S2::Tensor{3,dim}) where {dim}
+    idxS1(i, j) = compute_index(get_base(S1), i, j)
+    idxS2(i, j, k) = compute_index(get_base(S2), i, j, k)
+    exps = Expr(:tuple)
+    for k in 1:dim, j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, m))]) for m in 1:dim]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(m, j, k))]) for m in 1:dim]
+        push!(exps.args, reducer(ex1, ex2))
+    end
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return Tensor{3, dim}($exps)
+    end
+end
+
+@generated function LinearAlgebra.dot(S1::Tensor{3,dim}, S2::SecondOrderTensor{dim}) where {dim}
+    idxS1(i, j, k) = compute_index(get_base(S1), i, j, k)
+    idxS2(i, j) = compute_index(get_base(S2), i, j)
+    exps = Expr(:tuple)
+    for k in 1:dim, j in 1:dim, i in 1:dim
+        ex1 = Expr[:(get_data(S1)[$(idxS1(i, j, m))]) for m in 1:dim]
+        ex2 = Expr[:(get_data(S2)[$(idxS2(m, k))]) for m in 1:dim]
+        push!(exps.args, reducer(ex1, ex2))
+    end
+    quote
+        $(Expr(:meta, :inline))
+        @inbounds return Tensor{3, dim}($exps)
+    end
+end
+
 """
     dot(::SymmetricTensor{2})
 

src/utilities.jl

@@ -7,6 +7,8 @@ function tensor_create(::Type{Tensor{order, dim}}, f) where {order, dim}
         ex = Expr(:tuple, [f(i) for i=1:dim]...)
     elseif order == 2
         ex = Expr(:tuple, [f(i,j) for i=1:dim, j=1:dim]...)
+    elseif order == 3
+        ex = Expr(:tuple, [f(i,j,k) for i=1:dim, j=1:dim, k=1:dim]...)
     elseif order == 4
         ex = Expr(:tuple, [f(i,j,k,l) for i=1:dim, j=1:dim, k = 1:dim, l = 1:dim]...)
     end
@@ -80,6 +82,10 @@ function dcontract end
 @pure getreturntype(::typeof(dcontract), ::Type{<:SymmetricTensor{4, dim}}, ::Type{<:SecondOrderTensor{dim}}) where {dim} = SymmetricTensor{2, dim}
 @pure getreturntype(::typeof(dcontract), ::Type{<:SecondOrderTensor{dim}}, ::Type{<:Tensor{4, dim}}) where {dim} = Tensor{2, dim}
 @pure getreturntype(::typeof(dcontract), ::Type{<:SecondOrderTensor{dim}}, ::Type{<:SymmetricTensor{4, dim}}) where {dim} = SymmetricTensor{2, dim}
+@pure getreturntype(::typeof(dcontract), ::Type{<:Tensor{3,dim}}, ::Type{<:SecondOrderTensor{dim}}) where {dim} = Vec{dim}
+@pure getreturntype(::typeof(dcontract), ::Type{<:SecondOrderTensor{dim}}, ::Type{<:Tensor{3,dim}}) where {dim} = Vec{dim}
+@pure getreturntype(::typeof(dcontract), ::Type{<:Tensor{3,dim}}, ::Type{<:FourthOrderTensor{dim}}) where {dim} = Tensor{3,dim}
+@pure getreturntype(::typeof(dcontract), ::Type{<:FourthOrderTensor{dim}}, ::Type{<:Tensor{3,dim}}) where {dim} = Tensor{3,dim}
 
 # otimes
 function otimes end

test/test_ad.jl

@@ -117,6 +117,16 @@ S(C) = S(C, μ, Kb)
                 @test ∇∇(A -> T(1), A, :all)[2] ≈ 0*A
                 @test ∇∇(A -> T(1), A, :all)[3] == T(1)
             end
+            @testsection "3rd order" begin
+                δ(i,j) = (i==j ? 1 : 0)
+                #
+                f1(x::Vec) = x⊗x 
+                df1(x::Vec{d}) where d = Tensor{3,d}((i,j,k)-> δ(i,k)*x[j] + δ(j,k)*x[i]);
+                for dim in 1:3
+                    z = rand(Vec{dim})
+                    @test gradient(f1, z) ≈ df1(z)
+                end 
+            end
             end # loop T
         end # testsection
     end # loop dim
@@ -328,5 +338,4 @@ S(C) = S(C, μ, Kb)
 
     end
 
-
 end # testsection