Export rotation_tensor and add rotation implementations for 2D (#185)

This patch renames the internal rotation_matrix function to
rotation_tensor and exports it. In addition, this patch includes
implementations for rotate and rotation_tensor for 2D tensors, where it
is assumed that the axis of rotation is the out-of-plane axis.

Project.toml

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

src/Tensors.jl

@@ -22,7 +22,7 @@ export tdot, dott, dotdot
 export hessian, gradient, curl, divergence, laplace
 export @implement_gradient
 export basevec, eᵢ
-export rotate
+export rotate, rotation_tensor
 export tovoigt, tovoigt!, fromvoigt, tomandel, tomandel!, frommandel
 #########
 # Types #

src/math_ops.jl

@@ -289,17 +289,9 @@ Rotate a three dimensional tensor `x` around the vector `u` a total of `θ` radi
 
 # Examples
 ```jldoctest
-julia> x = Vec{3}((0.0, 0.0, 1.0))
-3-element Vec{3, Float64}:
- 0.0
- 0.0
- 1.0
+julia> x = Vec{3}((0.0, 0.0, 1.0));
 
-julia> u = Vec{3}((0.0, 1.0, 0.0))
-3-element Vec{3, Float64}:
- 0.0
- 1.0
- 0.0
+julia> u = Vec{3}((0.0, 1.0, 0.0));
 
 julia> rotate(x, u, π/2)
 3-element Vec{3, Float64}:
@@ -308,37 +300,79 @@ julia> rotate(x, u, π/2)
  6.123233995736766e-17
 ```
 """
-rotate(x::AbstractTensor, u::Vec{3}, θ::Number)
+rotate(x::AbstractTensor{<:Any,3}, u::Vec{3}, θ::Number)
 
 function rotate(x::Vec{3}, u::Vec{3}, θ::Number)
     ux = u ⋅ x
     u² = u ⋅ u
-    c = cos(θ)
-    s = sin(θ)
+    s, c = sincos(θ)
     (u * ux * (1 - c) + u² * x * c + sqrt(u²) * (u × x) * s) / u²
 end
 
-function rotation_matrix(u::Vec{3, T}, θ::Number) where T
+"""
+    rotate(x::AbstractTensor{2}, θ::Number)
+
+Rotate a two dimensional tensor `x` `θ` radians around the out-of-plane axis.
+
+# Examples
+```jldoctest
+julia> x = Vec{2}((0.0, 1.0));
+
+julia> rotate(x, π/4)
+2-element Vec{2, Float64}:
+ -0.7071067811865475
+  0.7071067811865476
+```
+"""
+rotate(x::AbstractTensor{<:Any, 2}, θ::Number)
+
+function rotate(x::Vec{2}, θ::Number)
+    s, c = sincos(θ)
+    return Vec{2}((c * x[1] - s * x[2], s * x[1] + c * x[2]))
+end
+
+@deprecate rotation_matrix rotation_tensor false
+
+"""
+    rotation_tensor(θ::Number)
+
+Return the two-dimensional rotation matrix corresponding to rotation of `θ` radians around
+the out-of-plane axis (i.e. around `(0, 0, 1)`).
+"""
+function rotation_tensor(θ::Number)
+    s, c = sincos(θ)
+    return Tensor{2, 2}((c, s, -s, c))
+end
+
+"""
+    rotation_tensor(u::Vec{3}, θ::Number)
+
+Return the three-dimensional rotation matrix corresponding to rotation of `θ` radians around
+the vector `u`.
+"""
+function rotation_tensor(u::Vec{3, T}, θ::Number) where T
     # See http://mathworld.wolfram.com/RodriguesRotationFormula.html
     u = u / norm(u)
     z = zero(T)
     ω = Tensor{2, 3}((z, u[3], -u[2], -u[3], z, u[1], u[2], -u[1], z))
-    return one(ω) + sin(θ) * ω + (1 - cos(θ)) * ω^2
+    s, c = sincos(θ)
+    return one(ω) + s * ω + (1 - c) * ω^2
 end
 
-function rotate(x::SymmetricTensor{2, 3}, u::Vec{3}, θ::Number)
-    R = rotation_matrix(u, θ)
+# args is (u::Vec{3}, θ::Number) for 3D tensors, and (θ::number,) for 2D
+function rotate(x::SymmetricTensor{2}, args...)
+    R = rotation_tensor(args...)
     return unsafe_symmetric(R ⋅ x ⋅ R')
 end
-function rotate(x::Tensor{2, 3}, u::Vec{3}, θ::Number)
-    R = rotation_matrix(u, θ)
+function rotate(x::Tensor{2}, args...)
+    R = rotation_tensor(args...)
     return R ⋅ x ⋅ R'
 end
-function rotate(x::Tensor{4, 3}, u::Vec{3}, θ::Number)
-    R = rotation_matrix(u, θ)
+function rotate(x::Tensor{4}, args...)
+    R = rotation_tensor(args...)
     return otimesu(R, R) ⊡ x ⊡ otimesu(R', R')
 end
-function rotate(x::SymmetricTensor{4, 3}, u::Vec{3}, θ::Number)
-    R = rotation_matrix(u, θ)
+function rotate(x::SymmetricTensor{4}, args...)
+    R = rotation_tensor(args...)
     return unsafe_symmetric(otimesu(R, R) ⊡ x ⊡ otimesu(R', R'))
 end

test/F64.jl

@@ -42,6 +42,7 @@ Base.convert(::Type{F64}, a::T) where {T <: Number} = F64(a)
 Base.acos(a::F64) = F64(acos(a.x))
 Base.cos(a::F64) = F64(cos(a.x))
 Base.sin(a::F64) = F64(sin(a.x))
+Base.sincos(a::F64) = (F64(sin(a.x)), F64(cos(a.x)))
 Base.precision(::Type{F64}) = precision(Float64)
 Base.floatmin(::Type{F64}) = floatmin(Float64)
 

test/test_ops.jl

@@ -203,8 +203,8 @@ end # of testsection
     @test AAT ⊡ (b ⊗ a) ≈ (@inferred dotdot(a, AA_sym, b))::Tensor{2, dim, T}
 end # of testsection
 
-if dim == 3
-    @testsection "rotation" begin
+@testsection "rotation" begin
+    if dim == 3
         x = eᵢ(Vec{3, T}, 1)
         y = eᵢ(Vec{3, T}, 2)
         z = eᵢ(Vec{3, T}, 3)
@@ -214,9 +214,11 @@ if dim == 3
         @test rotate(3*z, y, π) ≈ -3*z
         @test rotate(x+y+z, z, π/4) ≈ Vec{3}((0.0,√2,1.0))
 
-        @test rotate(a, b, 0) ≈ a
+        @test rotate(a, b, 0) ≈ a ≈ rotation_tensor(b, 0) ⋅ a
         @test rotate(a, b, π) ≈ rotate(a, b, -π)
+        @test rotate(a, b, π) ≈ rotation_tensor(b, π) ⋅ a rtol=1e-6
         @test rotate(a, b, π/2) ≈ rotate(a, -b, -π/2)
+        @test rotate(a, b, π/2) ≈ rotation_tensor(b, π/2) ⋅ a rtol=1e-6
 
         @test rotate(A, x, π) ≈ rotate(A, x, -π)
         @test rotate(rotate(rotate(A, x, π), y, π), z, π) ≈ A
@@ -240,14 +242,33 @@ if dim == 3
 
         v1, v2, v3, v4, axis = [rand(Vec{3,T}) for _ in 1:5]
         α = rand(T) * π
-        R = Tensors.rotation_matrix(axis, α)
+        R = rotation_tensor(axis, α)
         v1v2v3v4 = (v1 ⊗ v2) ⊗ (v3 ⊗ v4)
         Rv1v2v3v4 = ((R ⋅ v1) ⊗ (R ⋅ v2)) ⊗ ((R ⋅ v3) ⊗ (R ⋅ v4))
         @test rotate(v1v2v3v4, axis, α) ≈ Rv1v2v3v4
         v1v1v2v2 = otimes(v1) ⊗ otimes(v2)
         Rv1v1v2v2 = otimes(R ⋅ v1) ⊗ otimes(R ⋅ v2)
         @test rotate(v1v1v2v2, axis, α) ≈ Rv1v1v2v2
     end
+    if dim == 2
+        θ = T(pi/2)
+        zT = zero(T)
+        z = Vec{3}((zT, zT, one(T)))
+        r = 1:2
+        R = rotation_tensor(θ)
+        R3 = rotation_tensor(z, θ)
+        @test R ≈ R3[r, r]
+        a3 = Vec{3}((a[1], a[2], zT))
+        @test rotate(a, θ) ≈ R ⋅ a ≈ rotate(a3, z, θ)[r]
+        A3 = Tensor{2,3}((i,j) -> i < 3 && j < 3 ? A[i, j] : zT)
+        @test rotate(A, θ) ≈ rotate(A3, z, θ)[r, r]
+        A3s = SymmetricTensor{2,3}((i,j) -> i < 3 && j < 3 ? A_sym[i, j] : zT)
+        @test rotate(A_sym, θ) ≈ rotate(A3s, z, θ)[r, r]
+        AA3 = Tensor{4,3}((i,j,k,l) -> i < 3 && j < 3 && k < 3 && l < 3 ? AA[i, j, k, l] : zT)
+        @test rotate(AA, θ) ≈ rotate(AA3, z, θ)[r, r, r, r]
+        AA3s = SymmetricTensor{4,3}((i,j,k,l) -> i < 3 && j < 3 && k < 3 && l < 3 ? AA_sym[i, j, k, l] : zT)
+        @test rotate(AA_sym, θ) ≈ rotate(AA3s, z, θ)[r, r, r, r]
+    end
 end
 
 @testsection "tovoigt/fromvoigt" begin