Singular mode shapes around the sharp frictional contact interfaces of different materials with linear elastic behavior are derived semi-analytically using the asymptotic eigen-function expansion method. The characteristic equation which is not proposed in a short explicit expression in the literature is extracted in symbolic terms; then solved numerically and the singularity order is derived using Newton's iterative method for the first time in the literature.The extracted eigen-shapes of the high gradient stress fields and the corresponding displacement fields satisfy the compatibility of deformations and tractions on the contact interface with arbitrary frictional interface equations. These asymptotic mode shapes are then implemented into the finite element interpolation using Partition of Unity Finite Element Method (PUFEM). Since the additional enrichment degrees of freedom on the singular contact nodes have the physical meaning of generalized stress intensity factor (GSIF), the direct GSIFs obtained are more accurate and the contour integral manipulations need not be maintained in the post processing phase of the contact analysis.Finally, some engineering problems are solved which demonstrate the capacities of the approach concerning the accuracy of GSIF obtained from the direct-PUFEM proposed compared with that derived by reciprocal work contour integral method used in traditional finite element analyses.