This paper introduces a two-strain epidemic model with reaction-diffusion process, which is derived from the probabilistic viewpoint using master equations. This model encompasses situations where two distinct strains coexist and interact within a population. Using the diffuse domain method, we easily embed the model with a complex spatial domain into a large regular domain by reformulating it using a phase-field-like variable. The equivalence between the reformulated system and the original system of reaction-diffusion equations with Neumann boundary conditions are shown through both matched asymptotic analysis and numerical validations. The diffuse-domain-based epidemic model is numerically solved using a semi-implicit finite difference scheme, giving rise to an easy simulation of two-strain coexistence/extinction in complex geometric regions. Furthermore, we compute the basic reproduction numbers of both strains and analyze their dependencies on the shapes of domains and the parameters of the model.