Abstract: Simulation system and method for a reservoir represented by a model having a plurality of matrix nodes, a plurality of fracture nodes, a fracture zone and a fracture-free zone are described. One embodiment includes computer-implemented steps of, for each of a plurality of matrix nodes in interconnected zones, characterizing the matrix nodes as dual matrix nodes or non-dual matrix nodes. Dual matrix nodes may be treated as a block and eliminated, after which the coefficient matrix and RHS vector for the model are updated. In-fill patterns between and among fracture nodes and non-dual matrix nodes resulting from the elimination are determined. For non-dual matrix nodes, the nodes are merged to a fracture grid, after which the coefficient matrix and RHS vector for the model are updated accordingly. The resultant linear system is linearly solved, and backsolving is performed for any eliminated dual matrix nodes.