Open access peer-reviewed chapter

The Finite Element Method Applied to the Magnetostatic and Magnetodynamic Problems

Written By

Dang Quoc Vuong and Bui Minh Dinh

Submitted: 15 June 2020 Reviewed: 24 August 2020 Published: 13 November 2020

DOI: 10.5772/intechopen.93696

From the Edited Volume

Finite Element Methods and Their Applications

Edited by Mahboub Baccouch

Chapter metrics overview

576 Chapter Downloads

View Full Metrics


Modelling of realistic electromagnetic problems is presented by partial differential equations (FDEs) that link the magnetic and electric fields and their sources. Thus, the direct application of the analytic method to realistic electromagnetic problems is challenging, especially when modeling structures with complex geometry and/or magnetic parts. In order to overcome this drawback, there are a lot of numerical techniques available (e.g. the finite element method or the finite difference method) for the resolution of these PDEs. Amongst these methods, the finite element method has become the most common technique for magnetostatic and magnetodynamic problems.


  • finite element method
  • magnetostastics
  • magnetodynamics
  • Maxwell’s equations
  • weak formulations

1. Introduction

Mathematical modeling of realistic problems in the framework of electromagnetics leads to a set of partial derivates equations that have to be solved on a domain with complex geometry associated with boundary conditions and initial conditions. This complexity makes any analytical approach unpracticable. In the past (until 1960), people used experimentation (very expensive, sometimes destructive) or analogic simulation (lack of generality) to solve these problems. Since 1970, the growth of computer capabilities makes the numerical simulation a tool that is more and more used by the people interested in solving these complex problems. When using the computer, the continuous problem is represented with a finite number of degrees of freedom (d.o.f.). The continuous problem is then replaced by a discrete problem. There are a lot of numerical techniques available. We will see that the most common ones can be derived from the same general principle of weighted residuals.

A continuous formulation of a problem cannot generally be solved analytically and some numerical methods have to be used in order to obtain quantitative information about the solution. The unknown functions of a continuous problem belong to continuous function spaces which are usually of infinite dimensions, that is, those functions are usually described by an infinite number of parameters. The basis of any numerical method is to discretize such a problem in order to obtain a similar discrete problem, characterized by a finite number of unknowns which are called degrees of freedom. This discretization process consists of replacing the considered continuous function spaces by some discrete function spaces, whose dimensions are finite, and which are usually subspaces of them. Those spaces are also called approximation spaces and their elements are called approximation functions.

The function spaces are defined in a particular studied domain. If this one is discretized, that is, if it is defined as the union of geometric elements of simple shapes, and if the discrete function spaces are built in such a way that their functions are piecewise defined, then the approximation numerical method is called the finite element method (FEM). It is this kind of method we are interested in. We can thus see that the finite element method necessitates a double discretization: a discretization of some function spaces and a discretization of the studied geometric domain, which leads to a mesh.

Weak formulations are well adapted to the finite element method, which will appear in the following. Such formulations make use of several kinds of Green formulas.


2. Numerical technique

2.1 The Laplacian problem

The formalism used in the case of a Laplacian problem is sufficiently simple to be very understandable without lack of generality. The description of a Laplacian problem is presented now. Let us consider a bounded domain Ω and its boundary Γ=ΓhΓe (Figure 1).

Figure 1.

Studied domain Ω and its boundary Γ = ΓhΓe.

The Laplace equation has to be solved in Ω [1, 2, 3]:


where u is the unknown field defined at each point x (x, y, z) of the studied domain. The associated boundary conditions are respectively Dirichlet and Neumann conditions, that is


This diffusion equation describes a wide range of physical phenomena. The next table shows some of these phenomena.

Problemsuon Γhon Γe
ThermostaticsT (temperature)T fixedthermal flux fixed
ElectrostaticsV (voltage)V fixed (electrode)fixed flux of electrical displacement
Perfect fluidsψ (flow function)
I (current function)
MagnetostaticsBr (reduced potential)fixed magnetic flux densityfixed tangential magnetic field

A natural way to discretize the problem is to impose the error on the equation and on the boundary conditions weighted by a trial function w to be equal to zero, that is


Equations (1–3) are then meanly solved, the sense of the mean being the principle of the numerical method. In fact, the numerical method used (F.D.M, F.E.M or B.E.M) are directly related to the chosen trial functions.

2.2 Green formulas

The following notations are used for integration of products of scalar or vector fields over a volume Ω or on a surface Γ, where L2 and L2 are the spaces of square-summable scalar and vector functions [2, 3]:


A first relation of vectorial analysis


integrated in the domain Ω, after applying the divergence theorem, gives the Green formula said of kind grad-div in Ω, that is

ugradv+divuv=<v,n·u>Γ, uH1Ω,vH1ΩE5

where H1ΩandvH1Ω are function spaces built for scalar and vector fields, respectively.

