Using GLMakie to animate crank mechanism

General Usage
, ,
1

I am trying to animate a crank mechanism using the GLMakie package. I was following this youtube tutorial, Making animations and interactive applications in Makie.jl - YouTube , and here’s the code I wrote.

using GLMakie

crank = 3.3

fixed = 6.8

coupler = 2

output = 9.3

θ1 = (61/90)*pi

k1 = fixed/crank

k2 = fixed/output

k3 = ((output^2 - coupler^2 + output^2 + fixed^2)/(2 *crank* output))

P1 = [0;0]

p4x = fixed*-cos(θ1)

p4y = fixed*-sin(θ1)

function mysqrt(x)

if x < 0

    0.5

else

    sqrt(x)

end

end

θ2 = Array(0:0.05:2*pi)

function xycoordinates(θ2)

A = (1-k2)*cos(θ2) - k1 + k3

B = -2*sin(θ2)

C = k1 - (1 + k2)*cos(θ2 + k3)

b = mysqrt(B^2 - 4A*C)

θ4 = 2*(atan((-B - b)/2A, (-B + b)/2A))


p2x = crank*cos(θ2)

p2y = crank*sin(θ2)

p3x = output*-cos(θ4)

p3y = output*sin(θ4)

return p2x, p2y, p3x, p3y

end

function progress_for_one_step!(θ2)

step!(θ2)

u = θ2.u

return xycoordinates(u)

end

p2x, p3x, p2y, p3x = xycoordinates(θ2)

mechanism = Observable([Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)])

fig = Figure(); display(fig)

ax = Axis(fig[1,1])

lines!(ax, mechanism; linewidth=4, color = :purple)

ax.title = “crank mechanism”

function animstep!(θ2, mechanism)

p2x, p3x, p2y, p3x = progress_for_one_step!(θ2)

mechanism[] = [Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)]

end

for i in 1:1000

animstep!(θ2, mechanism)

sleep(0.001)

end

Almost every line of the code runs except the line “p2x, p3x, p2y, p3x = xycoordinates(θ2)” which prevents the plot from animating. I think the problem could be linked to my “progress_for_one_step” function and how my data is presented in the program. Could someone help me in identifying and correcting the error? Thank you.

2

Your function xycoordinates seems to be intended to operate on a scalar, but here

p2x, p3x, p2y, p3x = xycoordinates(θ2)

you are calling it with θ2 being an array. Therefore the first error I see is

ERROR: MethodError: no method matching cos(::Vector{Float64})
1 Like
3

The idea was to calculate p2x, p3x, p2y, p3x for all the values of θ2. Maybe a for loop would have worked but I don’t see how I can write the progress_for_one_step with it, and also the other following functions and the animations.

4

Then you probably have to adapt xycoordinates to use broadcasting? The first line should look something like

    A = (1-k2)*cos.(θ2) .- k1 .+ k3

for example.

Edit: or the simpler solution: broadcast xycoordinates like

p2x, p3x, p2y, p3y = xycoordinates.(θ2)
1 Like
5

This solves the method error problem. However when I print p2x, p3x, p2y or p3y it returns only a 4 element vector.

 p3y
(3.262944557188939, 0.49314583716287747, -8.191358441435042, -4.40359476836038)

Shouldn’t each return a 126 element vector?

It’s a bit confusing cause this p2x, p3x, p2y, p3y = xycoordinates.(θ2) returns the following,

