In this investigation, a numerical procedure for modeling sliding and nonsliding joint constraints for the B-spline thin plate element is developed for the large deformation analysis of multibody systems. A concept of intermediate reference coordinates proposed for the absolute nodal coordinate formulation is generalized for B-spline elements such that a wide variety of joint constraints can be modeled using existing joint constraint libraries already implemented in multibody dynamics codes. This procedure allows for modeling sliding joints for B-spline elements that requires a solution to moving boundary problems by introducing time-variant surface parameters in the B-spline parametric domain. Since surface parameters treated as knot variables in the basis function are defined in the entire parametric domain rather than the element domain, the location of the constraint definition point can be determined without knowing in which elements the sliding point is located. Furthermore, using the B-spline recurrence formula, control points used for describing the constraint equations can be systematically extracted. It is shown that many types of nonsliding joints fixed on the flexible body can also be modeled as a special case of the sliding joint formulation developed in this investigation, leading to a unified joint constraint formulation for B-spline elements. Several numerical examples are presented in order to demonstrate the use of the numerical procedure developed in this investigation.