Very simple PrismaticJoint numerically explodes for certain angles

[hpmv]

Setup: Two rigidbody cubes, one at (0, 0, 0) and one at (1, 0, 0), with a prismatic joint between them with local anchors (0, 0, 0) and the global common axis of (1, 0, 0). Both objects are exempt from gravity, so there's no force other than the joint itself, and with the way the positions are set up, the joint is already satisfied.

However, if the second rigidbody already has a rotation transform, even though the local axis for the second rigidbody is specified so that the global axis is still (1, 0, 0), for some rotations, the simulation explodes even in one step.

    #[test]
    fn test_prismatic_joint() {
        for i in 0..8 {
            let gravity = Vector3::new(0.0, -9.81, 0.0);
            let mut integration_parameters = IntegrationParameters::default();
            integration_parameters.dt = 0.016;
            let mut island_manager = IslandManager::new();
            let mut broad_phase = DefaultBroadPhase::new();
            let mut narrow_phase = NarrowPhase::new();
            let mut impulse_joint_set = ImpulseJointSet::new();
            let mut multibody_joint_set = MultibodyJointSet::new();
            let mut ccd_solver = CCDSolver::new();
            let mut query_pipeline = QueryPipeline::new();
            let mut physics_pipeline = PhysicsPipeline::new();
            let mut rigid_body_set = RigidBodySet::new();
            let mut collider_set = ColliderSet::new();

            let rotation = AngVector::new(0.0, std::f32::consts::PI / 2.0 * i as f32, 0.0);
            let local_axis_2 = Rotation::new(rotation).inverse() * Vector3::x_axis();
            assert!((Rotation::new(rotation) * local_axis_2).dot(&Vector3::x_axis()) > 0.999);

            let body1 = rigid_body_set.insert(
                RigidBodyBuilder::dynamic()
                    .translation(Vector3::new(0.0, 0.0, 0.0))
                    .gravity_scale(0.0)
                    .build(),
            );
            let body2 = rigid_body_set.insert(
                RigidBodyBuilder::dynamic()
                    .translation(Vector3::new(1.0, 0.0, 0.0))
                    .rotation(rotation)
                    .gravity_scale(0.0)
                    .build(),
            );

            let collider1 = ColliderBuilder::cuboid(1.0, 1.0, 1.0);
            let collider2 = ColliderBuilder::cuboid(1.0, 1.0, 1.0);

            collider_set.insert_with_parent(collider1, body1, &mut rigid_body_set);
            collider_set.insert_with_parent(collider2, body2, &mut rigid_body_set);

            let joint = PrismaticJointBuilder::new(Vector3::x_axis())
                .local_anchor1(Point3::new(0.0, 0.0, 0.0))
                .local_anchor2(Point3::new(0.0, 0.0, 0.0))
                .local_axis1(Vector3::x_axis())
                .local_axis2(local_axis_2)
                .contacts_enabled(false)
                .build();

            impulse_joint_set.insert(body1, body2, joint, true);

            physics_pipeline.step(
                &gravity,
                &integration_parameters,
                &mut island_manager,
                &mut broad_phase,
                &mut narrow_phase,
                &mut rigid_body_set,
                &mut collider_set,
                &mut impulse_joint_set,
                &mut multibody_joint_set,
                &mut ccd_solver,
                Some(&mut query_pipeline),
                &(),
                &(),
            );

            let translation1 = rigid_body_set.get(body1).unwrap().position().translation;
            let translation2 = rigid_body_set.get(body2).unwrap().position().translation;

            println!("Case {}:", i);

            println!(
                "  Body 1: {:.4} {:.4} {:.4}",
                translation1.x, translation1.y, translation1.z
            );
            println!(
                "  Body 2: {:.4} {:.4} {:.4}",
                translation2.x, translation2.y, translation2.z
            );
        }
    }

This selects 8 different angles to specify the same global axis, but some simulations are stable, but some are not:

Case 0:
  Body 1: 0.0000 0.0000 0.0000
  Body 2: 1.0000 0.0000 0.0000
Case 1:
  Body 1: 0.0000 0.0000 0.0000
  Body 2: 1.0000 -0.0000 -0.0000
Case 2:
  Body 1: 0.0000 -45756184.0000 0.0000
  Body 2: 0.0000 45756184.0000 0.0000
Case 3:
  Body 1: 0.2545 -1.0026 0.4699
  Body 2: 0.7455 1.0026 -0.4699
Case 4:
  Body 1: 0.0000 -0.0000 0.0000
  Body 2: 1.0000 0.0000 -0.0000
Case 5:
  Body 1: 0.0000 -0.0000 0.0000
  Body 2: 1.0000 0.0000 -0.0000
Case 6:
  Body 1: 0.7500 -23258734.0000 -22446404.0000
  Body 2: -0.5000 23258734.0000 22446404.0000
Case 7:
  Body 1: 0.0000 -7340032.0000 -75863659905024.0000
  Body 2: 0.0000 7340032.0000 75863659905024.0000

[hpmv]

actually i luckily figured this out... prismatic joint's "local axis" approach is not sound, as it arbitrarily derives the other two axes based on that one axis vector. If the resulting frames end up being unaligned then the solver goes nuts. If I use a GenericJointBuilder and specify sane frames, then this is all fine.

The docs mention that there's the ability to specify a local tangent axis but that is nowhere to be found in the code. If that were implemented (or automatically calculated correctly like the document says), then we should be good.