In electrical impedance tomography, one tries to recover the spatial conductivity distribution inside a body from boundary measurements of current and voltage. In many practically important situations, the object has known background conductivity but it is contaminated by inhomogeneities. The factorization method of Andreas Kirsch provides a tool for locating such inclusions. In earlier work, it has been shown, both theoretically and numerically, that the inhomogeneities can be characterized by the factorization technique if the input current can be controlled and the potential can be measured everywhere on the object boundary. However, in real-world electrode applications, one can only control the net currents through certain surface patches and measure the corresponding potentials on the electrodes. In this work, the factorization method is translated to the framework of the complete electrode model of electrical impedance tomography and its functionality is demonstrated through two-dimensional numerical experiments. Special attention is paid to the efficient implementation of the algorithm in polygonal domains.