1. That function is not aware of any joint limits (or velocity limits), it just computes the raw forward dynamics in the full state space of the robot.

For fixed joints, you should actually solve a constrained dynamics equation. In other words, solve for fixed-joint torques such that accelerations are 0. Let F be the k x n matrix that selects out the k fixed joints. Let tf be the fixed-joint torques and ta be the vector of joint torques for active joints, with 0's in the fixed joint positions. Then the overall torque is t = ta + F^{T tf, and the dynamics are B q'' + C + G = t. You would like to find q'' under the constraint that qf'' = F*q'' = 0. The forward dynamics solver gives q'' = B}-1 (t - C - G). So, tf needs to satisfy

F*B^-1 (ta + F^T tf - C - G) = 0

Expanding, F*B^-1 (ta - C - G) = - F B^-1 F^{T tf, or in other words, tf = - (F B}-1 F^T)-1 F B^-1 (ta - C - G). Replacing this back into the definition of q'', we get q'' = (I - B^-1 F^{T (F B}-1 F^T)-1 F) B^-1 (ta - C - G).

Note that the vector part B^-1 (ta - C - G) is exactly the result that you get from accelFromTorques.

(Come to think of it, it shouldn't be hard to modify the Newton-Euler (Featherstone) solver to account for fixed joints... but this should take a bit more work.)

2. The accelerations at virtual joints are the accelerations of the virtual linkage caused by the robot's link torques & gravity. (Note that the rotation accelerations need to be interpreted with respect to the Euler angles representing the current configuration, it is NOT an angular acceleration vector )

Thanks for the detailed reply, and I have one follow-up question. It seems like that the B matrix should be the one from the robot with only active joints, since the B matrices do differ when there are fixed joints. I just want to confirm this.

The B matrix is defined over all joints.