Slender space elements can be employed for the design of advanced structures, like deployable and morphing systems. Their analysis requires reliable and effective models, able to accurately reproduce the complex non-linear behaviour of such systems, often very sensitive to initial conditions, onset of instability etc. The contribution extends the invariant formulation for space slender beams based on an extended B´ezier interpolation presented in [1] to account for complex internal constraints like pivots, Hooke joints and others that are usually met in the structures object of the investigation. A G1-conforming formulation is adopted in order to deal with the exact Kirchhoff beam model based on the smallest rotation map proposed in [2]. The kinematic descriptors for the Kirchhoff beam model are the placement of the centroid curve and the rotation angle for the directors of the cross section. The centroid curve is interpolated with B´ezier polynomials and a spherical B´ezier interpolation is used for the rotation angle. In this way the spherical linear interpolation proposed in [3] is extended to more than two control variables for the rotation. In order to avoid locking phenomena and reduce the computational effort a two field mixed formulation is considered introducing discontinuous assumed internal forces, so that they can be condensed at the element level. A symmetric formulation is recovered consistently performing the variations of the relevant variables on the configuration manifold using the proper Levi-Civita connection. A path-following algorithm, generalized to the case that the configuration space is a manifold, is used. Examples are presented in order to show the reliability of the proposed finite element formulation.