126-element Vector{NTuple{4, Float64}}:
 (3.3, 0.0, -5.694607616035193e-16, -9.3)
 (3.2958758593033886, 0.16493125859323848, -3.575547441280739, -8.585188436787558)
 (3.283513745417485, 0.3294502749345329, -6.406050177980947, -6.741848494084563)
 (3.262944557188939, 0.49314583716287747, -8.191358441435042, -4.40359476836038)
 (3.234219706876097, 0.655608791623702, -9.059692957512048, -2.10046745168941)
 (3.1974109916451274, 0.8164330655399257, -9.299493352799107, -0.09707410182553187)
 (3.1526104141144997, 0.9752166819824206, -9.17138658677134, 1.541320237974818)
 (3.09992995239635, 1.1315627646029895, -8.854311906471532, 2.844496557023141)
 (3.0395012802095205, 1.2850805296185466, -8.455865787554524, 3.8714769510956564)
 (2.9714754377638335, 1.4353862625670597, -8.035875190729774, 4.6813149775478395)
 (2.89602245423823, 1.5821042773938698, -7.625577377449992, 5.32358616540795)
 (2.8133309227963688, 1.7248678554711754, -7.240040434941139, 5.837106689141233)
 (2.723607529201938, 1.863320162203617, -6.885502387781394, 6.251388395217235)
 ⋮
 (2.7545521899692273, -1.817262290572204, 7.001109501653403, 6.121639138813921)
 (2.8419349766330186, -1.6773209557475526, 7.3664180702933395, 5.676784724970269)
 (2.9222144059063537, -1.5331871920653966, 7.76130731451458, 5.123681173693625)
 (2.9951898210247645, -1.3852212588715342, 8.177158413287465, 4.429907480298212)
 (3.060678821455318, -1.2337929939397787, 8.594493631737825, 3.5531224597553637)
 (3.1185177188029507, -1.079280889070144, 8.974161142141433, 2.440170443817147)
 (3.1685619459462075, -0.9220711440564553, 9.242755549872747, 1.0302765867943209)
 (3.210686418379762, -0.7625567013879915, 9.271268734381575, -0.730462904519232)
 (3.2447858468605295, -0.6011362640979135, 8.852712283402793, -2.849471043420044)
 (3.2707750005759477, -0.4382132992133052, 7.7056524917506914, -5.2070067866650955)
 (3.2885889201766174, -0.2741950292977381, 5.566920703166823, -7.449791532966048)
 (3.298183080140848, -0.1094914146069375, 2.4258061282479093, -8.97805461267389)
6

I see, sorry. You probably want this

coords = xycoordinates.(θ2)
p2x = map(t -> t[1], coords)
p3x = map(t -> t[2], coords)
p2y = map(t -> t[3], coords)
p3y = map(t -> t[4], coords)
@show p2x p3x p2y p3y

then?

Edit: also found this nice solution

unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
p2x, p3x, p2y, p3y = unzip(xycoordinates.(θ2))
1 Like
7

Thanks I got the difference between the two. I thought the second method would work but all hit an error in this line
mechanism = Observable([Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)])
the error is:

ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float32
Closest candidates are:
  convert(::Type{T}, ::ColorTypes.Gray24) where T<:Real at C:\Users\ogwen\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:114
  convert(::Type{T}, ::ColorTypes.Gray) where T<:Real at C:\Users\ogwen\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:113
  convert(::Type{T}, ::Static.StaticFloat64{N}) where {N, T<:AbstractFloat} at C:\Users\ogwen\.julia\packages\Static\A5kpr\src\float.jl:26
  ...
Stacktrace:
 [1] macro expansion
   @ C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\util.jl:11 [inlined]
 [2] convert_ntuple
   @ C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\util.jl:7 [inlined]
 [3] Point{2, Float32}(x::Tuple{Vector{Float64}, Vector{Float64}})
   @ GeometryBasics C:\Users\ogwen\.julia\packages\GeometryBasics\WMp6v\src\fixed_arrays.jl:19
 [4] Point{2, Float32}(::Vector{Float64}, ::Vararg{Vector{Float64}})
   @ StaticArrays C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\convert.jl:4
 [5] top-level scope
   @ REPL[109]:1

It used to run past this point though, don’t know what changed.

8

I noticed that Point2f0 isn’t available for me. p2... and p3... are vectors here (as they should?) and you probably need some broadcasting again

mechanism = Observable([Point2f(0,0), Point2f.(p2x,p2y), Point2f.(p3x,p3y), Point2f(p4x,p4y)])

?

1 Like
9

Thanks, that got it through. There however seems to be some conversion issue with GLMakie. This line is meant to plot the lines but generates an error.

 lines!(ax, mechanism; linewidth=4, color = :purple)
ERROR: `Makie.convert_arguments` for the plot type Lines{Tuple{Vector{AbstractVector}}} and its conversion trait PointBased() was unsuccessful.

The signature that could not be converted was:
::Vector{AbstractVector}

