In the last week, I finalized the locate method for the point class, as well as the .rel() and .vel() methods. Just as relative angular velocity between adjacent ReferenceFrame instances is stored, so is the relative position and relative velocity of points. By doing it in this fashion, the partial velocities can be formed incrementally, and there are some nice efficiencies to be realized when it comes time to form the equations of motion by avoiding the need to take redundant dot products. This will also translate into natural definitions for intermediate variables in the final form of the equations of motion.
One subtlety that is encountered in the kinematic chain of points is the possibility for closed kinematic chains. Typical examples are four and six bar linkages, or a bicycle with both wheels on the ground. These closed kinematic loops introduce holonomic constraint equations which are generally not solvable directly, but can be dealt with in differentiated form. PyDy now has a function called kinematic_chain() that deals with such loops and automatically forms the constraint equations and the differentiated holonomic constraint equations, and its use is illustrated in examples/fourbarlinkage.py in the git repository.
I have neared the completion of my updated sin,cos,tan,cot,sec,csc functions that handle arbitrary rational shifts by pi, and am in the process of testing and making things work with the hyperbolic functions that already exist in Sympy. This has been very challenging to get all the behavior that Mathematica does, but I think I have finally converged upon a solution that will be functional. Once this behavior is done, my next task will be to implement the Fu et al. paper so that trigsimp works. I will be putting PyDy development on hold until this is complete because PyDy needs to have some good trig simplification routines for many things to simplify as they should.