FEniCS is an open source general purpose finite element solver. That means it incredibly general but also not as easy to start with as commercial finite element solvers; which are often tailored for specific applications. The purpose of the FEniCS electrostatics posts is to demonstrate how FEniCS can be used to solve electrostatics problems specifically. That means that I won’t focus too much on the FEniCS syntax or settings; these are nicely documented elsewhere, although I will highlight the main points. Importantly FEniCS can be controlled with python scripts.

This is my third electrostatics post. The first found the electric field between two infinite parallel conducting plates with a potential difference. This was achieved by making a rectangular mesh from within a FEniCS script, the boundary conditions were specified with python functions and the Laplace equation was solved for the space between the plates. The second post built on the first post by constructing a 2D geometry with Gmsh and specifying the boundaries within the Gmsh GUI. These boundaries were then imported into FEniCS allowing the Laplace equation to be solved without any functions being written to identify the boundaries.

This post will demonstrate how to solve a basic 3D electrostatic problem with FEniCS. Again the problem being solved is the Laplace equation, however this time the geometry will be a conducting sphere held inside a spherical shell, the space between these spheres will be empty space which is where the field will be calculated. This problem; like the previous examples, has an analytical solution which will allow the accuracy of the finite-element result to be assessed. What’s special about this post is that the geometry is going to be 3D, which allows general-purpose CAD software to be used to design the geometry; in this case I will use FreeCAD [1]. This means that the method is incredibly versatile, the exact same procedure can be used for much more complex geometries.

Let’s start!

Analytical Solution

The method to solving this electrostatic problem will be to use Gauss’ law in differential form

coupled with the fact that the curl of an electrostatic field is zero [2],

which means that the electric field is uniquely defined given only the divergence. This allows the electric field to be specified as the gradient of a scalar field known as the electric potential

where Φ is the potential. This can then be inserted into Gauss’ law to give

which is the Poisson equation with the “source” being particles with an electric charge.

For this particular problem there are no “sources” in the free space surrounding the conducting central sphere so the problem boils down to the Laplace equation

which has to be solved for certain boundary conditions. In this case spherical polar coordinate are the most appropriate, giving the Laplace operator the following form

which will now have to be solved. Due to the symmetry of the problem there will be no variation of Φ with φ or θ. This simplifies the Laplace equation to

using the product rule on this equation gives

which is an Euler-Cauchy equation; solutions can thus be looked up [3].

The solution is generally of the form

into which the boundary conditions must be input to find the constants c1 and c2. At the inner radius r=ri, Φ = Φi, at the outer radius r = ro, Φ = Φo, putting these in, solving for the constants and substituting them back into the general equation gives the following solution

which can then be used to find the electric field from the gradient

the final solution for this problem. A wxmaxima workbook is available for this derivation and can be found in the repository linked to at the end of the post.

FEniCS Solution

Geometry

In this example the geometry will be constructed with freeCAD, which is a free and open source computer aided design environment. It can be downloaded from the website [3] and works on all major computing platforms at the time of writing. The whole process of making the geometry with can be viewed here, alternatively I’ve included a full set of instructions.

Step 1 – Opening Part Environment

After opening FreeCAD a blank screen should appear, as shown in Figure 1.

Switch to part mode, which is a mode that allows the construction of simple 3D geometries with ease; see Figure 2. If designing a more complex geometry than that of this post, then you might have to go down the route of drawing 2D sketches and extruding etc.

Step 2 – Drawing the Outer Sphere

This sphere will represent the vacuum inside which FEniCS will solve, the outer surface will be where out conducting shell sits. To add the sphere click the add sphere button which is located next to the cursor in Figure 3, a sphere should be displayed. This will have some default dimensions which will be changed in the next step.

Select the sphere object in the object list and look at the Data tab where the dimensions of the sphere can be changed. In this case the only thing that has to be set is the radius to 100 mm and the centre to be at the origin. The sphere can also be given a name if desired, these settings are shown in Figure 4.

Step 3 – Adding the Inner Sphere

The central sphere must be added, although again it is only the surface of this sphere which is important for applying boundary conditions. To do this click the add sphere button again. This time the radius should be smaller than that of the outer sphere, in this example 10 mm was used as shown in Figure 5.

Step 4 – Subtracting the Spheres

