A finite element computer program EPAS-J1 was developed to calculate the stress intensity factors of three-dimensional cracks. In this program, the stress intensity factor is determined by the virtual crack extension method together with the distorted elements allocated along the crack front. The program also includes the connection elements based on the Lagrange multiplier concept to connect such different kinds of elements as the solid and shell elements, or the shell and beam elements. For the structure including the three-dimensional surface cracks, the solid elements are employed only at the neighborhood of the surface crack, while the remainder of the structure is modeled by the shell or beam elements since the crack singularity is very local. Computer storage and computational time can be highly reduced with the application of the above modeling for calculation of the stress intensity factors of the three-dimensional surface cracks, because the three-dimensional solid elements are required only around the crack front. Several numerical analyses were performed by the EPAS-J1 program. At first, the accuracies of the connection element and the virtual crack extension method were confirmed using the simple structures. Compared with other techniques of connecting different kinds of elements such as the tying method or the method using an anisotropic plate element, it is found that the present connection element provides better results than the others. It is also found that the virtual crack extension method provided the accurate stress intensity factor. Furthermore, the results are also presented for the stress intensity factor analyses of cylinders with longitudinal or circumferential surface cracks using the combination of various kinds of elements together with the connection elements.