Field data suggest that hydraulic fractures under in situ conditions are three-dimensional and non-planar in shale gas reservoirs. In this study, a fully coupled deformation, fracture growth, and fluid flow model is developed to simulate this complex phenomenon with a focus on the reduction of the total computational effort. An extended finite element method with high-order tip enrichment functions is used to simulate the rock deformation while the finite element method for the laminar flow in fractures. The explicit fracture representation is adopted to describe the non-planar 3D fractures. Schemes for identifying the fluid elements, selecting the enriched nodes, and subdividing the enriched elements are presented. The Newton-Raphson scheme is used to solve the coupled equilibrium and fluid flow equations in an element-by-element fashion. Within each Newton-Raphson iteration step, only enriched degrees of freedom are involved in the solution process by using the reduction technique to further reduce the computational effort. A displacement extrapolation method considering high-order terms is proposed to extract stress intensity factors from the displacement field. Several examples are performed to demonstrate the applicability and effectiveness of the proposed approach. The effectiveness of MPI parallel implementation of the proposed method is also investigated.