Start PlanarMechanics module

[AbdulrhmnGhanem]

Inspired by https://github.com/dzimmer/PlanarMechanics

I'll add all the basic components from the PlanarMechanics Modelica library. Are you OK with adding it to the standard library?

[codecov]

Codecov Report

Attention: Patch coverage is 0.29412% with 339 lines in your changes are missing coverage. Please review.

Project coverage is 12.98%. Comparing base (35a41f7) to head (8838212).
Report is 5 commits behind head on main.

❗ Current head 8838212 differs from pull request most recent head 611a124. Consider uploading reports for the commit 611a124 to get more accurate results

Files	Patch %	Lines
src/Mechanical/PlanarMechanics/sensors.jl	0.00%	263 Missing ⚠️
src/Mechanical/PlanarMechanics/joints.jl	0.00%	45 Missing ⚠️
src/Mechanical/PlanarMechanics/components.jl	0.00%	31 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main     #220       +/-   ##
===========================================
- Coverage   54.12%   12.98%   -41.15%     
===========================================
  Files          48       52        +4     
  Lines        1648     1987      +339     
===========================================
- Hits          892      258      -634     
- Misses        756     1729      +973     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ChrisRackauckas]

I think this idea is great! Down the line we may move this out to a separate SciML/PlanarMechanics to reduce the loading time of the standard library, but we can do any reformatting like that down the line so I wouldn't worry about that right now. For now, just getting more modules is good and we'd be happy to help maintain them.

[YingboMa]

There is @test_nowarn for testing instantiation.

[AbdulrhmnGhanem]

Thanks, I didn't know it 😅 pretty neat!

I will replace it with a proper test case like the Pendulum and Free body, it's only temporary for now.

[AbdulrhmnGhanem]

Is there a way I can do this in MTK, i.e., asserting the number of connections?

assert(cardinality(frame_a) > 0, "Connector frame_a must be connected at least once");

[baggepinnen]

Is there a way I can do this in MTK, i.e., asserting the number of connections?

Not at the moment, but it would be a nice feature to have

[AbdulrhmnGhanem]

Not at the moment, but it would be a nice feature to have

Could you hint where I can start to add this in MTK?

[baggepinnen]

Could you hint where I can start to add this in MTK?

I'm not sure TBH, but this feature is a debugging feature, the PR here should be good to go without it.

Do you feel that the PR has reached a state where it's ready for review?

[AbdulrhmnGhanem]

Do you feel that the PR has reached a state where it's ready for review?

Sorrowfully no, but I will resume working on it and get it ready as soon as possible.

[AbdulrhmnGhanem]

@baggepinnen, @YingboMa, could you take a look at the pendulum model below?

ModelingToolkitStandardLibrary.jl/test/Mechanical/planar_mechanics.jl

Lines 35 to 61 in 7051dbe

	@testset "Pendulum" begin
	# https://github.com/dzimmer/PlanarMechanics/blob/743462f58858a808202be93b708391461cbe2523/PlanarMechanics/Examples/Pendulum.mo
	@named ceiling = Fixed()
	@named rod = FixedTranslation(rx = 1.0, ry = 0.0)
	@named body = Body(m = 1, I = 0.1)
	@named revolute = Revolute()
	connections = [
	connect(ceiling.frame, revolute.frame_a),
	connect(revolute.frame_b, rod.frame_a),
	connect(rod.frame_b, body.frame),
	]
	@named model = ODESystem(connections,
	t,
	[],
	[],
	systems = [body, revolute, rod, ceiling])
	sys = structural_simplify(model)
	@test_broken length(states(sys)) == 2
	unset_vars = setdiff(states(sys), keys(ModelingToolkit.defaults(sys)))
	prob = ODEProblem(sys, unset_vars .=> 0.0, tspan, []; jac = true)
	sol = solve(prob, Rodas5P())
	@test_broken SciMLBase.successful_retcode(sol)
	end

I have been stuck with it for quite long time so far.

Why would such a system have 7 states after simplification? it's unlikely to be a problem with structural_simplify but I ran out of debugging tricks.

[baggepinnen]

