The Finite Element Method Applied to the Magnetostatic and Magnetodynamic Problems

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.


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.

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).
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  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.

Green formulas
The following notations are used for integration of products of scalar or vector fields over a volume Ω or on a surface Γ, where L 2 and L 2 are the spaces of squaresummable 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 where H 1 Ω ð Þ and v ∈ H 1 Ω ð Þ 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 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).

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).

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  3 [5,6].
with behavior relations of materials.
where j is called the imaginary unit, b is the magnetic induction (T), e is the electric field (V/m), j s 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 a i via the term: Combining (15 and 16) into (14), one has curl (e þ ∂ t aÞ ¼ 0, that leads to the presentation of an electric scalar potential ν through e ¼ À∂ t a À grad υ. 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 H 0 e curl, Ω ð Þis a function space defined on Ω containing the basis functions for a as well as for the test function a 0 (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.

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] 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 boundary condition (15). Indeed, by applying the Green formula of type grad-div to it, we get.
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].
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 b 1 and b 2 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.(b 1 -b 2 )|Σ = 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 h s , such as curl h s = j, and a reaction component h r , then curl h r = 0 and h r is therefore of the form h r = À 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.

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.
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.
In this case, for any function u regular enough, one can define a unique interpolation uK, called P K -interpolant, such as.
The set Σ K is said P K -unisolvant. Proof: Each function p ∈ P K can be written as a linear combination of functions of a base of P K , denoted {pi, 1 ≤ i ≤ n K }, that is p ¼ where the pi, 1 ≤ i ≤ nK, are called base functions.
As the functionals ϕj, 1 ≤ j ≤ n K , are linear, we have And, as ϕj(p) = 0, 1 ≤ j ≤ n K , leads to p 0, the determinant of the matrix Φ (Φ ji = ϕ j (p j ), 1 ≤ i, j ≤ n K ) is not equal to zero; indeed the solution of the corresponding system must be identically equal to zero (i.e., ai = 0, 1 ≤ i ≤ n K ). Consequently, the system has a unique solution (a j , 1 ≤ i ≤ n K ).

Degrees of freedom
The interpolation of a function u, in the space P K and in K, is given by the expression [4] u K ¼ where the n K coefficients a j associated with the base functions p j ∈ P K 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 ≤ n K , 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 u K ∈ P K is expressed by where the coefficients ϕj(u) = ϕj(u K ), 1 ≤ j ≤ n K , are called degrees of freedom.

Finite element spaces
A finite element space X h can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh M h of the domain Ω as well as the knowledge of the finite element (K, P K , Σ K ) associated with each domain K ∈ M h [6] Given a function u defined in Ω, regular enough, its interpolant u h ∈ X h is uniquely defined such as [6]: • The restriction u h j K , that is, the form of u h in the geometric element K, belongs to the space P K ; • The restriction Finite element spaces A finite element space X h can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh M h of the domain Ω as well as the knowledge of the finite element (K, P K , Σ K ) associated with each domain K∈ M h Given a function u defined in Ω, regular enough, its interpolant u h ∈ X h is uniquely defined such as [6]: • The restriction u h j K , that is, the form of uh in the geometric element K, belongs to the space P K ; • The restriction u h j K 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 M h 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.
The finite element space X h 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∈ M h , that is It is also possible to define the base functions p h,j , 1 ≤ i ≤ Dh, of the space Xh from the base functions of the spaces P K , ∀ 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 ≤ D h , are well defined, it can be associated a function u h , called X h -interpolant, defined by

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].
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.

Base functions of spaces S i
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 n i = {i}, is associated the function The finite dimensional space generated by all sn i 's is denoted by S 0 .
With edge e ij = {i, j}, is associated the vector field where N F,m n 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   edge {m, n} in common. Consequently, field s eij is zero in all the elements non adjacent to edge e ij .
The vector field space generated by all s e is denoted by S 1 .
With facet f = f ijk = {i, j, k} = {q 1 , q 2 , q 3 } or f = f ijkl = {i, j, k, l} = {q 1 , q 2 , q 3 , q 4 }, is associated the vector field where #Nf is the number of nodes of facet f, a f = 2 if #N f = 3, a f = 1 if #N f = 4, and the list of q i 's is made circular by setting q 0 q#N f and q#N f + 1 q 1 . Field s f is zero in all the elements non adjacent to facet f.
Vector fields s fijk(l) 's generate the space S 2 .
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: sn i is equal to 1 at node ni, and to 0 at other nodes; the circulation of s eij is equal to 1 along edge eij, and to 0 along other edges; the flux of s fijk(l) is equal to 1 across facet s fijk(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 6 ¼ j. These properties show up various kinds of functionals and involve that functions s n , s e , s f , s v 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.

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] gradP F,mn ¼ grad 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 N F,m n . 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).
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 s eij . The fact that its circulation along edge eij is equal to 1 needs some calculation to be proved.
The vector field p q c gradP F,q c q cþ1 Â gradP F,q c q cÀ1 (44) which appears in expression (35) of s f , 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.
The flux of this cross product, and in consequence the one of (44), is then equal to 0 across these facets. The term pq c 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 s f across facet f is then the only one to differ from zero (Figure 11).

Degrees of freedom
The expression of a field in the base of a space S i -S 0 or S 3 for a scalar field, S 1 or S 2 for a vector field-gives scalar coefficients, called degrees of freedom. Fields ϕ ∈ S 0 , h ∈ S 1 , j ∈ S 2 and ρ ∈ S 3 can be expressed as [4,6,9] Φ ¼ X j ¼ The degrees of freedom ϕ n , h e , j f 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, Figure 11. Geometric interpretation of the facet function s f (36) [6].
as mentioned in the definition of finite elements, are thus respectively pointwise evaluations, line, surface and volume integrals.

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 s e and for the normal component of s f . 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 S 1 can represent vector fields like magnetic field h or electric field e whose tangential components are continuous across interfaces between materials, and those of S 2 can represent fields like induction field b or current density field j whose normal components are continuous across interfaces between these materials.

Spaces S i 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, X e ∈ E i n, e ð Þs e ¼ grads n (49) 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.
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.

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] x ¼ where x∈K, u∈Kr; this element is isoparametric.

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.

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).
Edge-node incidence matrix ð54Þ Node i∈N Nodal base function p i (u, v, w) = s i (u, v, w)  Facet-edge incidence matrix Volume-facet incidence matrix ð56Þ Figure 14.
Geometric entities defined on a tetrahedron of type I [6].

Reference 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).

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

Node i∈N
Nodal base function p i (u, v, w) = s i (u, v, w)  Table 6.
Notation of the facets of the hexahedron of type I.
Facet-edge incidence matrix Volume-facet incidence matrix

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).

Node i∈N
Nodal base function p i (u, v, w) = s i (u, v, w)

Finite Element Methods and Their Applications
Facet-edge incidence matrix Volume-facet incidence matrix

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].  Table 9.
Notation of the facets of the prism of type I.  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:26 MS=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.
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 19.
Modeling of TEAM problem 7: Coil and conducting plate [10]. 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).

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.  [4].