An integrated methodology has been developed for computer simulation of electromagnetic scattering from large, nonperiodic, two-dimensional layouts of advanced photomasks (masks with optical proximity correction and phase-shifting masks). The domain decomposition method consists of three steps: First, by virtue of the linearity of the Kirchhoff-Fresnel diffraction integral, the mask layout is decomposed into a set of constituent single-opening masks. Second, the rigorous electromagnetic simulation of each three-dimensional structure from the set of these single-opening masks is circumvented and, instead, the result for the scattered field is synthesized based on two two-dimensional rigorous electromagnetic simulations that model the mask geometry in two cross-sectional planes. Subsequently, compact equivalent source models are used to describe the scattered fields on a reference plane. These models are constructed in such a way as to minimize the error in the part of the diffraction spectrum that passes through the projection system, allowing accurate and efficient image simulation. The normalized mean square error of the near scattered field is typically a fraction of 1% and speed-up factors for the total simulation time in excess of 400 (compared with the rigorous mask model) are achieved. The use of a look-up table approach facilitates orders of magnitude of further speed improvement.