Another relation of vectorial analysis


integrated in the domain Ω, after applying the divergence theorem, gives the Green formula said of kind curl-curl in Ω, that is

ucurlvcurluv=<u×n,v>Γ, .u,vH1ΩE7

Note that the surface integral term appearing in this last formula can take the following similar forms:


It is possible to define a generalized Green formula by


where L and L* are differential operators of order n which act respectively on functions u and v defined in Ω¯, with Ω¯ = ΩΓ; Q is a bi-linear function of u and v. The operator L* is called the dual operator of L. It can easily be seen that formulas (6) and (7) are particular cases of (8).

2.3 Weak formulations

Consider a partial differential problem of the form [4].


where L is a differential operator of order n, B is an operator which defines a boundary condition, f and g are functions respectively defined in Ω and on its boundary Γ, and u is an unknown function from a function space U and defined in Ω¯, that is, u ∈ U(Ω¯). Note that f can eventually depend on u.

Problems (9 and 10) constitute what is called a classical formulation, or strong formulation. A function u ∈ U(Ω¯) which verifies this problem is called a classical solution, or strong solution. Particularly, as L is of order n, the function u has to be n–1 times continuously differentiable, that is, u ∈ Cn–1(Ω).

A weak formulation of problem (9) is defined as having the generalized form.


where L* is the dual operator of L, defined by the generalized Green formula (8), Qg is a linear form in v which depend on g, and the space V (Ω) is a space of test functions which has to be defined according to the operator L* and particularly according to the boundary condition (9 and 10). A function u which satisfies this equation for any test function v ∈ V (Ω) is called a weak solution.

The generalized Green formula (8) can be applied to formulation (11) in order to get L instead of L*, which usually consists of performing an integration by parts. It is then possible to find again, thanks to a judicious choice of test functions, the equations and relations of the classical formulation of the problem, that is, Eq. (9) and boundary condition (11).

It is often easy to check that a classical solution is also a weak solution. Nevertheless, it is not always straightforward that a weak solution is also a classical solution because it has to be regular enough in order to be defined at the classic sense.

One mathematical advantage of weak formulations is that they usually enable to prove the existence of a solution easier than classical formulations do. The solution has then to be proved to be regular enough to be also a classical solution. Another advantage of weak formulations is that they are well adapted to a discretization using finite elements and then to a numerical solution, which is not the case with classical formulations.

In some cases, it is possible to define a minimization problem similar to the weak formulation (11).

2.4 A weak formulation for the magnetodynamic problem

In order to illustrate the notion of weak formulation, consider the magnetodynamic problem, limited to the domain Ω, with boundary ∂Ω = Γ = ΓhΓe (Figure 2), whose equations and material relations are written in Euclidean space R3[5, 6].


with behavior relations of materials.

Figure 2.

Studied domain Ω and its boundary.


where j is called the imaginary unit, b is the magnetic induction (T), e is the electric field (V/m), js is the current density (A/m2), h is the magnetic field (A/m), μ is the relative permeability and σ is the electric conductivity (S/m). From the Eq. (13), the field b can be obtained from a magnetic vector potential ai via the term:


Combining (15 and 16) into (14), one has curl (e+ta)=0, that leads to the presentation of an electric scalar potential ν through e=tagradυ.

By starting from the Ampere’s law (12), the weak form of magnetic vector potential is written as [4, 6].


Combining the magnetic vector potential a and the electrical field e, one has


