add rotation (#108)

docs/src/man/other_operators.md

@@ -144,6 +144,12 @@ where ``\lambda_i`` are the eigenvalues and ``\mathbf{v}_i`` are the correspondi
 eig
 ```
 
+## Rotations
+
+```@docs
+rotate
+```
+
 ## Special operations
 
 For computing a special dot product between two vectors $\mathbf{a}$ and $\mathbf{b}$ with a fourth order symmetric tensor $\mathbf{C}$ such that $a_k C_{ikjl} b_l$ there is `dotdot(a, C, b)`. This function is useful because it is the expression for the tangent matrix in continuum mechanics when the displacements are approximated by scalar shape functions.

src/ContMechTensors.jl

@@ -12,6 +12,8 @@ export minortranspose, majortranspose, isminorsymmetric, ismajorsymmetric
 export tdot, dott, dotdot
 export hessian#, gradient
 export basevec, eᵢ
+export rotate
+
 
 @deprecate extract_components(tensor) Array(tensor)
 

src/math_ops.jl

@@ -283,3 +283,42 @@ julia> trace(dev(A))
         $Tt($exp)
     end
 end
+
+# http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
+"""
+Rotate a three dimensional vector `x` around another vector `u` a total of `θ` radians.
+
+```julia
+rotate(x::Vec{3}, u::Vec{3}, θ::Number)
+```
+
+**Example:**
+
+```jldoctest
+julia> x = Vec{3}((0.0, 0.0, 1.0))
+3-element ContMechTensors.Tensor{1,3,Float64,3}:
+ 0.0
+ 0.0
+ 1.0
+
+julia> u = Vec{3}((0.0, 1.0, 0.0))
+3-element ContMechTensors.Tensor{1,3,Float64,3}:
+ 0.0
+ 1.0
+ 0.0
+
+julia> rotate(x, u, π/2)
+3-element ContMechTensors.Tensor{1,3,Float64,3}:
+ 1.0
+ 0.0
+ 6.12323e-17
+```
+"""
+function rotate(x::Vec{3}, u::Vec{3}, θ::Number)
+    ux = u ⋅ x
+    u² = u ⋅ u
+    c = cos(θ)
+    s = sin(θ)
+    (u * ux * (1 - c) + u² * x * c + sqrt(u²) * (u × x) * s) / u²
+end
+

test/test_ops.jl

@@ -233,5 +233,23 @@ end # of testset
     AAT = Tensor{4, dim, T}((i,j,k,l) -> AA_sym[i,l,k,j])
     @test AAT ⊡ (b ⊗ a) ≈ dotdot(a, AA_sym, b)
 end # of testset
+
+@testset "rotation" begin
+    x = eᵢ(Vec{3}, 1)
+    y = eᵢ(Vec{3}, 2)
+    z = eᵢ(Vec{3}, 3)
+
+    @test rotate(z, z, rand()) ≈ z
+    @test rotate(2*z, y, π/2) ≈ 2*x
+    @test rotate(3*z, y, π) ≈ -3*z
+    @test rotate(x+y+z, z, π/4) ≈ Vec{3}((0.0,√2,1.0))
+
+    a = rand(Vec{3})
+    b = rand(Vec{3})
+    @test rotate(a, b, 0) ≈ a
+    @test rotate(a, b, π) ≈ rotate(a, b, -π)
+    @test rotate(a, b, π/2) ≈ rotate(a, -b, -π/2)
+end # of testset
+
 end
 end # of testset