ic to guesses

Author: baggepinnen

src/fancy_joints.jl

@@ -331,7 +331,7 @@ This joint aggregation can be used in cases where in reality a rod with spherica
         render = render, [description="Whether or not to render the joint in animations"]
     end
 
-    pars = collect_all(pars)
+    # pars = collect_all(pars)
 
     vars = @variables begin
         f_rod(t), [description="Constraint force in direction of the rod (positive on frame_a, when directed from frame_a to frame_b)"]
@@ -340,14 +340,14 @@ This joint aggregation can be used in cases where in reality a rod with spherica
         # e3_ia(t)[1:3], [description="Unit vector perpendicular to eRod_ia and e2_ia, resolved in frame_ia"]
         f_b_a1(t)[1:3], [description="frame_b.f without f_rod part, resolved in frame_a (needed for analytic loop handling)"]
         eRod_a(t)[1:3], [description="Unit vector in direction of rRod_a, resolved in frame_a (needed for analytic loop handling)"]
-        rRod_0(t)[1:3] = rRod_ia, [description="Position vector from frame_a to frame_b resolved in world frame"]
-        rRod_a(t)[1:3] = rRod_ia, [description="Position vector from frame_a to frame_b resolved in frame_a"]
-        (constraint_residue(t) = 0), [description="Constraint equation of joint in residue form: Either length constraint (= default) or equation to compute rod force (for analytic solution of loops in combination with Internal.RevoluteWithLengthConstraint/PrismaticWithLengthConstraint)"]
+        rRod_0(t)[1:3], [guess = rRod_ia, description="Position vector from frame_a to frame_b resolved in world frame"]
+        rRod_a(t)[1:3], [guess = rRod_ia, description="Position vector from frame_a to frame_b resolved in frame_a"]
+        (constraint_residue(t)), [guess = 0, description="Constraint equation of joint in residue form: Either length constraint (= default) or equation to compute rod force (for analytic solution of loops in combination with Internal.RevoluteWithLengthConstraint/PrismaticWithLengthConstraint)"]
         f_b_a(t)[1:3], [description="frame_b.f resolved in frame_a"]
         f_ia_a(t)[1:3], [description="frame_ia.f resolved in frame_a"]
         t_ia_a(t)[1:3], [description="frame_ia.t resolved in frame_a"]
         n2_a(t)[1:3], [description="Vector in direction of axis 2 of the universal joint (e2_ia), resolved in frame_a"]
-        (length2_n2_a(t) = 1), [description="Square of length of vector n2_a"]
+        (length2_n2_a(t)), [guess = 1, description="Square of length of vector n2_a"]
         length_n2_a(t), [description="Length of vector n2_a"]
         e2_a(t)[1:3], [description="Unit vector in direction of axis 2 of the universal joint (e2_ia), resolved in frame_a"]
         e3_a(t)[1:3], [description="Unit vector perpendicular to eRod_ia and e2_a, resolved in frame_a"]

test/test_fourbar.jl

@@ -11,7 +11,7 @@ D = Differential(t)
 world = Multibody.world
 
 systems = @named begin
-    j1 = Revolute(n = [1, 0, 0], w0 = 5.235987755982989, state_priority=10.0, radius=0.1f0) # Increase state priority to ensure that this joint coordinate is chosen as state variable
+    j1 = Revolute(n = [1, 0, 0], w0 = 5.235987755982989, state_priority=10, radius=0.1f0) # Increase state priority to ensure that this joint coordinate is chosen as state variable
     j2 = Prismatic(n = [1, 0, 0], s0 = -0.2)
     b1 = BodyShape(r = [0, 0.5, 0.1], radius=0.03)
     b2 = BodyShape(r = [0, 0.2, 0], radius=0.03)
@@ -38,7 +38,6 @@ connections = [connect(j2.frame_b, b2.frame_a)
                connect(b0.frame_b, j2.frame_a)
                ]
 @named fourbar2 = System(connections, t, systems = [world; systems])
-fourbar2 = complete(fourbar2)
 ssys = multibody(fourbar2)
 
 prob = ODEProblem(ssys, [], (0.0, 1.4399)) # The end time is chosen to make the animation below appear to loop forever
@@ -49,7 +48,7 @@ sol = solve(prob, FBDF(autodiff=true));
 
 
 systems = @named begin
-    j1 = Revolute(n = [1, 0, 0], w0 = 5.235987755983, state_priority=12.0, radius=0.1f0) # Increase state priority to ensure that this joint coordinate is chosen as state variable
+    j1 = Revolute(n = [1, 0, 0], w0 = 5.235987755983, phi0=nothing, state_priority=12, radius=0.1f0) # Increase state priority to ensure that this joint coordinate is chosen as state variable
     j2 = Prismatic(n = [1, 0, 0], s0 = -0.2)
     b1 = BodyShape(r = [0, 0.5, 0.1], radius=0.03)
     b2 = BodyShape(r = [0, 0.2, 0], radius=0.03)
@@ -67,9 +66,8 @@ connections = [connect(j2.frame_b, b2.frame_a)
 ]
 
 @named model = System(connections, t, systems = [world; systems])
-model = complete(model)
-ssys = multibody(model)
-prob = ODEProblem(ssys, [], (0.0, 1.4399)) # The end time is chosen to make the animation below appear to loop forever
+ssys = multibody(model, allow_symbolic=true)
+prob = ODEProblem(ssys, [], (0.0, 1.4399), guesses = [ssys.j1.phi => 0]) # The end time is chosen to make the animation below appear to loop forever
 sol2 = solve(prob, FBDF(autodiff=true)) # 3.9x faster than above
 # plot(sol2, idxs=[j2.s]) # Plot the joint coordinate of the prismatic joint (green in the animation below)
 

test/test_wheels.jl

@@ -14,6 +14,7 @@ using LinearAlgebra
                             I_long = 0.12,
                             x0 = 0.2,
                             z0 = 0.2,
+                            angles = [0,0,0],
                             der_angles = [0, -5, -1])
     end
 
@@ -31,22 +32,9 @@ end
 @named worldwheel = WheelInWorld()
 
 defs = Dict([
-
-])
-
-ssys = multibody(worldwheel)
-guesses = Dict([
-    # collect(ssys.wheel.wheeljoint.angles) .=> [0,0,0];
-    ssys.wheel.wheeljoint.v_0[2] => 0.00001;
-# (ssys.wheel.wheeljoint.f_wheel_0)[3]  => 0
-# (ssys.wheel.wheeljoint.vContact_0)[3]  => 0
-# (ssys.wheel.wheeljoint.delta_0)[2]  => 0
-# (ssys.wheel.wheeljoint.f_wheel_0)[2]  => 0
-# (ssys.wheel.wheeljoint.f_wheel_0)[1]  => 0
-# (ssys.wheel.wheeljoint.delta_0)[1]  => 0
 ])
 
-prob = ODEProblem(ssys, [], (0, 4); guesses, missing_guess_value = MissingGuessValue.Random(Random.GLOBAL_RNG))
+prob = ODEProblem(ssys, defs, (0, 4); guesses, missing_guess_value = MissingGuessValue.Random(Random.GLOBAL_RNG))
 @test prob[collect(worldwheel.wheel.wheeljoint.der_angles)] == prob[collect(worldwheel.wheel.wheeljoint.der_angles)]
 
 sol = solve(prob, Tsit5(), abstol=1e-8, reltol=1e-8)