A Lagrangian vorticity-based method is presented for simulating two-way phase interaction in a two-phase flow with heavy particles. The flow is computed by solving the vorticity transport equation, including the particle-induced vorticity source, and the mass conservation equation for particle concentration on separate sets of fluid control points and particle control points, respectively. The fluid control points are advected with the local fluid velocity, plus a diffusion velocity for viscous problems to account for the spread of the vorticity support via diffusion, while the particle control points are advected by solution of the Lagrangian particle momentum equation. The particle concentration and vorticity transport equations are evaluated using volumeaveraged particle velocity and contact force fields, obtained by a weighted average over nearby particle control points. One novel feature of the numerical method is the scheme for calculation of the particle-induced vorticity source using a "moving least-square" differentiation scheme across the two sets of control points. Another feature of the method is its ability to absorb the vorticity generated by particle forces through an adaptive scheme for generation of new fluid control points. Test calculations with a vortex patch filled with particles show that the numerical results compare well with the results obtained both by a traditional finite-difference method and by an asymptotic approximation valid for small Stokes numbers. Other features of the numerical method are demonstrated for calculations involving a particle cloud falling under gravity and a two-phase mixing layer flow.