Why would such a system have 7 states after simplification?

What are these 7 variables?

These equations in Revolute look suspicious

phi ~ ifelse(constant_phi === nothing, phi, constant_phi)

You are adding the equation phi ~ phi here, I'm not sure if that trips MTK up?

The constant_phi argument is not documented, so I'm not sure what it's supposed to accomplish. If you keep the angle constant, it's not really a joint any more?

[AbdulrhmnGhanem]

What are these 7 variables?

julia> states(sys)
7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 rod₊frame_b₊y(t)
 rod₊frame_b₊yˍt(t)
 rod₊frame_a₊phi(t)
 rod₊frame_b₊fx(t)
 rod₊frame_a₊phiˍt(t)
 rod₊frame_a₊phiˍtt(t)
 body₊ryˍtt(t)

julia> equations(sys)
7-element Vector{Equation}:
 Differential(t)(rod₊frame_b₊y(t)) ~ rod₊frame_b₊yˍt(t)
 Differential(t)(rod₊frame_b₊yˍt(t)) ~ body₊ryˍtt(t)
 0 ~ rod₊frame_b₊y(t) - rod₊frame_a₊y(t) - rod₊rx*sin(rod₊frame_a₊phi(t)) - rod₊ry*cos(rod₊frame_a₊phi(t))
 0 ~ rod₊frame_b₊j(t) + body₊frame₊j(t)
 0 ~ rod₊frame_b₊xˍtt(t) - rod₊frame_a₊xˍtt(t) + rod₊rx*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍtt(t) + rod₊ry*rod₊frame_a₊phiˍtt(t)*cos(rod₊frame_a₊phi(t)) + rod₊rx*cos(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2) - rod₊ry*sin(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2)
 0 ~ rod₊frame_b₊yˍt(t) - rod₊frame_a₊yˍt(t) - rod₊rx*cos(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍt(t) + rod₊ry*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍt(t)
 0 ~ -rod₊frame_a₊yˍtt(t) + body₊ryˍtt(t) - rod₊rx*rod₊frame_a₊phiˍtt(t)*cos(rod₊frame_a₊phi(t)) + rod₊ry*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍtt(t) + rod₊rx*sin(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2) + rod₊ry*cos(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2)

You are adding the equation phi ~ phi here, I'm not sure if that trips MTK up?

Removing it didn't make a difference.

The constant_phi argument is not documented, so I'm not sure what it's supposed to accomplish. If you keep the angle constant, it's not really a joint any more?

You're right, I needed constant_ω to make a revolute move with a constant angular velocity and added constant_phi, and constat_tau just in case. Is there a better way to do this?

ModelingToolkitStandardLibrary.jl/test/Mechanical/planar_mechanics.jl

Line 79 in abc9d7f

@named revolute = Revolute(constant_ω = ω)

[baggepinnen]

It looks like MTK has chosen rod variables as state

 rod₊frame_a₊phi(t)
 rod₊frame_a₊phiˍt(t)

these are equivalent to the revolute joint angle and angular velocity and should be sufficient as state variables.

You're right, I needed constant_ω to make a revolute move with a constant angular velocity and added constant_phi, and constat_tau just in case. Is there a better way to do this?

Something like this component?

https://github.com/modelica/ModelicaStandardLibrary/blob/master/Modelica/Mechanics/Rotational/Sources/ConstantSpeed.mo

i.e., a constant-speed source

[baggepinnen]

The problem here might be the same problem that we have with the multibody library. Multibody requires JuliaSimCompiler to simpify properly. I tried the pendulum example with JuliaSimCompiler but ran into a bug, so I cannot verify that this is indeed the case

julia> sys = structural_simplify(JuliaSimCompiler.IRSystem(model))
ERROR: BoundsError: attempt to access 61-element Vector{Vector{Int64}} at index [62]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] 𝑑neighbors

[AbdulrhmnGhanem]

i.e., a constant-speed source

Yes, this is better, I'll add it.

[AbdulrhmnGhanem]

The problem here might be the same problem that we have with the multibody library.

I merged the main branch into this branch, just in case this bug was fixed in the latest version.