The boundary element method was newly implemented into an optical lithography simulator so that it could evaluate rigorously the topological effects of 2dimensional phase-shifting masks. Both transparent and periodic boundary conditions were applied for the method, and the continuity conditions were used for treating interface nodes. The accuracy of the module developed for simulating aerial images was verified by comparison with analytic solutions and published results. In addition, it was found that our simulator would be more efficient than the conventional method based on the rigorous coupled wave analysis in views of the convergence and the calculation speed. Finally, the optimal design of two phase-shifting masks was performed.