Abstract: Research into numerical computing problems involving coupled elastoplastic damage of complex materials not only needs to select the appropriate constitutive model to forecast damage, but also needs a robust and effective constitutive integration algorithm. Firstly, the basic steps of setting up the elastoplastic damage constitutive model are expounded within the framework of thermodynamics and continuum mechanics. Secondly, a finite element solving program of the constitutive model, which is based on Lemaitre's isotropic hardening coupled elastoplastic damage constitutive model, the corresponding constitutive integration algorithm-fully implicit return mapping algorithm, and the consistent tangent modulus, is compiled using C ++ language in Visual 6.0 environment. When solving the full return mapping algorithm, it is possible to develop a simplified form which can be solved by iterating a non-linear equation, so improving computational efficiency is relatively simple. Finally, the validity of the application program is demonstrated through a numerical example of a cylindrical notched bar, and the calculation results are graphically displayed by a self-developed interface program. The research results show the effectiveness of the algorithm and the correctness of program, proving that the development process of damage of the ductile material can be simulated. Boundary value problems of finite elements can be solved well, which provides a basis for research and design for structures which considers the damage characteristics of ductile materials.