where He0curlΩis a function space defined on Ω containing the basis functions for a as well as for the test function a (at the discrete level, this space is defined by edge FEs; notations (·, ·) and < ·, · > are respectively a volume integral in and a surface integral of the product of their vector field arguments.

2.5 A weak formulation for the magnetostatic problem

The magnetostatic problem is considered as a simplification of the magnetodynamic formulation where all time dependent phenomena are neglected. In a same way, by starting from the Ampere’s law (12), this initial form of the problem is its classical formulation.

Consider the Green formula of type grad-div in Ω (5) applied to the field b and to a scalar field ϕ’ to be defined, that is [6, 7, 8]


If the space Φ(Ω) is defined such as


then the last term of Eq. (20) is reduced to < n.b, ϕ’ > Γe and is equal to zero if condition (15) is taken into account. Moreover, the second term of this equation is equal to zero because of Eq. (16). Eq. (20) can then be reduced to.


This last form is called a weak formulation of the problem. It has been established starting from a Green formula but it can be considered now as an a priori posed form whose enclosed information can be deduced.

In fact, weak formulation (22) contains both Eq. (12) and boundarycondition (15). Indeed, by applying the Green formula of type grad-div to it, we get.


This equation is verified for any test function ϕ’ ∈ Φ(Ω) and thus, particularly, for any function ϕ’ whose value is equal to zero on Γ, that is, ϕ’ ∈ Φ0 (Ω) because Φ0(Ω) ⊂ Φ(Ω). Therefore, it comes that (div b, ϕ’) = 0 for any function ϕ’ of this kind and, consequently, that div b = 0 in Ω, that is, Eq. (12) is satisfied. Then, Eq. (23) is reduced to nbϕ>Γ= 0, and by considering now all the functions ϕ’ ∈ Φ(Ω) without any restriction, that is, which can vary freely on Γe, it comes that nbΓe= 0, that is, that condition (13) is satisfied.

It is possible to obtain more information from the weak formulation, particularly as far as the transmission conditions on surfaces inside Ω are concerned. Consider for that two subdomains Ω1 and Ω2 of Ω separated by an interface Σ (Figure 3) [7].

Figure 3.

Interface Σ between two subdomains Ω1 and Ω2.

Let us apply the Green formula of type grad-div (5) to the fields b and ϕ’ successively in both subdomains Ω1 and Ω2. By taking into account that div b = 0 in Ω, and thus particularly in Ω1 and Ω2, then by summing the obtained relations, we get the relation [6, 7].


where b1 and b2 represent the field b on both sides of Σ in the respective domains Ω1 and Ω2. Considering the test functions ϕ’ whose support is Ω1∪Ω2 and which are equal to zero on _(Ω1∪Ω2), it remains from (24) the well known transmission condition n.(b1b2)∣Σ = 0. Note that the first term of (24) vanishes thanks to Eq. (22) indeed, the domain of integration Ω1 ∪ Ω2 can be extended to Ω thanks to the chosen test functions.

The way to establish a weak formulation of Eq. (13) has been described and the richness of the information enclosed in such a formulation has been shown up. Using a similar procedure, a weak formulation associated with Eq. (12) can be established, but we will proceed differently in order to keep some classical equations.

If the field h is decomposed into a given source component hs, such as curl hs = j, and a reaction component hr, then curl hr = 0 and hr is therefore of the form hr = − grad ϕ (if Ω is simply connected). This consists of satisfying Eq. (15) classically. Taking into account the behavior law (15), we can write b = μ (hs − grad ϕ) and put this last expression in (24) to obtain.


This formulation contains the whole problem (12 and 13). The potential ϕ is the unknown and all the other fields can be deduced from ϕ thanks to the equations which have been kept on a classical form. It appears that the potential ϕ belongs to the same space of the test functions or at least to a space Φr (Ω) which is parallel to it, that is, where the boundary condition relative to ϕ on Γh is not necessarily homogeneous, that is, ϕ∣Γh = constant. Note that this boundary condition on Γh is implicitly taken into account in the space Φ(Ω).

Weak formulation (25) can be considered as a system of an infinite number of equations with an infinite number of unknowns. It will be seen in the following how such a problem can be approximated to lead to a numerical solution. This approximation will constitute the phase of discretization.

A similar minimization problem

It is possible to define a minimization problem associated with (25). For that, let us define the functional [2, 3].


and let us pose the following minimization problem:


The physical materials are considered having linear magnetic behavior, but the following can be generalized easily for nonlinear materials.

By stationarizing functional (25) in relation to ϕ, it can be verified that (25) is obtained. It can also be verified that the solution ϕ of (25) minimizes this functional. Indeed, let us suppose that ϕ ∈ Φr(Ω) is solution of (25) and let us consider any function ψ ∈ Φr(Ω); then let us pose η = ψ − ϕ, which implies η ∈ Φ(Ω); we have


and thus.


As the second term of this sum is positive or equal to zero and the third term is equal to zero, because of (25), it comes that W(ψ) and W(ϕ).

Formulations (25) and (27) are then similar. Note that W(ϕ) is the magnetic coenergy and that the problem actually consists of minimizing this coenergy.

If the continuous function spaces are replaced by discrete spaces, and if the considered test functions are limited to these spaces, then the information inside a weak formulation will only be satisfied approximately, or weakly.

The basis of the discretization of weak formulations can be illustrated for the above magnetostatic problem, whose weak formulation is (25), that is


with ϕ ∈ Φ(Ω). The space Φ(Ω) has to be replaced by a discrete space Φh(Ω) which is a subset of it, that is, Φh(Ω) ⊂ Φ(Ω). This space has a finite dimension, denoted N, and can then be defined by N linearly independent base functions. The principle is then to look for the function ϕ in Φh (Ω), which consists of determining N unknown parameters. This function will be only an approximation of the exact solution ϕ ∈ Φ (Ω). The more the functions of Φh (Ω) approximate well those of Φ (Ω), the higher the quality of the approximation is. Each test function ϕ’ will lead to an equation of the form (28) and, as the number of equations and unknowns has to be the same, N linearly independent test functions have to be chosen. This choice can be made on the base functions of Φh (Ω) and the method is called the Galerkine method. Such base functions are defined thanks to finite elements.


3. Finite elements

3.1 Definition of a finite element

A finite element is defined by the three element set (K, PK, ΣK) where [2, 6, 7]:

  • K is a domain of space called a geometric element (usually of simple shape, that is, a tetrahedron, a hexahedron or a prism);

  • PK is a function space of dimension nK, defined in K with base functions;

  • ΣK is a set of nK degrees of freedom represented by nK linear functionals ϕi, 1 ≤ i ≤ nK, defined in space PK;

moreover, any function u ∈ PK must be defined uniquely by the degrees of freedom of ΣK, which defines the unisolvance of the finite element (K, PK, ΣK).

The role of a finite element is to interpolate a field in a function space of finite dimension. Several finite elements can be defined on the same geometric element and then, under certain conditions, can form mixed finite elements. Figure 4 shows the various spaces which occur in the definition of a finite element; the definition of the subspace of points κ ⊂ K is actually associated with the definition of the functionals.

Figure 4.

Spaces associated with a finite element (K, PK, ΣK) [6].

For the most commonly used finite elements, the degrees of freedom are associated with nodes of K and the functionals ϕi are reduced to functions of the coordinates in K; these elements are called nodal finite elements. Nevertheless, the above definition is more general thanks to the freedom let in the choice of the functionals. It will be shown that these can be, in addition to nodal values, integrals along segments, on surfaces or in volumes; the subspace of points κ ⊂ K (Figure 4) is then respectively a point, a segment, a surface or a volume.

3.2 Unisolvant finite element

The finite element (K, PK, ΣK) is unisolvant if [6].


In this case, for any function u regular enough, one can define a unique interpolation uK, called PK-interpolant, such as.


The set ΣK is said PK - unisolvant.


Each function p ∈ PK can be written as a linear combination of functions of a base of PK, denoted {pi, 1 ≤ i ≤ nK}, that is


where the pi, 1 ≤ i ≤ nK, are called base functions.

As the functionals ϕj, 1 ≤ j ≤ nK, are linear, we have


And, as ϕj(p) = 0, 1 ≤ j ≤ nK, leads to p ≡ 0, the determinant of the matrix Φ (Φji = ϕj(pj), 1 ≤ i, j ≤ nK) is not equal to zero; indeed the solution of the corresponding system must be identically equal to zero (i.e., ai = 0, 1 ≤ i ≤ nK). Consequently, the system


has a unique solution (aj, 1 ≤ i ≤ nK).

3.3 Degrees of freedom

The interpolation of a function u, in the space PK and in K, is given by the expression [4]


where the nK coefficients aj associated with the base functions pj ∈ PK can be determined thanks to relations (26), that is, thanks to the solution of the linear system.


provided that the function u is sufficiently regular for the ϕj(u), 1 ≤ j ≤ nK, to exist.

This solution is simplified to the maximum if we define the functionals so that


where δi,j is the Kronecker symbol, that is


The matrix of the system is then the unit matrix and the solution is


In this case, the interpolation uK ∈ PK is expressed by


where the coefficients ϕj(u) = ϕj(uK), 1 ≤ j ≤ nK, are called degrees of freedom.

3.4 Finite element spaces

A finite element space Xh can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh Mh of the domain Ω as well as the knowledge of the finite element (K, PK, ΣK) associated with each domain K ∈ Mh [6]

Given a function u defined in Ω, regular enough, its interpolant uh ∈ Xh is uniquely defined such as [6]:

  • The restriction uhK, that is, the form of uh in the geometric element K, belongs to the space PK;

  • The restriction Finite element spaces

A finite element space Xh can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh Mh of the domain Ω as well as the knowledge of the finite element (K, PK, ΣK) associated with each domain K∈ Mh

Given a function u defined in Ω, regular enough, its interpolant uh ∈ Xh is uniquely defined such as [6]:

  • The restriction uhK, that is, the form of uh in the geometric element K, belongs to the space PK;

  • The restriction uhK is entirely determined by the knowledge of the set of values ΣK(u) of the degrees of freedom of the function u - this is a consequence of the unisolvance;

Some continuous conditions have to be ensured across the interfaces between geometric elements, which is the property of conformity.

A mesh Mh of the studied domain Ω is defined as a collection of geometric elements which have in common either a facet, or an edge, or a node, or nothing (Figure 5). The elements cannot overlap each other.

Figure 5.

Mesh of a part of a two-dimensional domain Ω.

The finite element space Xh has a finite dimension, denoted Dh. It can be characterized by a set of degrees of freedom Σh linked up to the sets ΣK,∀K∈ Mh, that is


It is also possible to define the base functions ph,j, 1 ≤ i ≤ Dh, of the space Xh from the base functions of the spaces PK, ∀ K∈ Mh. Those have to verify the relations


similar to relations (32). They are actually piecewise defined and their supports are as “small” as possible, that is, are constituted by a limited number of geometric elements.

Then, with any function u regular enough so that the degrees of freedom ϕh,j(u), 1 ≤ j ≤ Dh, are well defined, it can be associated a function uh, called Xh-interpolant, defined by


4. Construction of a sequence of finite element spaces

4.1 Geometric elements

A mesh of a domain is considered which is built with a collection of geometric elements which can be tetrahedra (4 nodes), hexahedra (8 nodes) and prisms (6 nodes) (Figure 6) [4, 6, 9].

Figure 6.

Collection of different geometric elements [6].

These elements are called volumes and their vertices represent nodes. The sets of nodes, edges, facets and volumes of this mesh are denoted by N, E, F and V, respectively. Their sizes are #N, #E, #F and #V.

The i-th node of the mesh is denoted by ni or {i}. The edges and facets can be defined with ordered sets of nodes. An edge is denoted by eij or {i, j}, a triangular facet by fijk or {i, j, k}, and a quadrangular facet by fijkl or {i, j, k, l}. These geometric entities are shown in Figure 7.

Figure 7.

Geometric entities: Node, edge and facets (i, j, k, l ∈ N) [6].

4.2 Base functions of spaces Si

Consider the function pi(x) of coordinates of point x and relative to node ni, which is equal to 1 at this node, varies continuously in geometric elements having this node in common, and becomes equal to 0 in other elements without discontinuity (Figure 6). This function is nothing else than the base function, relative to node ni, of the function space of nodal finite elements built on the considered geometric elements. The function subspaces associated with each of the finite elements have respective dimensions 4, 8 or 6, for tetrahedra, hexahedra and prisms [4, 6, 9].

With node ni = {i}, is associated the function


The finite dimensional space generated by all sni’s is denoted by S0.

With edge eij = {i, j}, is associated the vector field


where NF,mn¯ is the set of nodes which belong to the facet of the geometrical element including evaluation point x, and including node m but not node n; such a facet is uniquely defined for three-edge-per-node elements. Its determination is shown in Figure 8, where either a triangular or a quadrangular facet is involved, and where shown edges belong to the geometric element including point x. Directions of dotted edges can be modified in order to schematize either a tetrahedron, a hexahedron or a prism. The defined set of nodes comes into view as being either {{m}, {o}, {p}} or {{m}, {o}, {p}, {q}}, respectively. The set NF,mn¯ depends on point x, thus on elements. Particularly, it is empty (no node) in elements which have not edge {m, n} in common. Consequently, field seij is zero in all the elements non adjacent to edge eij.

Figure 8.

Determination of the facet associated with NF,mn¯ [6].

The vector field space generated by all se is denoted by S1.

With facet f = fijk = {i, j, k} = {q1, q2, q3} or f = fijkl = {i, j, k, l} = {q1, q2, q3, q4}, is associated the vector field


where #Nf is the number of nodes of facet f, af = 2 if #Nf = 3, af = 1 if #Nf = 4, and the list of qi’s is made circular by setting q0 ≡ q#Nf and q#Nf + 1 ≡ q1. Field sf is zero in all the elements non adjacent to facet f.

Vector fields sfijk(l)'s generate the space S2.

With volume v, is associated the function sv, equal to 1/vol(v) on v and 0 elsewhere. The space S3 is generated by these functions.

Some developments give the following results: sni is equal to 1 at node ni, and to 0 at other nodes; the circulation of seij is equal to 1 along edge eij, and to 0 along other edges; the flux of sfijk(l) is equal to 1 across facet sfijk(l), and to 0 across other facets; and the volume integration of sv is equal to 1 over volume v, and to 0 over other volumes; that is


where δij = 1 if i = j and δij = 0 if i ≠ j.

These properties show up various kinds of functionals and involve that functions sn, se, sf, sv form bases for the spaces they generate. They are then called nodal, edge, facet and volume base functions. The associated finite elements are called nodal, edge, facet and volume finite elements.

4.3 Geometric interpretation of edge and facet functions

A geometric interpretation of edge and facet functions may be helpful to verify some of their properties. The vector field [4, 6, 9]


involved in both expressions (31) and (32), should be analyzed at first. The continuous scalar field,


has the characteristic of being equal to 1 at every point on the facet associated with NF,mn¯. This is a property of the nodal base functions. Therefore, vector field (42) is orthogonal to this facet at every point on it (Figure 9).

Figure 9.

Geometric interpretation of the edge function se(35) [6].

The vector field which is the product of pm and (35),


is considered now. This field is said to be associated with edge {m, n}. As far as the function pm is concerned, it is equal to 0 on all the edges of the geometric element including point x, except those which are incident to node {m}. Therefore, the circulation of (43) is equal to 0 along all the edges except emn; field (43) is either simply equal to zero on them, or orthogonal to them (Figure 9). The combination of two fields of form (43) associated with edges {j, i} and {i, j}, as in (35), leads to a vector field which has the same properties as (43) (Figure 9), and has consequently the announced properties of seij. The fact that its circulation along edge eij is equal to 1 needs some calculation to be proved.

The vector field


which appears in expression (35) of sf, is considered now. Both gradients in (44) are shown in Figure 10. Each one is orthogonal to its associated facet and, therefore, their cross product (i.e., a × b in Figure 10) is parallel to both these facets.

Figure 10.

Vector field a × b involved in sf(33) [6].

The flux of this cross product, and in consequence the one of (44), is then equal to 0 across these facets. The term pqc in (44) enables the flux of (44) to be equal to zero across all other facets except facet f. The summation in (44) keeps the same property. The flux of sf across facet f is then the only one to differ from zero (Figure 11).

Figure 11.

Geometric interpretation of the facet function sf(36) [6].

4.4 Degrees of freedom

The expression of a field in the base of a space Si –S0 or S3 for a scalar field, S1 or S2 for a vector field– gives scalar coefficients, called degrees of freedom. Fields ϕ ∈ S0, h ∈ S1, j ∈ S2 and ρ ∈ S3 can be expressed as [4, 6, 9]


The degrees of freedom ϕn, he, jf and ρv are thus, respectively, values at nodes, circulations along edges, fluxes across facets or volume integrals, of the associated fields. This is a consequence of the base functions. The associated linear functionals, as mentioned in the definition of finite elements, are thus respectively pointwise evaluations, line, surface and volume integrals.

4.5 Continuity of base functions across facets

It can be proved that the function sn is continuous across facets. The same holds true for the tangential component of se and for the normal component of sf. As for function sv, it is discontinuous. This property, called conformity, allows to take exactly into account interface conditions for fields used in the modeling of physical problems. For example, in electromagnetic problems, vector fields of S1 can represent vector fields like magnetic field h or electric field e whose tangential components are continuous across interfaces between materials, and those of S2 can represent fields like induction field b or current density field j whose normal components are continuous across interfaces between these materials.

4.6 Spaces Si form a sequence

The notion of incidence is first defined [4, 6, 9]:

The incidence of node n in edge e, denoted by i (n, e), is equal to 1 if n is the extremity of e, −1 if n is the origin of e, and 0 if n does not belong to e.

Next, the incidence of edge e in facet f is denoted by i(e, f). If e belongs to f, and if the ordered set of nodes of e appears as a direct subset in the circular set of nodes of f, then it is equal to 1. It is equal to −1 in the case of an inverse subset. If e does not belong to f, it is equal to 0.

Finally, the incidence of facet f in volume v is denoted by i(f, v). If f belongs to v, and if the normal to f, whose direction is given by the ordered set of nodes of f (right-hand rule), is outer to v, then it is equal to 1. It is equal to −1 in the case of an inner normal. If f does not belong to v, it is equal to 0.

Thanks to this notion, the following equalities can be proved,


The following inclusions are then verified,


Therefore, the spaces Si, i = 0 to 3, form a sequence, that can be schematized by the diagram in Figure 12.

Figure 12.

The sequence of spaces Si.

These spaces can then constitute approximation spaces for some continuous spaces Fi, i = 0 to 3, which contain scalar and vector fields associated with electromagnetic fields. The associated finite elements can then be called mixed elements.

All the established properties of base functions are valid for any collection of considered geometric elements, that is, for any mixing of tetrahedra, hexahedra and prisms.


5. Practical information about finite elements

5.1 Isoparametric elements

An isoparametric element is a finite element whose nodal base functions, which enable the interpolation of scalar fields, are also used to parametrize the associated geometric element. The base functions are usually piecewise defined, in each of the geometric elements which cover the studied domain, and some continuity conditions have to be satisfied at the interfaces between elements. Then, there will be no discontinuity of the interpolated scalar fields, nor of the coordinates after transformation from the reference elements towards the real ones. Such base functions are said to be conformal.

Consider a nodal finite element (K, PK, ΣK). If NK is the set of nodes of K, whose coordinates are xi, i ∈ N, and if the pi(u), i ∈ NK, are its base functions expressed in the coordinates u of the reference element Kr associated with K, then the parametrization of K (i.e., x = x(u)) is given by [6]


where x∈K, u∈Kr; this element is isoparametric.

5.2 Reference elements

We define here the reference elements which are associated with the considered geometric elements, that is, with tetrahedra, hexahedra and prisms. Nodal, edge, facet and volume finite elements are defined in these geometric elements.

5.2.1 Reference tetrahedron of type I

The reference tetrahedron of type I is an element with 4 nodes whose coordinates are given in Figure 13. The associated geometric entities, as well as their notation, are shown in Figure 13. The nodal and edge base functions of this element are given in Tables 1 and 2. Table 3 shows the notation of facets. The incidence matrices are given by (53), (54) and (55) (Figure 14).

Figure 13.

Reference tetrahedron of type I [6].

Nodal base function
pi (u, v, w) = si (u, v, w)
11 − u − v − w

Table 1.

Nodal base functions of the tetrahedron of type I.

Edgee = {i, j}se(u), u = (u, v, w)
1121 − v –wuu
213v1 − u –wv
314ww1 − u –v
423– vu0
524– w0u
6340– wv

Table 2.

Notation of the edges of the tetrahedron of type I and associated edge base functions (se).

Facetf = {i, j, k}

Table 3.

Notation of the facets of the tetrahedron of type I.

Figure 14.

Geometric entities defined on a tetrahedron of type I [6].

Edge-node incidence matrix


Facet-edge incidence matrix


Volume-facet incidence matrix


5.2.2 Reference hexahedron of type I

Figure 15.

Reference hexahedron of type I [6].

Figure 16.

Geometric entities defined on a hexahedron of type I [4].

Nodal base function
pi (u, v, w) = si (u, v, w)
1(1 − u) (1 − v) (1 − w) / 8
2(1 + u) (1 − v) (1 − w) / 8
3(1 + u) (1 + v) (1 − w) / 8
4(1 − u) (1 + v) (1 − w) / 8
5(1 − u) (1 − v) (1 + w) / 8
6(1 + u) (1 − v) (1 + w) / 8
7(1 + u) (1 + v) (1 + w) / 8
8(1 − u) (1 + v) (1 + w) / 8

Table 4.

Nodal base functions of the hexahedron of type I.

Edgee = {i, j}se(u), u = (u, v, w)
112(1 − v) (1 − w) / 800
2140(1 − u) (1 − w) / 80
31500(1 − u) (1 − v) / 8
4230(1 + u) (1 − w) / 80
52600(1 + u) (1 − v) / 8
634–(1 + v) (1 − w) / 800
73700(1 + u) (1 + v) / 8
84800(1 − u) (1 + v) / 8
956(1 − v) (1 + w) / 800
10580(1 − u) (1 + w) / 80
11670(1 + u) (1 + w) / 80
1278–(1 + v) (1 + w) / 800

Table 5.

Notation of the edges of the hexahedron of type I and associated edge base functions (se).

Facetf = {i, j, k, l}

Table 6.

Notation of the facets of the hexahedron of type I.

The reference hexahedron of type I is an element with 8 nodes whose coordinates are given in Figure 15. The associated geometric entities, as well as their notation, are shown in Figure 16. The nodal and edge base functions of this element are given in Tables 4 and 5. Table 6 shows the notation of facets. The incidence matrices are given by (56), (57) and (58).

Edge-node incidence matrix


Facet-edge incidence matrix


Volume-facet incidence matrix


5.2.3 Reference prism of type I

The reference prism of type I is an element with 6 nodes whose coordinates are given in Figure 17. The associated geometric entities, as well as their notation, are shown in Figure 18. The nodal and edge base functions of this element are given in Tables 7 and 8. Table 9 shows the notation of facets. The incidence matrices are given by (59), (60) and (61).

Figure 17.

Reference prism of type I [6].

Figure 18.

Geometric entities defined on a prism of type I [6].

Nodal base function
pi (u, v, w) = si (u, v, w)
1(1 − u − v) (1 − w) / 2
2u (1 − w) / 2
3v (1 − w) / 2
4(1 − u − v) (1 + w) / 2
5u (1 + w) / 2
6v (1 + w) / 2

Table 7.

Nodal base functions of the prism of type I.

Edgee = {i, j}se(u), u = (u, v, w)
112(1 − v) (1 − w) / 2u (1 − w) / 20
213v (1 − w) / 2(1 − u) (1 − w) / 20
31400(1 − u − v) / 2
423– v (1 − w) / 2u (1 − w) / 20
52500u / 2
63600v / 2
745(1 − v) (1 + w) / 2u (1 + w) / 20
846v (1 + w) / 2(1 − u) (1 + w) / 20
956– v (1 + w) / 2u (1 + w) / 20

Table 8.

Notation of the edges of the prism of type I and associated edge base functions (se).

Facettef = {i, j, k (, l)}

Table 9.

Notation of the facets of the prism of type I.

Edge-node incidence matrix


Facet-edge incidence matrix


Volume-facet incidence matrix


6. Applications

The practical test problem is a 3-D model based on the benchmark problem 19 of the TEAM workshop including a stranded inductor (coil) and an aluminum plate (Figure 19) [10].

Figure 19.

Modeling of TEAM problem 7: Coil and conducting plate [10].

The coil is excited by a sinusoidal current which generates the distribution of time varying magnetic fields around the coil (Figure 20). The relative permeability and electric conductivity of the plate are μr,plate=1,σr,plate=35.26MS/m, respectively. The source of the magnetic field is a sinusoidal current with the maximum ampere turn being 2742AT. The problem is tested with two cases of frequencies of the 50 Hz and 200 Hz.

Figure 20.

The 3-D mesh model with edge elements of the coil and conducting plate, and the limited boundary [4] (left), and distribution of magnetic flux density generated by the excited sinusoidal current in the coil, with μr,plate=1,σr,plate=35.26MSmand f = 50 Hz.

The 3-D dimensional mesh with edge elements is depicted in Figure 21 (left). The distribution of magnetic flux density generated by the excited electric current in the coil is pointed out in Figure 21 (right). The computed results on the of the z-component of the magnetic flux density along the lines A1-B1 and A2-B2 (Figure 19) is checked to be close to the measured results for different frequencies of exciting currents (already proposed by authors in [10]) are shown in Figure 21. The mean errors between calculated and measured methods [10] on the magnetic flux density are lower than 10%.

Figure 21.

The comparison of the calculated results with the measured results on magnetic flux densities at y = 72 mm, with μr,plate=1,σr,plate=35.26MSmand different frequencies [4].

The y component of the varying of the eddy current losses with different frequencies (50 Hz and 200 Hz) along the lines A3-B3 and A4-B4 (Figure 19) is shown in Figure 22. The computed results are also compared with the measured results as well [7]. The obtained results from the theory modeling are quite similar as what measured from the measurements. The maximum error near the end of the conductor plate on the eddy currents between two methods are below 20% for both cases (50 Hz and 200 Hz).

Figure 22.

The comparison of the calculated results with the measured results at z = 19 mm, with μr,plate=1,σr,plate=35.26MSmand different frequencies [4].


7. Conclusions

In the 3D computation of the magnetic flux density and eddy current, thanks to the set of Maxwell’s equations, it has been successfully developed for two weak formulations, where the discretization of the fields is performed by Whitney edge elements [2, 3, 8]: magnetostatic formulation and magnetodynamic formulation. The developments of the method is validated on the actual problem (TEAM problem 7) [10]. The numerical error between simlated and measured results on the magnetic flux densities and eddy current is lower than 10%. This is also proved that there is a very good validation between two methods. The results have been achieved by a detailed study of the magnetodynamic formulation.


  1. 1. Gerard Meunier, et al,. “The Finite Element Method for Electromagnetic Modeling”, John Wiley & Sons Inc, 2008
  2. 2. A. Bosavit, “Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism”, IEEE Proceedings, Vol. 135, Pt. A, no. 8, pp. 493–499, 1998
  3. 3. A. Bosavit, “Electromagnetisme en vue de la modelisation”, Collection Mathematiques et Applications, Springer-Verlag, 1991
  4. 4. Vuong Dang Quoc and Christophe Geuzaine “Using edge elements for modeling of 3-D magnetodynamic problem via a subproblem method”, Science and Technology Development Journal. ; 23(1): 439–445
  5. 5. S. Koruglu, P. Sergeant, R.V. Sabarieqo, Vuong. Q. Dang, M. De Wulf “Influence of contact resistance on shielding efficiency of shielding gutters for high-voltage cables,” IET Electric Power Applications, Vol.5, No.9, (2011), pp. 715–720
  6. 6. Quoc VD. Modeling of Electromagnetic Systems by Coupling of Subproblems with Application to Thin Shell Models [Ph. D Thesis]. 2013. Belgium: University of Liege
  7. 7. Vuong Q. Dang, P. Dular R.V. Sabariego, L. Krähenbühl, C. Geuzaine, “Subproblem approach for modelding multiply connected thin regions with an h-conformal magnetodynamic finite element formulation”, in EPJ AP., vol. 63, no. 1, 2013
  8. 8. P. Dular, W. Legros, A. Nicolet, “Coupling of local and global quantities in various finite element formulations and its application to electrostatics, magnetostatics and magnetodynamics”, IEEE Transactions on Magnetics, vol. 34, no. 5, pp. 3078–3081, 1998
  9. 9. P. Dular, J.-Y.Hody, A. Nicolet, A. Genon, W. Legros, “Mixed finite elements associated with a collection of tetrahedra, hexahedra and prisms”, IEEE Transactions on Magnetics, vol. 30, no. 5, pp. 2980–2983, 1994
  10. 10. Kovacs G, Kuczmann M. Solution of the TEAM workshop problem No.7 by the finite Element Method. In: Approved by the International Compumag Society Board at Compumag. 2011

Written By

Dang Quoc Vuong and Bui Minh Dinh

Submitted: 15 June 2020 Reviewed: 24 August 2020 Published: 13 November 2020