A Mathematical Model for Electrostatics as a Literate FreeFem++ Program

using FreeFem++-js 17.1

Summary The purpose of this document is to demonstrate literate programming with FreeFem++-js. It is a modified version of subsection 9.1.2 "Electrostatics" from the FreeFem++ Manual. In this version, FreeFem++ commands are operational and generate numerical results, including pictures, directly into the document. Some commands are modifiable to view the effect of different parameter values.

We assume that there is no current and a time independent charge distribution. Then the electric field \(\vec E\) satisfies

\begin{equation} \label{eq:Maxwell} \mathrm{div}\vec E=\rho/\epsilon,\quad \mathrm{curl}\vec E=0 \end{equation}

where \(\rho\) is the charge density and \(\epsilon\) is called the permittivity of free space. From the second equation in \eqref{eq:Maxwell}, we can introduce the electrostatic potential such that \(\vec E=-\nabla \phi\). Then we have Poisson's equation \(-\Delta \phi=f\), \(f=-\rho/\epsilon\). We now obtain the equipotential line which is the level curve of \(\phi\), when there are no charges except conductors \(\{C_i\}_{1,\cdots,K}\). Let us assume \(K\) conductors \(C_1,\cdots,C_K\) within an enclosure \(C_0\). Each one is held at an electrostatic potential \(\varphi_i\). We assume that the enclosure \(C_0\) is held at potential 0. In order to know \(\varphi(x)\) at any point \(x\) of the domain \(\Omega\), we must solve

$$-\Delta \varphi =0\quad \textrm{in }\Omega,$$

where \(\Omega\) is the interior of \(C_0\) minus the conductors \(C_i\), and \(\Gamma\) is the boundary of \(\Omega\), that is \(\sum_{i=0}^N C_i\). Here \(g\) is any function of \(x\) equal to \(\varphi_i\) on \(C_i\) and to 0 on \(C_0\). The boundary equation is a reduced form for:

$$\varphi =\varphi _{i}\;\text{on }C_{i},\;i=1...N,\varphi =0\;\text{on }C_{0}.$$


First we give the geometrical informations. Here they are in mathematical and in FreeFem++ form. The FreeFem++ definition of $C_1$ and $C_2$ can be modified to view the effect of different geometries.

\(C_0=\{(x,y);\; x^2+y^2=5^2\}\) is a circle with center at (0,0) and radius 5,

border C0(t=0,2*pi) { x = 5 * cos(t); y = 5 * sin(t); }

\(C_1=\{(x,y):\; \frac{1}{0.3^2}(x-2)^2+\frac{1}{3^2}y^2=1\}\) first elliptical hole,

\(C_2=\{(x,y):\; \frac{1}{0.3^2}(x+2)^2+\frac{1}{3^2}y^2=1\}\) second elliptical hole.

Let \(\Omega\) be the disk enclosed by \(C_0\) with the elliptical holes enclosed by \(C_1\) and \(C_2\). Note that \(C_0\) is described counterclockwise, whereas the elliptical holes are described clockwise, because the boundary must be oriented so that the computational domain is to its left.

mesh Th = buildmesh(C0(60)+C1(-50)+C2(-50));
plot(Th,ps="electroMesh.eps"); // see figure [electroMesh]

P1 FE-space, unknown and test function :

fespace Vh(Th,P1);
Vh uh,vh;

Definition of the problem :

problem Electro(uh,vh) = int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) // bilinear
+ on(C0,uh=0) // boundary condition on $C_0$
+ on(C1,uh=1) // +1 volt on $C_1$
+ on(C2,uh=-1); // -1 volt on $C_2$

Solve the problem :

plot(uh,ps="electro.eps",wait=true); // see figure [electro]

Figure [electroMesh]:
Disk with two elliptical holes
Figure [electro]:
Equipotential lines, where \(C_1\) is located in the right hand side