and $P=F(\psi)-\frac{1 }{2}|\nabla\psi|^2$ is a stationary solution to the Euler equations, where $(\nabla\psi)^\bot:=(\frac{\partial\psi}{\partial
x_2},-\frac{\partial\psi}{\partial x_1}), F(t)=\int_0^tf(s)ds$. Moreover, the velocity $\mathbf{v}$ is irrotational on the set where $f(\psi)=0$.

Set $q=-\psi_0$ and $u=\psi-\psi_0$, then $u$ satisfies the following boundary value problem $$\label{1.3}
\begin{cases}
-\Delta u=f(u-q),\quad & x\in\Omega,\\
u=0,\quad &x\in\partial\Omega.
\end{cases}$$ In addition, if we suppose that $\inf_\Omega q>0$ and $f(t)=0,~t\leq0$, the vorticity set $\{x: f(\psi)>0\}$ is bounded away from the boundary.

The motivation to study is to justify the weak formulation for point vortex solutions of the incompressible Euler equations by approximating these solutions with classical solutions.

Marchioro and Pulvirenti [@MP] have approximated these solutions on finite time intervals by considering regularized initial data for the vorticity. On the other hand, the stationary point vortex solutions can also be approximated by stationary classical solutions. See e.g. [@BF1; @N; @SV; @T; @Y1] and the references therein. It is worth pointing out that the above approximations can just give explanation for the formulation to single point vortex solutions. In this paper, we will show that multi-point vortex solutions can be approximated by stationary classical solutions. There are many results for problem on the existence and asymptotic behavior of solutions under various assumptions. In [@Ba; 