A finite element procedure is developed for calculating the order and mode of singularities at 3D vertices in anisotropic linear elastic solids and composites. It is an extension of the method of Bazant and Estenssoro. although it differs from the latter in some essential aspects. It is based on a variational principle derived from the statement of virtual work on the surface of the unit sphere surrounding the singularity. The sphere is divided into six-node spherical triangles. The singularities of the spherical coordinates at the poles are avoided by coordinate transformations. The three matrices of the quadratic eigenvalue problem are explicitly evaluated and used to advantage. Real and complex eigenvalues and eigenvectors arc calculated by inverse treppen iteration. With a relatively small mesh, accuracy is to two decimal places. The method is capable of solving any 3D and any 2D vertex singularity problem with reasonable accuracy and without any assumptions with regard to the behavior of displacements in the neighborhood of line singularities. It is quite robust and stable. Several examples and applications to practical problems are given. A procedure to handle the incompressible case is discussed.