diff --git a/Project.toml b/Project.toml index a4fc989f..4666f9e6 100644 --- a/Project.toml +++ b/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" diff --git a/src/Tensors.jl b/src/Tensors.jl index 8bc8bc58..6088d80d 100644 --- a/src/Tensors.jl +++ b/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 # diff --git a/src/math_ops.jl b/src/math_ops.jl index a0beb976..c8f274d9 100644 --- a/src/math_ops.jl +++ b/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 diff --git a/test/F64.jl b/test/F64.jl index 2e5118f2..d639bbab 100644 --- a/test/F64.jl +++ b/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) diff --git a/test/test_ops.jl b/test/test_ops.jl index d143adab..c457fdce 100644 --- a/test/test_ops.jl +++ b/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,7 +242,7 @@ 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 @@ -248,6 +250,25 @@ if dim == 3 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