Gerrit Mur, retired, but not quite
Summary
I look back at the discussion about edge elements for electromagnetics. Not much progress seems to have been made since I published 'the Fallacy' and the points against edge elements I made in that paper remain to be refuted. In this paper I shortly describe the situation and propose a new finite-element method for computing electromagnetic fields. Unfortunately I do not have the means any more to put them to the test.
1. What were the problems?
a. The problems with edge elements
In my paper 'The fallacy of edge elements' I argued against the use of edge elements for computing electromagnetic fields. The main points against edge elements are listed below, for a full discussion I refer to my paper.
1. Contrary to what is stated by the advocates of edge elements these elements do allow spurious solutions, they even do soby definition.
2. The use of edge elements is very inefficient because of generating a unnecessarily large number of unknowns.
3. The representation of the field by edge elements has a poor condition as compared with Cartesian bases.
The good thing about edge-element methods is that, despite the fact that edge elements do allow spurious solutions, methods using them do not seem to suffer very much from it in practice and it is inefficiency that is their main disadvantage.
a.1 Why do edge elements not usually produce wrong results?
Edge elements do allowspurious, i.e. erroneous, solutions, there is no doubt about that, but methods using them do not usually produce errors, how can we understand that? The answer to this question is found by studying the properties of 1: the resulting system of equations together with 2: the solution vector and 3: the method of solution.
First we note that the system of equations is going to be solved iteratively. (The reason for this is, of course, that only iterative methods are practical for practical problems, but remember that the matrix is highly singular because of which a direct solution would be impossible anyhow.) In an iterative solution we have a starting solution vector and this vector is iteratively updated until a sufficiently accurate result is obtained. About this process there are two observations we can make here.
1. For the initial solution all elements of the solution vector are set to zero. Note that this initial 'solution' is 100% incorrect but because all coefficients in the expansion of the field are zero the normal continuity conditions across the interfaces between elements are satisfied exactly and no error is made in the normal fluxes. It are the normal fluxes where edge elements freely allow errors.
2. The iterative solution proceeds by manipulating the solution vector using the system matrix and possibly, for speeding up the process, an approximate inverse of this matrix. Neither the system matrix, nor its approximate inverse, contains any information about the normal continuity and consequently the errors in the normal fluxes remain unaltered i.e. zero.
b. The problem with nodal elements
The well known nodal elements, I prefer to refer to them as Cartesian elements, are very efficient and provide an optimum condition of the representation of the field. The problem with Cartesian elements is, however, that they have three unknowns (the three components of either the electric or magnetic field vector) at each node whereas four conditions that apply to the simplicial star of this node (3 for the curl and 1 for the divergence) are imposed upon them. Hence we have more independent equations than unknowns and problems must be expected, and they are in fact found.
2. What are the options?
Thinking about our
problem we can start
from two sides:
1. Use
edge elements and try to
remove the problems with normal fluxes by adding additional equations.
In this
way we may be able to exclude the possibility of spurious solutions but
we add
to the complexity of the large system matrix, i.e. we increase the
already
existing
inefficiency of the method. This must be a dead end.
2.
Use
Cartesian elements and add,
at each node, one degree of freedom making the number of conditions
applying at
each node (4) equal to the number of unknowns (4). The question then
becomes: 'How to add a meaningful additional degree of freedom
at each
node?'. It is that question that we want to answer below.
When
deciding how to add an
additional degree of freedom at a node we first note that the curl
equations
already provide sufficient information for computing the 3-component
vector. So
it is the modelling of the divergence that we need additional solution spacefor.
A hint to
solve our problem is
given by edge elements that do not have this type of problem, with edge
elements we have a much to
wide solution space since the normal fluxes across
interfaces are left free
to jump.
3. My proposal
Consider,
for
simplicity of the argument, a node in a regular rectangular
3-dimensional mesh
that consists of rectangular bricks, each of these bricks being
subdivided in
tetrahedra. Let us, again for simplicity of the argument,
assume that the
origin of our coordinate system coincides with this node and that we
aim at
computing the electric field strength E using
first order, i.e. linear, elements and that the medium is locally
homogeneous. We now know that we have 4 equations for the 3 unknowns, Ex,
Ey and Ez,at the centre of this simplicial
star.
For
solving our shortage of
unknowns we now propose to widen our solution space by splittingone of them, e.g. Exbut this choice is
arbitrary, in two unknowns viz Ex+,
applying to the
part of the simplicial star above the plane x=0,
and Ex-applying to the part below this plane. In this way we allow
for a
discontinuity in Ex across
the plane x=0, and
consequently a possible divergence across this 'interface'. Our fourth
equation, modelling the divergence, will now ensure that the divergence
in
the
simplicial star of this node is set to 0 (or to any
other value
required, an option not available with edge elements). In case we omit
this
fourth
equation we obtain a situation comparable with not prescribing normal
flux conditions at edge-element interfaces and the difference between
the
final
results Ex+and Ex- will
most probably be 0.
Applying this splitting to each node of the mesh we have four independent equations and four unknowns at each node throughout the domain of computation and reliable results can be obtained with an optimum efficiency. With linear elements we have 4, or perhaps 6, unknowns per node, with equivalent linear edge elements we would have tens of unknowns related to each node, the exact number depending on the mesh.
4. Some additional notes
1. Note
that the choice for
splitting up Ex is
arbitrary. Each of the other components
could have been chosen with the same result.
2. For
symmetry reasons one could
opt for splitting up all three components of the electric field
strength. This
gives rise to two more unknowns but we can solve this symmetrically by
adding
the two equations Ex+- Ex-= Ey+- Ey-=Ez+- Ez-.
This choice may be useful for the modelling
of the field near corners and points where interfaces intersect.
3. It is easily understood that our proposal can straightforwardly be generalized to the application at general meshes.
4. At arbitrarily oriented interfaces between different media the Cartesian basis of the representation of the field may locally be rotated so as to easily accommodate the local continuity conditions. Another option is to use edge elements along (intersecting) interfaces, but only there, like we did in the FEMAX packages. In this way we combine the best things of both worlds, simplicity and efficiency in homogeneous regions and flexbility along interfaces and their intersections. Note that the number of elements related to interfaces (2D) is, for practical problems, much smaller than the number of elements in the homogeneous subdomains (3D) because of which the efficiency of our method is not seriously reduced by the introduction of edge elements along interfaces.
5. Final conclusion
We have proposed a new and very efficient way of representing the field in a finite element method for the direct computation of electric or magnetic field strength distributions i.e. without the use of methods like vector potentials and/or methods based on the exclusive use of edge elements.
On a personal level:
When my health forced me to retire I was not anymore able to deal with the challenges I saw in the finite element modelling of electromagnetic fields. With the above thoughts I have closed the book.