Abstract Modelling fluid flow through a porous rock mass containing fractures is a common issue in many applications in civil engineering, mining engineering, and extraction of energy resources such as petroleum, gas, and geothermal heat. Fractures in these applications normally play the dominant role in conducting the fluid through the media but the contribution from the porous rock cannot be ignored, particularly when the relative conductivity is significantly high. In this case, flow through both the fractures and the porous rock must be considered simultaneously and their interactions could potentially be very complex. The Dual Porosity/Permeability Model (DPM) is commonly used to address the problem, however it has in general the limitation of over-simplified representation of fractures and fracture network within the rock mass. In this paper, for the solution of the fluid flow in such media, a numerical modelling approach based on the application of the Boundary Element Method (BEM) on the fluid flow in the block in combination with the Finite Element Method (FEM) on the fractures is presented, where fractures are represented explicitly in the model as discretised boundaries. A code dedicated to this approach was developed and is presented in this paper, together with results from illustrative examples demonstrating the effectiveness of the approach. The proposed method has the additional benefit of reducing computational costs, which is particularly useful for cases with large-scale fracture networks embedded in a conductive rock matrix. The source code is available for downloading at the link: https://github.com/cx-adelaide/BEMFEM_FlowSim.