I have some molecules in different positions of a plane. Each molecule has a number which is shown by $\hat{f}$. I want to find an continues field is written as:
$\int |\nabla \cdot f|^2 d^2r$.
Integration is the whole plane. I need to discretize formula and sum over different molecules. But I don't know how to do that. Could anyone help me?
$\newcommand{\vect}[1]{{\bf #1}}$
I will assume that you have $N$ molecules each at a location $\vect{r}_j\in \Omega$, we can define the density at an arbitrary location $\vect{r}$ as the weighted sum of the masses of the $N$ molecules. Of course, this weight is distance dependent: the farther a molecule is, the smaller is its contribution
$$ \rho(\vect{r}) = \sum_{j = 1}^N m_j W_h(|\vect{r} - \vect{r}_j|) \tag{1} $$
where $W_h\in C^{(1)}(\Omega)$ is then a decaying function of $r = |\vect{r}|$ with finite support, in this case is of size $h$. Below I give a possible choice for $W_h$, in the mean time, note that it must satisfy
$$ \sum_j m_j = M(\Omega) = \int_\Omega{\rm d}\Omega ~\rho(\vect{r}) = \sum_jm_j\int_\Omega{\rm d}\Omega W_h(|\vect{r} - \vect{r}_j|) ~~~\Rightarrow~~~ 2\pi \int_{-\infty}^{+\infty}{\rm d}r ~rW_h(r) = 1 \tag{2} $$
where I assume that $\Omega \subset \mathbb{R}^2$, but clearly you can extend this to the dimensions you want. Now, with equation (1) you can calculate the density associated to any molecule $i$
$$ \rho(\vect{r}_i) = \sum_{j = 1}^N m_j W_h(|\vect{r}_i - \vect{r}_j|) \tag{3} $$
The size of the kernel is usually calculated as the distance $h$ that encloses the $k$-th nearest neighbor of the particle $i$. From my experience this is the slowest part of the algorithm, but can be done really efficiently in $\mathcal{O}(N\log N)$.
In general any field can be written as
$$ f(\vect{r}) = \sum_j m_j \frac{f_j}{\rho_j}W_h(|\vect{r} - \vect{r}_j|) \tag{4} $$
In particular $f = \{n_x, n_y\}$. Now it is a matter of calculating integrals, just note that
$$ \nabla f(\vect{r}) = \sum_j m_j \frac{f_j}{\rho_j}\nabla W_h(|\vect{r} - \vect{r}_j|) = \sum_j m_j \frac{f_j}{\rho_j}\frac{\vect{r} - \vect{r}_j}{|\vect{r} - \vect{r}_j|}W_h'(|\vect{r} - \vect{r}_j|) \tag{5} $$
I will leve the rest to you, $\nabla \times \vect{n} = (\ldots)$
This is one popular choice $$ W_h(x) = \frac{10}{7\pi h^2}\left\{\begin{array}{ll} 1 - (3/2)z^2(1 - z^2/2) & 0 \leq z = x/h < 1 \\ (2 - z)^3/4 & 1 \leq z = x/h < 2 \\ 0 & z = x/h \geq 2 \end{array} \right. $$
If you're familiar with python you can try pysph, it implements all the equations above in a fairly transparent way