A direct numerical simulation method based on the Maxwell stress tensor and a fictitious domain method has been developed to solve flows with suspended paramagnetic particles. The numerical scheme enables us to take into account both hydrodynamic and magnetic interactions between particles in a fully coupled manner. Particles are assumed to be non-Brownian with negligible inertia. Rigid body motions of particles in 2D are described by a rigid-ring description implemented by Lagrange multipliers. The magnetic force, acting on the particles due to magnetic fields, is represented by the divergence of the Maxwell stress tensor, which acts as a body force added to the momentum balance equation. Focusing on two-dimensional problems, we solve a single-particle problem for verification. With the magnetic force working on the particle, the proper number of collocation points is found to be two points per element. The convergence with mesh refinement is verified by comparing results from regular mesh problems with those from a boundary-fitted mesh problem as references. We apply the developed method to two application problems: two-particle interaction in a uniform magnetic field and the motion of a magnetic chain in a rotating field, demonstrating the capability of the method to tackle general problems. In the motion of a magnetic chain, especially, the deformation pattern at break-up is similar to the experimentally observed one. The present formulation can be extended to three-dimensional and viscoelastic flow problems.