From 961099d6a61aafc6c8c0e55441a3e8e270cadaef Mon Sep 17 00:00:00 2001 From: Taylor McDonnell Date: Tue, 17 May 2022 14:41:29 -0600 Subject: [PATCH 1/2] fix rotation parameter scaling jacobians --- src/element.jl | 6 +- src/math.jl | 208 ++++++++++++++++++++++++----------------------- test/runtests.jl | 4 +- 3 files changed, 113 insertions(+), 105 deletions(-) diff --git a/src/element.jl b/src/element.jl index fa37b1418..96b9f01e0 100644 --- a/src/element.jl +++ b/src/element.jl @@ -2300,7 +2300,7 @@ time-marching simulation. Ct_θ1, Ct_θ2, Ct_θ3 = C_θ1', C_θ2', C_θ3' Cdot_θ1, Cdot_θ2, Cdot_θ3 = get_C_t_θ(θ, θdot) Ctdot_θ1, Ctdot_θ2, Ctdot_θ3 = Cdot_θ1', Cdot_θ2', Cdot_θ3' - Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(C, θ) + Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(θ) Ctdot_θdot1, Ctdot_θdot2, Ctdot_θdot3 = Cdot_θdot1', Cdot_θdot2', Cdot_θdot3' # element strain and curvature @@ -2721,7 +2721,7 @@ residual equations with respect to the time derivatives of the state variables CtCab = Ct*Cab # rotation matrix derivatives - Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(C, θ) + Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(θ) Ctdot_θdot1, Ctdot_θdot2, Ctdot_θdot3 = Cdot_θdot1', Cdot_θdot2', Cdot_θdot3' # element linear and angular momentum @@ -2782,7 +2782,7 @@ by the scaling parameter `gamma`. CtCab = Ct*Cab # rotation matrix derivatives - Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(C, θ) + Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 = get_C_t_θdot(θ) Ctdot_θdot1, Ctdot_θdot2, Ctdot_θdot3 = Cdot_θdot1', Cdot_θdot2', Cdot_θdot3' # element linear and angular momentum diff --git a/src/math.jl b/src/math.jl index 10985b8eb..bde336f85 100644 --- a/src/math.jl +++ b/src/math.jl @@ -55,21 +55,37 @@ parameter allows deflections greater than 360 degrees. θnorm = sqrt(θ'*θ) - if θnorm <= 4 - scaling = one(θnorm) - else - half_rotations = θnorm/4 - irot, m = floor(half_rotations/2), half_rotations%2 - if m <= 1 - scaling = 1 .- 8*irot/θnorm - else - scaling = 1 .- 8*(irot+1)/θnorm - end + scaling = 1 + + m = div(div(θnorm, 4) + 1, 2) # number of times to add/subtract + + if θnorm > 4 + scaling -= 8*m/θnorm end return scaling end +@inline function rotation_parameter_scaling_θ(θ) + + return rotation_parameter_scaling_θ(rotation_parameter_scaling(θ), θ) +end + +@inline function rotation_parameter_scaling_θ(scaling, θ) + + return ifelse(isone(scaling), zero(θ), (1 - scaling)*θ/(θ'*θ)) +end + +@inline function rotation_parameter_scaling_θ_θ(θ) + + return rotation_parameter_scaling_θ_θ(rotation_parameter_scaling(θ), θ) +end + +@inline function rotation_parameter_scaling_θ_θ(scaling, θ) + + return ifelse(isone(scaling), (@SMatrix zeros(eltype(θ),3,3)), (1 - scaling)/(θ'*θ)*(I - 3*θ*θ'/(θ'*θ))) +end + """ get_C(θ) @@ -84,6 +100,56 @@ Returns the transformation matrix `C` given the three angular parameters in `θ` return wiener_milenkovic(c) end +""" + get_C_θ([C, ] θ) + +Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with +respect to each of the rotation parameters in `θ`. +""" +@inline get_C_θ(c) = get_C_θ(get_C(c), c) + +@inline function get_C_θ(C, θ) + + scaling = rotation_parameter_scaling(θ) + scaling_θ = rotation_parameter_scaling_θ(scaling, θ) + + c = scaling*θ + c_θ = scaling_θ*θ' + scaling*I3 + + c0 = 2 - c'*c/8 + c0_c = -c/4 + + tmp = 1/(4-c0)^2 + tmp_c = 2/(4-c0)^3*c0_c + + C_c1 = tmp_c[1]*C/tmp + tmp*(@SMatrix [ + 2*c0*c0_c[1] + 2*c[1] 2*(c[2] + c0_c[1]*c[3]) 2*(c[3] - c0_c[1]*c[2]); + 2*(c[2] - c0_c[1]*c[3]) 2*c0*c0_c[1] - 2*c[1] 2*(c0_c[1]*c[1] + c0); + 2*(c[3] + c0_c[1]*c[2]) -2*(c0_c[1]*c[1] + c0) 2*c0*c0_c[1] - 2*c[1] + ] + ) + + C_c2 = tmp_c[2]*C/tmp + tmp*(@SMatrix [ + 2*c0*c0_c[2] - 2*c[2] 2*(c[1] + c0_c[2]*c[3]) -2*(c0_c[2]*c[2] + c0); + 2*(c[1] - c0_c[2]*c[3]) 2*c0*c0_c[2] + 2*c[2] 2*(c[3] + c0_c[2]*c[1]); + 2*(c0_c[2]*c[2] + c0) 2*(c[3] - c0_c[2]*c[1]) 2*c0*c0_c[2] - 2*c[2] + ] + ) + + C_c3 = tmp_c[3]*C/tmp + tmp*(@SMatrix [ + 2*c0*c0_c[3] - 2*c[3] 2*(c0_c[3]*c[3] + c0) 2*(c[1] - c0_c[3]*c[2]); + -2*(c0_c[3]*c[3] + c0) 2*c0*c0_c[3] - 2*c[3] 2*(c[2] + c0_c[3]*c[1]); + 2*(c[1] + c0_c[3]*c[2]) 2*(c[2] - c0_c[3]*c[1]) 2*c0*c0_c[3] + 2*c[3] + ] + ) + + C_θ1 = C_c1*c_θ[1,1] + C_c2*c_θ[2,1] + C_c3*c_θ[3,1] + C_θ2 = C_c1*c_θ[1,2] + C_c2*c_θ[2,2] + C_c3*c_θ[3,2] + C_θ3 = C_c1*c_θ[1,3] + C_c2*c_θ[2,3] + C_c3*c_θ[3,3] + + return C_θ1, C_θ2, C_θ3 +end + """ get_C_t([C, ] θ, θ_t) @@ -96,9 +162,10 @@ respect to the scalar parameter `t`. `θ_t` is the derivative of the angular par @inline function get_C_t(C, θ, θ_t) scaling = rotation_parameter_scaling(θ) + scaling_θ = rotation_parameter_scaling_θ(scaling, θ) c = scaling*θ - c_t = scaling*θ_t + c_t = scaling_θ'*θ_t*θ + scaling*θ_t c0 = 2 - c'*c/8 c0dot = -c'*c_t/4 @@ -120,50 +187,6 @@ respect to the scalar parameter `t`. `θ_t` is the derivative of the angular par return Cdot end -""" - get_C_θ([C, ] θ) - -Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with -respect to each of the rotation parameters in `θ`. -""" -@inline get_C_θ(c) = get_C_θ(get_C(c), c) - -@inline function get_C_θ(C, θ) - - scaling = rotation_parameter_scaling(θ) - - c = scaling*θ - - c0 = 2 - c'*c/8 - c0_θ = -c/4 - - tmp = 1/(4-c0)^2 - tmp_θ = 2/(4-c0)^3*c0_θ - - C_θ1 = tmp_θ[1]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0_θ[1] + 2*c[1] 2*(c[2] + c0_θ[1]*c[3]) 2*(c[3] - c0_θ[1]*c[2]); - 2*(c[2] - c0_θ[1]*c[3]) 2*c0*c0_θ[1] - 2*c[1] 2*(c0_θ[1]*c[1] + c0); - 2*(c[3] + c0_θ[1]*c[2]) -2*(c0_θ[1]*c[1] + c0) 2*c0*c0_θ[1] - 2*c[1] - ] - ) - - C_θ2 = tmp_θ[2]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0_θ[2] - 2*c[2] 2*(c[1] + c0_θ[2]*c[3]) -2*(c0_θ[2]*c[2] + c0); - 2*(c[1] - c0_θ[2]*c[3]) 2*c0*c0_θ[2] + 2*c[2] 2*(c[3] + c0_θ[2]*c[1]); - 2*(c0_θ[2]*c[2] + c0) 2*(c[3] - c0_θ[2]*c[1]) 2*c0*c0_θ[2] - 2*c[2] - ] - ) - - C_θ3 = tmp_θ[3]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0_θ[3] - 2*c[3] 2*(c0_θ[3]*c[3] + c0) 2*(c[1] - c0_θ[3]*c[2]); - -2*(c0_θ[3]*c[3] + c0) 2*c0*c0_θ[3] - 2*c[3] 2*(c[2] + c0_θ[3]*c[1]); - 2*(c[1] + c0_θ[3]*c[2]) 2*(c[2] - c0_θ[3]*c[1]) 2*c0*c0_θ[3] + 2*c[3] - ] - ) - - return scaling*C_θ1, scaling*C_θ2, scaling*C_θ3 -end - """ get_C_t_θ(θ, θ_t) @@ -181,45 +204,18 @@ transformation matrix `C` with respect to `θ`. end """ - get_C_t_θdot([C, ] θ) + get_C_t_θdot(θ) Calculate the derivative of the time derivative of the Wiener-Milenkovic transformation matrix `C` with respect to each of the time derivatives of `θ`. """ -get_C_t_θdot +@inline function get_C_t_θdot(θ) -@inline get_C_t_θdot(θ) = get_C_t_θdot(get_C(θ), θ) - -@inline function get_C_t_θdot(C, θ) - - scaling = rotation_parameter_scaling(θ) - - c = scaling*θ - - c0 = 2 - c'*c/8 - c0dot_θdot = -c/4 - tmp = 1/(4-c0)^2 - tmpdot_θdot = 2/(4-c0)^3*c0dot_θdot - - C_θdot1 = tmpdot_θdot[1]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0dot_θdot[1] + 2*c[1] 2*(c[2] + c[3]*c0dot_θdot[1]) 2*(c[3] - c[2]*c0dot_θdot[1]); - 2*(c[2] - c[3]*c0dot_θdot[1]) 2*c0*c0dot_θdot[1] - 2*c[1] 2*(c[1]*c0dot_θdot[1] + c0); - 2*(c[3] + c[2]*c0dot_θdot[1]) 2*(-c[1]*c0dot_θdot[1] - c0) 2*c0*c0dot_θdot[1] - 2*c[1]] - ) - - C_θdot2 = tmpdot_θdot[2]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0dot_θdot[2] - 2*c[2] 2*(c[1] + c[3]*c0dot_θdot[2]) 2*(-c[2]*c0dot_θdot[2] - c0); - 2*(c[1] - c[3]*c0dot_θdot[2]) 2*c0*c0dot_θdot[2] + 2*c[2] 2*(c[3] + c[1]*c0dot_θdot[2]); - 2*(c[2]*c0dot_θdot[2] + c0) 2*(c[3] - c[1]*c0dot_θdot[2]) 2*c0*c0dot_θdot[2]- 2*c[2]] - ) - - C_θdot3 = tmpdot_θdot[3]*C/tmp + tmp*(@SMatrix [ - 2*c0*c0dot_θdot[3] - 2*c[3] 2*(c[3]*c0dot_θdot[3] + c0) 2*(c[1] - c[2]*c0dot_θdot[3]); - 2*(-c[3]*c0dot_θdot[3] - c0) 2*c0*c0dot_θdot[3] - 2*c[3] 2*(c[2] + c[1]*c0dot_θdot[3]); - 2*(c[1] + c[2]*c0dot_θdot[3]) 2*(c[2] - c[1]*c0dot_θdot[3]) 2*c0*c0dot_θdot[3] + 2*c[3]] - ) - - return scaling*C_θdot1, scaling*C_θdot2, scaling*C_θdot3 + Cdot_θdot1 = ForwardDiff.derivative(cdot1 -> GXBeam.get_C_t(θ, SVector(cdot1, 0, 0)), 0) + Cdot_θdot2 = ForwardDiff.derivative(cdot2 -> GXBeam.get_C_t(θ, SVector(0, cdot2, 0)), 0) + Cdot_θdot3 = ForwardDiff.derivative(cdot3 -> GXBeam.get_C_t(θ, SVector(0, 0, cdot3)), 0) + + return Cdot_θdot1, Cdot_θdot2, Cdot_θdot3 end """ @@ -245,29 +241,35 @@ end get_Q_θ(Q, θ) Calculate the derivative of the matrix `Q` with respect to each of the rotation -parameters in `c`. +parameters in `θ`. """ @inline get_Q_θ(θ) = get_Q_θ(get_Q(θ), θ) @inline function get_Q_θ(Q, θ) scaling = rotation_parameter_scaling(θ) + scaling_θ = rotation_parameter_scaling_θ(scaling, θ) c = scaling*θ + c_θ = scaling_θ*θ' + scaling*I3 c0 = 2 - c'*c/8 - c0_θ = -c/4 + c0_c = -c/4 tmp = 1/(4-c0)^2 - tmp_θ = 2/(4-c0)^3*c0_θ + tmp_c = 2/(4-c0)^3*c0_c - Q_θ1 = tmp_θ[1]*Q/tmp + tmp*(-c[1]/2*I - 2*tilde(e1) + (e1*c' + c*e1')/2) + Q_c1 = tmp_c[1]*Q/tmp + tmp*(-c[1]/2*I - 2*tilde(e1) + (e1*c' + c*e1')/2) - Q_θ2 = tmp_θ[2]*Q/tmp + tmp*(-c[2]/2*I - 2*tilde(e2) + (e2*c' + c*e2')/2) + Q_c2 = tmp_c[2]*Q/tmp + tmp*(-c[2]/2*I - 2*tilde(e2) + (e2*c' + c*e2')/2) - Q_θ3 = tmp_θ[3]*Q/tmp + tmp*(-c[3]/2*I - 2*tilde(e3) + (e3*c' + c*e3')/2) + Q_c3 = tmp_c[3]*Q/tmp + tmp*(-c[3]/2*I - 2*tilde(e3) + (e3*c' + c*e3')/2) - return scaling*Q_θ1, scaling*Q_θ2, scaling*Q_θ3 + Q_θ1 = Q_c1*c_θ[1,1] + Q_c2*c_θ[2,1] + Q_c3*c_θ[3,1] + Q_θ2 = Q_c1*c_θ[1,2] + Q_c2*c_θ[2,2] + Q_c3*c_θ[3,2] + Q_θ3 = Q_c1*c_θ[1,3] + Q_c2*c_θ[2,3] + Q_c3*c_θ[3,3] + + return Q_θ1, Q_θ2, Q_θ3 end """ @@ -297,14 +299,20 @@ rotation parameters in `θ`. @inline function get_Qinv_θ(θ) scaling = rotation_parameter_scaling(θ) + scaling_θ = rotation_parameter_scaling_θ(scaling, θ) c = scaling*θ + c_θ = scaling_θ*θ' + scaling*I3 + + Qinv_c1 = -c[1]/8*I3 + 1/2*tilde(e1) + 1/8*(e1*c' + c*e1') + Qinv_c2 = -c[2]/8*I3 + 1/2*tilde(e2) + 1/8*(e2*c' + c*e2') + Qinv_c3 = -c[3]/8*I3 + 1/2*tilde(e3) + 1/8*(e3*c' + c*e3') - Qinv_θ1 = -c[1]/8*I3 + 1/2*tilde(e1) + 1/8*(e1*c' + c*e1') - Qinv_θ2 = -c[2]/8*I3 + 1/2*tilde(e2) + 1/8*(e2*c' + c*e2') - Qinv_θ3 = -c[3]/8*I3 + 1/2*tilde(e3) + 1/8*(e3*c' + c*e3') + Qinv_θ1 = Qinv_c1*c_θ[1,1] + Qinv_c2*c_θ[2,1] + Qinv_c3*c_θ[3,1] + Qinv_θ2 = Qinv_c1*c_θ[1,2] + Qinv_c2*c_θ[2,2] + Qinv_c3*c_θ[3,2] + Qinv_θ3 = Qinv_c1*c_θ[1,3] + Qinv_c2*c_θ[2,3] + Qinv_c3*c_θ[3,3] - return scaling*Qinv_θ1, scaling*Qinv_θ2, scaling*Qinv_θ3 + return Qinv_θ1, Qinv_θ2, Qinv_θ3 end """ diff --git a/test/runtests.jl b/test/runtests.jl index 3423f8918..6df9f734b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,8 +10,8 @@ const RNG = MersenneTwister(1234) @testset "Math" begin - c = rand(RNG, 3) - cdot = rand(RNG, 3) + c = 1e3*rand(RNG, 3) + cdot = 1e3*rand(RNG, 3) # get_C_θ C_θ1, C_θ2, C_θ3 = GXBeam.get_C_θ(c) From 5cf824cf496bd20b5878a5b7a7596729e5ce7506 Mon Sep 17 00:00:00 2001 From: Taylor McDonnell Date: Tue, 17 May 2022 14:42:03 -0600 Subject: [PATCH 2/2] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a99c6fd6e..e442e1167 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GXBeam" uuid = "974624c9-1acb-4ad6-a627-8ac40fc27a3e" authors = ["Taylor McDonnell and contributors"] -version = "0.3.2" +version = "0.3.3" [deps] ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"