Let us have a mesh:

n=60
h = 4/(n-1)
x = seq(-2,2,len=n)
p = expand.grid(x=x,y=x)
r = sqrt(p$x**2 + p$y**2)

Let us now define some helper functions:

d_1 = c(2:n,1)
d_2 = c(n,2:n-1)
dx = function(a) {a = matrix(a,n); as.vector((a[d_1,] - a[d_2,])/(2*h))}
dy = function(a) {a = matrix(a,n); as.vector((a[,d_1] - a[,d_2])/(2*h))}
lens = function(v) sqrt(rowSums(v**2))
plot.vec = function(p,v,scale=1,...) {
  plot(p,type="n",asp=1,...)
  segments(p[,1],p[,2],p[,1]+v[,1]*scale,p[,2]+v[,2]*scale)
}
plot.field = function(x,...) filled.contour(matrix(x,n), color.palette = rainbow, asp=1, ...)

Let us now define a \(\alpha\) field:

theta = h/4
a = (1-r)/theta
alpha = exp(a)/(1+exp(a))
plot.field(alpha,main=expression(alpha))

Let us have the gradient of \(\alpha\):

v = cbind(dx(alpha),dy(alpha))
plot.vec(p,v,0.02,main=expression(nabla * alpha))

And normal vector:

lv = sqrt(rowSums(v**2))
nv = v/lv
plot.vec(p,nv,h,main="n")

Finally the curvature \(\kappa\):

kappa = dx(nv[,1])+dy(nv[,2])
kappa[is.nan(kappa)] = 0
plot.field(kappa, main=expression(kappa))

Let us now compute the force field:

f = kappa * v
plot.vec(p,f,0.02,main=expression(f == kappa * nabla * alpha))

One can clearly see that it have a non-zero rotation:

rot = dy(f[,1]) - dx(f[,2])
plot.field(rot, main=expression(rot ~ f))

Let us now solve an equation: \[\Delta\phi = \text{div} f\] So that \(\text{grad}\phi\) will be the best potential field approximating \(f\):

div = dx(f[,1]) + dy(f[,2])
plot.field(div, main=expression(div ~ f))

Let us use the fourier transform. For convoltion we have: \(\mathcal{F}(f\star g) =\mathcal{F}(f)\cdot\mathcal{F}(g)\). The finite difference operator of laplacian: \[(\Delta u)_{i\;j} \simeq \frac{u_{i-1\;j}+u_{i+1\;j}+u_{i\;j-1}+u_{i\;j+1}-4u_{i\;j}}{h^2}\] can be seen as convolution with a field: \[m = \frac{1}{h^2}\left[\begin{array}{ccc} & 1 & \\ 1 & -4 & 1 \\ & 1 & \end{array}\right]\] Thus in the fourier space the operation of taking a laplasian of \(u\) is a multiplication of the fourier transform of \(u\) by the fourier transform of this field \(m\). Finally this operation can be reversed and we can divide the fourier transform of a given field to solve a Poisson equation:

div = matrix(div,n) #make a matrix for 2D FFT
m = matrix(0,n,n)
m[1,1] = -4
m[2,1] = 1
m[1,2] = 1
m[n,1] = 1
m[1,n] = 1
m = m / (h^2)
vf = fft(div)/fft(m)
vf[1] = 0
phi = Re(fft(vf,inverse = TRUE))/length(vf) # reverse FFT
plot.field(phi,main=expression(phi))

Let us now see the gradient of \(\phi\):

dp = cbind(dx(phi),dy(phi))
plot.vec(p, dp, 0.02, main=expression(nabla*phi))

Let us now see the difference between \(f\) and \(\text{grad}\phi\):

f2 = dp - f
plot.vec(p,f2,0.02, main=expression(nabla*phi - f))

Apart form the high values at the interface it is hard to see the vectors. We will normalize the high vectors

maxlen = h
sel = lens(f2)>maxlen
f2[sel] = f2[sel]/lens(f2)[sel]*maxlen
plot.vec(p,f2,1)