At the moment there are two spheres, however we need to have a single object which will represent the whole space of this problem. If we subtract the inner sphere from the outer sphere then the resulting shape will have one volume; a “hollow” sphere, and two surfaces – the inner and outer shells. This is ideal as the boundary conditions can be applied to the surfaces and the field will be solved inside the volume.

In order to subtract the inner sphere from the outer a boolean operation must be performed. Select the outer sphere from the object list, and then the inner sphere (the order is important) then click the subtract operation; indicated by the cursor in Figure 6.

After performing this operation the object won’t look different as only the central part has been removed. By using a clipping plane (View > clipping plane) it becomes possible to see inside the object; as shown in Figure 7.

Step 5 – Export the Geometry

The geometry is now complete and the file just has to be exported into a format that Gmsh can understand. Select the object, then click file > export, choose a filename and the step file format can been used; as shown in Figure 8. I also advise saving the FreeCAD file itself as this will make it easier to modify in future.

Mesh

Now that the geometry has been made a mesh must now be generated for the hollow sphere. The whole process is described in this video, and each step is described below.

Step 1 – Import Geometry

Open GMsh and import the geometry which was exported from FreeCAD. To do this click file > open and locate the step file. If surfaces are visible then the page should show the outer and inner spheres as in Figure 9.

Step 2 – Identify Boundaries and Subdomains

The geometry consists of one volume and two surfaces, these will be identified in GMsh so that FEniCS can automatically identify the boundaries and subdomains. This process is made easier by being able to view only certain surfaces in GMsh. To do this use Tools > visibility > Tree browser, selecting an object and clicking apply; only the selected object will remain visible (see the GIF below for an example).

Make only the inner sphere visible so that the window resembles Figure 10. Add a new surface group (modules > geometry > physical groups > add > surface), call it “inner” and give it the number 1. Select the inner surface, as shown in Figure 10, then press the letter e on the keyboard to finish adding surfaces to this group.

Next make only the outer sphere visible, as shown in Figure 11. Add another surface group, but this time call it “outer” and give it the number 2; as shown in Figure 12. Select the outer surface; as shown in Figure 11, and then press “e” on the keyboard to finish adding surfaces to this group. All surfaces have now been identified so press “q” on the keyboard to stop adding surfaces.

The last object to add to a group is the volume, make the only volume visible from the Tree Browser and in the Physical Groups menu click add volume. Select the single yellow sphere that appears. Give it the label vacuum and the number 3 as shown in Figure 14

Step 3 – Improving Mesh Density

All of the necessary steps have been completed, however the default mesh is quite coarse and does not properly describe the geometry; see Figure 15. To improve the mesh some transfinite lines will be specified on the inner and outer spheres.

Add a new transfinite line (mesh > define > transfinite > line) and select the line surrounding the inner sphere as shown in Figure 16. Enter 15 points for this line, as shown in Figure 17 and press “e” on the keyboard to end the selection.

Next add the line surrounding the outer sphere as a transfinite line, and choose 30 points; see Figure 18.

Step 4 – Generate the Mesh

Press the 3D button (mesh > 3D) to generate a mesh, from the outside this will resemble Figure 19. To view the inside and check the inner sphere has been appropriately meshed use a clipping plane (tools > clipping plane > mesh); see Figure 20. Click the save button to save this mesh.

Step 5 – Convert the mesh to the FEniCS Format

In a terminal window move to the directory where the mesh is saved, then use the command (replace input/ output with appropriate filenames)

dolfin-convert input.msh output.xml

the output for this command for me is shown in Figure 21. The program dolfin-convert is packaged with FEniCS. If this has completed successfully three xml files should be present geo.xml, geo_physical_region.xml and geo_facet_region.xml, and in this case the mesh is now ready to be used with FEniCS.

Solution

From now the code is exactly the same as the code from the previous post the inner sphere will be set at 10 V and the outer sphere will be set at 1V. The boundary conditions are set using the numbers used when creating the physical groups in gmsh. The only difference is that in the 2D problem the boundaries were set on lines, in this 3D example the boundaries will be the surfaces. The full code is listed below.

