1. A method for performing simulation of a reservoir represented by a model having a plurality of matrix nodes, a plurality of fracture nodes, at least one fracture zone and at least one fracture-free zone, the method comprising computer-implemented steps of: for each of a plurality of matrix nodes comprising interconnected zones, characterizing the matrix nodes; for a block comprised of one or more dual matrix nodes, eliminating the block and updating a coefficient matrix and RHS vector for the model accordingly; determining in-fill patterns between fracture nodes, including virtual fracture nodes, that will result from the elimination; for one or more virtual fracture nodes, merging the virtual fracture nodes to a fracture grid and updating the coefficient matrix and RHS vector for the model accordingly; solving the resultant linear system using a linear solver; and backsolving for the eliminated dual matrix nodes.