The sintering of thin films is widely used for surface coatings and because of its technological importance has generated extensive research interest. During the sintering process, the thin film is constrained by the substrate, which generates considerably high levels of stresses. Crackings are often observed and are considered as one of the major problems of the surface coating technique. This paper has proposed a new numerical method in order to tackle the traditional difficulties in simulating multi-crackings during constrained sintering. Main features of the present method include: (i) the material data is represented by an anisotropic constitutive law, (ii) a new numerical scheme is developed for the crack initialization and growth based on the material point method, (iii) the 3D viscous film shrinkage model is solved by using a dynamic FE scheme, and (iv) the random nature of the initial green body density is represented by statistical variabilities. It is shown that the model proposed by the present paper is capable for the nucleation and propagation of multi-cracks in a straightforward manner. Cracking patterns are shown to be consistent with experimental understandings. The focus of the paper is on the numerical issues and demonstrating the capacity of the model.