Regularized forms of the traction and tangent derivative boundary integral equations of elasticity are derived for the case of closed regions. The hypersingular and strongly singular integrals of the displacement gradient representation are regularized independently, through identities of the fundamental solution and its various derivatives, before the boundary integral equations are formed. Besides the displacements and the tractions, only the tangential derivatives of the displacements evaluated at the singular point appear in the regularized equations making them well suited for numerical treatment. The regularization of the hypersingular integrals demands that the displacement components have Holder continuous first derivatives at the singular point. Consistent with this requirement, the regularization of the strongly singular integrals is effective if the tractions and the unit vectors normal and tangent to the surface are continuous at that location.;Higher order elements for two and three dimensional elastostatic problems are implemented through the coincident collocation of regularised forms of the displacement and the tangent derivative equations. The nodal values of the displacements, the fractions and their tangential derivatives are used as the degrees of freedom associated with the functional representation of the boundary variables. The tangential derivatives of the displacements and the tractions at the functional nodes are directly recovered from the boundary solution with comparable accuracy as the primative variables. Hence, the nodal values of the stress components are directly obtained through Hooke's law and need not be determined in a post processing manner. Several numerical examples demonstrate the advantages of the higher order elements versus the conventional ones. In two dimensions, four degrees of freedom per node Hermitian elements are used for functional interpolation only on those portions of the boundary where the gradients are high and quadratic Lagrangian elements are employed for the remaining parts of the modelled region. In three dimensions, nine degrees of freedom per node, incomplete quartic elements are employed for the approximation of the displacements and the tractions. Finally, the methodology presented here is general and can be extended to other problems amenable to a boundary integral formulation.