The author deals with a stochastic impulse control problem in an open bounded domain \(\Omega\subset\mathbb R^N\). The state of the controlled system is given by the solution \(X_t^x\) of the equation \[ dX_t^x=b(X_t^x,\alpha_t)\,dt+\sigma(X_t^x,\alpha_t)\,dW_t, \quad t\in]\theta_i\theta_{i+1}[,\quad X_{\theta_i}^x=X_{\theta_i^-}^x+\xi_i,\quad X_{0}^x=x\in\Omega, \] where \(W_t\) is a \(D\)-dimensional Brownian motion, \(\alpha_t\) is an adapted process with values in a compact set \(\mathcal A\subset\mathbb R^p\), \(\theta_i\) is an increasing sequence of stopping times (both with respect to the filtration generated by the Brownian motion \(W_t\)), \(\xi_i\) is a sequence in \(\mathbb R_+^N\) and \(b,\sigma\) are continuous functions on \(\overline\Omega\times\mathcal A\) taking values respectively in \(\mathbb R^N\) and in the space of \(N\times D\) real matrices. The control is the triplet \(A=(\alpha_t,\theta_i,\xi_i)\). For any fixed \(x\) and \(A\), the cost function is \[ \mathcal T(x,A)\colon=\mathbb E\left[ \int_0^{\tau_x}f(X_t^x,\alpha_t)e^{-\lambda t}\,dt+ \sum_{i\leq\tau_x}(c(\xi_i)+k)e^{-\lambda\theta_i}+ {\mathbb I}_{\tau_x<\infty}\varphi(X_t^x)e^{-\lambda\tau_x} \right], \] where \(f\) and \(c\) are real continuous functions, \(\lambda\) and \(k\) are positive numbers, and \(\tau_x\) stands for the first exit time from \(\overline\Omega\) of the trajectory \(X^x_t\). The value function \(U\) is of the type \[ U(x):=\inf_{A}\mathcal T(x,A). \] The dynamic programming principle relates this optimization problem to a Hamilton-Jacobi-Bellman equation, which now stands in the quasi-variational inequality \[ \max\{H(x,u,Du,D^2u),u-Mu\}=0\quad\text{on}\;\Omega, \] coupled with the Dirichlet boundary condition \[ u=\varphi\quad\text{on}\;\partial\Omega. \] Here, \(H\) is the (possibly degenerate) elliptic operator \[ H(x,u,p,X)\colon=\sup_{\alpha\in\mathcal A}\left\{ -\frac{1}{2}\text{tr}[a(x,\alpha)X]-b(x,\alpha)p+\lambda u-f(x,\alpha) \right\} \] for any \(x\in\overline\Omega\), \(p\in\mathbb R^n\), \(u\in\mathbb R\), and \(X\in{\mathcal S}^N\) (the space of \(N\times N\) symmetric matrices), where \(a(x,\alpha)=\sigma(x,\alpha) \sigma^T(x,\alpha)\), M is a nonlocal operator of the form \[ Mu(x):=\inf\{u(x+\xi)+c(\xi)+k : \xi\in\mathbb R^N_+,x+\xi\in\overline\Omega\}. \] The boundary condition in the viscosity solutions sense does not identify a unique solution. The discontinuities of solutions at the boundary of the domain play an essential role and cannot be removed. The author proposes a selection criterion which, enforcing the information coming from the boundary, selects the value function among all possible viscosity solutions. In addition, a monotone iterative scheme which approximates the value function is produced.