To estimate entire elastic-plastic behaviors of cracked bodies, fully plastic solutions are utilized with linear elastic solutions in the engineering approach. Some numerical algorithms such as the Selective Reduced Integration/Penalty Function (SRI/PF) method have been developed and utlized to calculate various two-dimensional fully plastic solutions. However, only a few three-dimensional solutions have been obtained because of their numerical instability caused by the interaction among crack-tip singularity, material nonlinearity and incompressibility. This paper describes a new finite element algorithm for three-dimensional fully plastic solutions. The algorithm is basically classified into the mixed formulations. By introducing an artificial viscosity term to the governing equations, static crack problems are converted into quasi-nonsteady ones, which are solved using the fractional step method. The conversion makes the algorithm stable even in the analyses of complex crack geometries though it would need a number of iterations. In the analyses, mixed interpolation tetrahedral elements are also employed from a viewpoint of high quality mesh generation for three-dimensional cracked geometries. Numerical accuracy of the present algorithm is clearly demonstrated through the analyses of the three-dimensional fully plastic solutions of center cracked plates.