# Simple, easy to install and use Python FEM solver (and example) for 2D cylindrical Laplace equation

In my earlier question Interpolating the gradient of a cylindrically symmetric potential field expected obey the Laplace equation, especially near/across r=0 I give a simple example of a Jacobi relaxation script - the Laplace equation is solved for a cylindrically symmetric electrostatic lens system. The symmetry reduces the problem to 2 dimensions.

I am reading and learning about FEM and rather than write my own I would like to start with an existing package to solve the same cylindrical lens problem (ideally in 2D for speed) with all Dirichlet boundary conditions (i.e. I define the potential on all surfaces and boundaries).

I have seen PySfe, pyGIMLi, Fenics and Firedrake and scikit-fem Most if not all of these have good examples for Laplace in 2D but these are cartesian examples.

What I need is

1. A recommendation for a package that is easy to install (I'm not a developer, I can spell "pip" and "conda" (barely) but not much more) and will have some way to generate solutions I can compare directly with my Jacobi relaxation work.
2. Advice how to modify the (likely present) Laplace example for the problem as 2D cylindrical, i.e. $$\phi(r, z)$$ and $$\theta$$-independent.

Ideally I'll set up the problem in an IDE and run from a command line. I'm using macOS with an Anaconda installation.

• Please provide the actual equations that you would like to solve and the necessary boundary conditions instead of a reference to your previous experimentation with this problem. Commented Feb 8 at 15:11
• @Konstantinos if I'm not mistaken (it's 2:30 AM here) Laplace's equation expressed in cylindrical coordinates with no $\theta$ dependence will be $$\frac{\partial^2 \phi}{\partial r^2} + \frac{1}{r} \frac{\partial \phi}{\partial r} + \frac{\partial^2 \phi}{dz^2} = 0.$$
– uhoh
Commented Feb 8 at 18:34

You did not provide all the details so I made some assumptions. First, assuming symmetry in $$\theta$$, the Laplace equation in cylindrical coordinates reads $$\partial_r^2 u + \frac{1}{r} \partial_r u + \partial^2_z u = 0.$$ Let $$\Omega = (0, 1)^2$$ where the variables are $$(r,z)$$. The symmetry boundary is $$\Gamma_0 = \{ (r, z) : r = 0 \}$$ and the remaining part is $$\Gamma_1 = \partial \Omega \setminus \Gamma_0$$. We consider the boundary conditions $$\partial_n u|_{\Gamma_0} = 0, \quad u|_{\Gamma_1} = g$$ for given boundary data $$g$$.

For testing, an analytical solution is given by the Bessel function of the first kind $$J_0$$, see a question on another forum. In particular, $$J_0(r)(\exp(z) + \exp(-z))$$ satisfies the above Laplace equation when its restriction onto $$\Gamma_1$$ is used as $$g$$.

We define the function spaces $$V_z = \{ w \in H^1(\Omega) : w|_{\Gamma_1} = z \}$$ for a given $$z$$. The weak formulation reads: find $$u \in V_g$$ such that $$\int_\Omega -\partial_r u \partial_r v + r^{-1} \partial_r u v - \partial_z u \partial_z v = 0 \quad \forall v \in V_0.$$

After pip install scikit-fem==9.0.1 you can implement the weak formulation as follows:

from skfem import *
from skfem.helpers import *

@BilinearForm
def laplace(u, v, w):
r, z = w.x
return - ur * vr + 1 / r * ur * v - uz * vz


Then you can solve the problem with the boundary data from the Bessel function:

from scipy.special import jv

m = MeshTri().refined(4)
basis = Basis(m, ElementTriP1())
A = laplace.assemble(basis)

y = basis.project(lambda x: jv(0, x[0]) * (np.exp(x[1]) + np.exp(-x[1])))
x = solve(*condense(A, 0 * y, x=y, D=basis.get_dofs({'top', 'right', 'bottom'})))



If you have matplotlib installed, running the above code will result in the following figure:

You can also plot the difference x - y to see that the error is small:

I was a bit suspicious of the larger error at $$r=0$$ so I tried different mesh refinements to see that it converges as expected. I concluded that the error is probably due to the term with $$1/r$$ which warrants for more integration points or more elements near $$r=0$$.

(Note: I am the maintainer of scikit-fem.)

• This is exactly what I need to get started. I spent yesterday at the library reading books and viewing online material but found myself saying "OK I believe you, very rigorous, but how do I do it?" A quick installation question; yesterday I'd already used pip install scikit-fem and pip install meshio (I have matplotlib from my anaconda installation) and while your example runs fine, skfem.__version__ yields 8.1.0 If I use your command for 9.0.1 and add  -U should that bring me up to date?
– uhoh
Commented Feb 9 at 2:20
• Yes, that should work.
– knl
Commented Feb 9 at 4:54