Makie needs to convert all plot input arguments to types that can be consumed by the backends (typically Arrays with Float32 elements).
You can define a method for `Makie.convert_arguments` (a type recipe) for these types or their supertypes to make this set of arguments convertible (See http://makie.juliaplots.org/stable/documentation/recipes/index.html).

Alternatively, you can define `Makie.convert_single_argument` for single arguments which have types that are unknown to Makie but which can be converted to known types and fed back to the conversion pipeline.
10

According to this error, I think p2... and p3... should be arrays of Float32 type.

11

We can fix this with splatting

mechanism = Observable([Point2f(0,0), Point2f.(p2x,p2y)..., Point2f.(p3x,p3y)..., Point2f(p4x,p4y)])
1 Like
12

I just looked up the script for the demo: I believe we are doing it wrong…

1 Like
13

I suspected so. My interpretation must have been wrong. Especially in the implementation of the xycoordinates and progress_for_one_step functions.

14

In the first step, we’re supposed to initialize the simulation in a stepping manner, but since his problem is different from what we’re solving I don’t quite understand how he achieves that. I can’t draw the boundary between when he is just solving the diff eqn, or when he is initializing it in a stepping manner.
In our case θ2 is the known variable from which we’re analytically solving for θ4 and hence p2... and p3...

15

OK, the following program does something

using GLMakie

crank = 3.3
fixed = 6.8
coupler = 2
output = 9.3
θ1 = (61/90)*pi
k1 = fixed/crank
k2 = fixed/output
k3 = ((output^2 - coupler^2 + output^2 + fixed^2)/(2 *crank* output))
P1 = [0;0]
p4x = fixed*-cos(θ1)
p4y = fixed*-sin(θ1)

function mysqrt(x)
    if x < 0
        0.5
    else
        sqrt(x)
    end
end

θ2 = Array(0:0.05:2*pi)

function xycoordinates(θ2)
    A = (1-k2)*cos(θ2) - k1 + k3
    B = -2*sin(θ2)
    C = k1 - (1 + k2)*cos(θ2 + k3)
    b = mysqrt(B^2 - 4A*C)
    θ4 = 2*(atan((-B - b)/2A, (-B + b)/2A))
    p2x = crank*cos(θ2)
    p2y = crank*sin(θ2)
    p3x = output*-cos(θ4)
    p3y = output*sin(θ4)
    return p2x, p2y, p3x, p3y
end

fig = Figure(); display(fig)
ax = Axis(fig[1,1])
ax.title = "crank mechanism"

p2x, p3x, p2y, p3y = xycoordinates(θ2[1])
mechanism = Observable([Point2f(0,0), Point2f(p2x,p2y), Point2f(p3x,p3y), Point2f(p4x,p4y)])
lines!(ax, mechanism; linewidth=4, color = :purple)
for i in 2:length(θ2)
    global p2x, p3x, p2y, p3x = xycoordinates(θ2[i])
    mechanism[] = [Point2f(0,0), Point2f(p2x,p2y), Point2f(p3x,p3y), Point2f(p4x,p4y)]
    sleep(0.001)
end

Edit: fixed some mistakes…

1 Like
16

Yes!!!
This is exactly what I wanted. Thank you so much, I’m deeply grateful for your help. Could you mind explaining the logic behind how you did it?

17

Unfortunately there was not much logic to it, only experience and experimentation…

2 Likes
18

@goerch

I am trying the code but it only show axis. I was hoping to see real crank there…

Capture d’écran_2022-07-15_15-20-44
Capture d’écran_2022-07-15_15-20-44794×612 9.62 KB

Maybe the code or the package installed in my Julia is the one at fault?

19

Hi, I just wanted to add that this approach can also be used to record videos of moving objects using GLMakie. I have been searching for a way to do that for several hours. And I made the following adaptation of your code to record a video of a moving circle. I did not find any other way to do this with GLMakie, as the old plots would always stick around. Another advantage of this approach is that it is possible to have parts of the plot that are not moving only be plotted once, and then have a moving part on top.

using GLMakie

# GLMakie
fig = Figure(resolution = (1920, 1080))
ax1 = fig[1, 1] = Axis(fig,
                        aspect=1.0,
                        limits = (-2, 2, -2, 2))

# What to plot
frames = 100 # frames
t_range = range(0.0, 10.0, length=frames) # time vector
radius = sin.(t_range) .+ 1.0 # changing radius
theta = range(0.0 ,2 * pi ,length=150) # angle

# circle function
function circleFun(theta,radius)
    # define circle (parametric)
    a = radius * cos.( theta )
    b = radius * sin.( theta )
    return a, b
end

# set everything up
a, b = circleFun(theta,radius[1]) # call function
mechanism = Observable([Point2f(a[i],b[i]) for i in eachindex(a)]) # make mechanism
lines!(ax1, mechanism, linewidth=5, color=:red) # plot initial
frame_iterator = 1:1:frames # frame iterator

# record video of moving circle
record(fig, "movingCircle.mp4", frame_iterator; framerate = 12) do i
    global a, b = circleFun(theta,radius[i])
    mechanism[] = [Point2f(a[j],b[j]) for j in eachindex(a)]
end