Creating a stimulated reservoir domain (SRD) is a necessity to effectively extract natural gas from shale reservoirs. After hydraulic fracturing, a shale reservoir has three distinct domains: SRD, non-stimulated reservoir domain (NSRD) and hydraulic fractures (HFs). In previous studies, the property contrasts and interactions between different domains are often not fully considered. In this study, a fully coupled multi-domain and multi-physics model is developed to incorporate these complexities. The shale reservoir is characterized as an assembly of three distinct components: organic kerogen, inorganic matrix and HF. Furthermore, the kerogen and inorganic matrix are defined as dual-porosity-dual-permeability media, while the HF is simplified as a 1-D cracked medium. Particularly, the inorganic matrix has different properties in each of the SRD and NSRD to reflect the stimulation effect of hydraulic fracturing on the near-HF matrix. Under this framework, a series of partial differential equations (PDEs) fully coupled by mass transfer and mechanical deformation relations were derived to define various processes in the shale reservoir. These PDEs were numerically solved by the finite element method. The proposed model is validated against analytical solutions and verified against gas production data from the field. Sensitivity analyses reveal: (1) that both size and internal structure of the SRD significantly affect gas extraction by improving SRD properties and in creating low-pressure zones around the NSRD; (2) that the NSRD determines the sustainability of gas production; and (3) that the change of mechanical properties in one domain affects the evolution of transport properties in the entire reservoir.