import fenics as fn # Import the mesh, identify the subdomains and boundaries # and then make a function space and a dx element. mesh = fn.Mesh('geo.xml') markers = fn.MeshFunction("size_t", mesh, 'geo_physical_region.xml') boundaries = fn.MeshFunction('size_t', mesh, 'geo_facet_region.xml') dx = fn.Measure('dx', domain=mesh, subdomain_data=markers) V = fn.FunctionSpace(mesh, 'CG', 1) # Set the boundary conditions using the boundary numbers # specified in Gmsh. inner_boundary = fn.DirichletBC(V, fn.Constant(10.0), boundaries, 1) outer_boundary = fn.DirichletBC(V, fn.Constant(1), boundaries, 2) bcs =[inner_boundary, outer_boundary] # Solve the Poisson equation with the source set to 0 # (the Laplace equation) via the L = fn.Constant('0')... line u = fn.TrialFunction(V) v = fn.TestFunction(V) a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx L = fn.Constant('0') * v * fn.dx u = fn.Function(V) fn.solve(a == L, u, bcs) # Find the electric field by taking the negative of the # gradient of the electric potential. Project this onto # the function space for evaluation. electric_field = fn.project(-fn.grad(u)) # output the results to two separate files potentialFile = fn.File('output/potential.pvd') potentialFile << u e_fieldfile = fn.File('output/e_field.pvd') e_fieldfile << electric_field

Results

As before the results files can be viewed using paraview, or plotted with matplotlib depending on your FEniCS set-up. To plot with matplotlib use

fn.plot(u) fn.plot(electric_field)

however for 3D problems I have personally found paraview to be more useful. The potential solution is shown in Figure 22, and the electric field is overlaid in Figure 23.

by eye these results appear sensible, however they can also be compared with the analytical solution. After solving the objects u and electric_field become functions which can be evaluated at any position by specifying x, y and z coordinates. If the field is a scalar this will return a single number, if the field is a vector then the x, y and z components will be returned for example

In[15]: electric_field(15, 23, 60) Out[15]: array([ 0.00513091, 0.00848525, 0.02177435])

note that the fn.project() function must be used for this to work; does not need to be manually performed on the solution u. The field can then be evaluated along a line using the python map function which runs a given function for an input list and returns an array of the results. For example

x = np.linspace(10, 99, 1000) y = np.zeros(len(x)) z = np.zeros(len(x)) coords = np.dstack((x, y, z))[0] u_line = np.array(list(map(u, coords)))

will evaluate u from 10 to 99 in 1000 steps; numpy (np) has been used for generating the coordinates. The result of this is shown in Figure 24, and is overlaid with the analytical answer; the percent difference is also shown.

The electric field can also be compared with the analytical solution by using the map command. By choosing to evaluate the function over an on-axis line only a single component need be compared which has been done and is shown in Figure 25.

The solution could be improved by increasing the mesh density, and since Gmsh has a command line interface the whole process of mesh refinement could be automated with a specific convergence criteria. Alternatively the order of the FunctionSpace can be increased by changing 1 to 2 in the line

V = fn.FunctionSpace(mesh, 'CG', 1)

which gives the results shown in Figures 26 and 27. This will be the final solution of this post.

Conclusion

The electric field for two concentric spheres has been calculated with FEniCS which; with the settings described, achieved around 5% agreement with the analytical solution. Importantly a sphere was chosen for this example so that an analytical solution would exist for comparison, but the sphere could just as easily be replaced with a geometry for which there is no clear analytical solution.

All files used for this example can be found in this git repository, including a wxmaxima file with the analytical solution, the FreeCAD file and exported geometry, the Gmsh geometry file and mesh, a complete python script for the solution, a jupyter notebook which could be used instead of the python script and the output files.

I like this example because it really could allow FEniCS to be used for something quite practical. It is also incredibly easy: the python code is the same as for a 2D problem, and all boundaries are identified using a usable GUI tool.

In future posts I’ll try to extend this to include linear materials as well as charge distributions – rather than just perfect conductors as boundary conditions. I’d also like to demonstrate some magnetostatics and to use these techniques for something RF related; perhaps finding the characteristic impedance of a uniform transmission line. If you’ve had a go let me know how it went and whether you’ve had any problems or found anything difficult – I can try to improve my instructions. If there is anything you’d like me to try and cover just let me know and I’ll give it a go; I’m learning to use FEniCS as I go with these posts! I can be contacted via the comments on this page or Mastodon ( @comphys@qoto.org ).

References

[1] – FreeCAD – https://www.freecadweb.org/

[2] – Griffiths, David J. Introduction to Electrodynamics. 4 edition. Cambridge, United Kingdom ; New York, NY: Cambridge University Press, 2017.

[3] – Boas, Mary L. Mathematical Methods in the Physical Sciences. 3rd edition. Hoboken, NJ: John Wiley & Sons, 2005.