Comparison Between NEM and FEM in 2-D Magnetostatics Using an Error Estimator

This communication deals with a comparison between two methods of discretization: the well known finite element method and the natural element method that is a meshless method. An error estimator, based on the nonverification of the constitutive law, is used. This estimation has been applied to two examples: a device with permanent magnets and a variable reluctance machine.


I. INTRODUCTION
N OWADAYS, the finite element method (FEM) is well established and allows solving many problems in applied physics. However, difficulties can appear with the FEM when the system includes moving parts. Actually, a remeshing step can be required to avoid a significant distortion of elements or its overlapping. Recently, novel methods have been proposed to circumvent, or at least alleviate, the main difficulties related to the mesh constraints. The most usual mesh constraints are related to the need of efficient remeshing procedures, which in the 3-D case constitutes a difficulty nowadays under active investigation, as well as to the field projection from the old mesh to the new mesh that induces an inevitable numerical diffusion whose unfavorable impact on the solution after numerous resmeshing steps is widely accepted. It was in this context that the so-called "meshless methods" were introduced one decade ago. Their applications in several domains experienced a spectacular progression in the last decade, after proving its ability for avoiding, or at least for alleviating significantly, the remeshing requirements. The natural element method [1] (from now on referred as NEM) belongs to this new and vast family of meshless methods and its interest compared with the vast majority of meshless techniques concern its interpolation property that allows enforcing essential boundary conditions in an easy way, as in the finite element method. We recall that the imposition of essential boundary conditions is the main issue in the other meshless methods. On the other hand, it has been proved that the solution accuracy does not depend on the quality of the subjacent triangulation ("meshless character"). Thus, we could affirm, crudely speaking, that the NEM is a kind of FEM in which the triangles could be distorted without a significant impact on the solution accuracy.
The NEM, which was firstly applied in the framework of solid and fluid mechanics [2], [3], has already been successfully tested in electromagnetic systems for accounting, in a simple manner, the existence of moving parts [4]. However, in the electromagnetic framework the NEM and FEM solution accuracy comparison has never been addressed. This is a question needing a clear response and consequently justifying or disappointing the use of meshless approaches in this field. In this communication, we compare both methods NEM and FEM by using an error estimator based on the non verification of the constitutive law on two examples [5], [6]. The first example concerns two permanent magnets in opposition and the second example a Variable Reluctance Machine.

II. PROBLEM DESCRIPTION AND ERROR ESTIMATOR
The magnetostatic problem, defined in the domain D, can be written as and (1) and (2) is the magnetic flux density, is the magnetic field, and is the current density that is assumed to be known in the domain . is bounded by the surface where represents the outward normal unity vector. Moreover, and verify the constitutive relationship given by: with the coercitive magnetic field and the permeability. In the sequel, the problem is assumed to be linear. Let's call "admissible" fields and that verify (1) and (2) respectively. It is obvious that an admissible field couple verifying the constitutive relationship (3) is the exact solution denoted . We consider now an admissible couple and we define a constitutive law error 'e': (4) with It can be proven [5] that this error e is equal to the distance between an average field and the exact solution : with (5) Therefore, the scalar e is a numerical error estimator for the admissible solution. To obtain the admissible field couple, several methods have been already proposed. The more accurate method consists in solving both potential formulations. The vector potential formulation ( and on ) yields an admissible induction field . The scalar potential formulation ( and on with such that and on ) yields an admissible magnetic field . Then, using (4), an estimation of the numerical error can be calculated. To evaluate the global solution quality, let us also define a global relative error estimator (6) III. DISCRETIZATION In 2-D, both potential formulations derived from the model defined in (1), (2) and (3), can be solved by discretizing the associated Galerkin's weak form. For this purpose we need to define the approximation to be considered for both the trial and the test functions. This approximation is usually built from the values that the unknown fields have at nodes distributed on the domain D. A shape function ( is the position vector) is associated to each node . Thus, both potentials can be approximated from a linear combination of such shape functions (7) The only difference between FEM and NEM lies on the expression of the shape functions that in both methods are defined in a different manner. The ones related to the FEM are well known and therefore they do not need additional explanations. On the contrary, as the ones related to the NEM being less known, a brief introduction seems to be suitable. The NEM shape functions have been introduced by Sibson and are calculated using the Voronoï diagram [1]. We are illustrating the concepts related to the construction of such shape functions by considering the simple distribution of five nodes depicted in Fig. 1. Some basic constructions, as the Delaunay triangulation and the associated Voronoi diagram, that serve to define the NEM shape functions in 2-D will be illustrated.

A. Voronoï Diagram
We define the Voronoï diagram associated to a set of nodes as the geometric subdivision of the  plane into regions , called first order Voronoï cells, each one related to a distinct node , such that (8) with the coordinates of the node and the distance between node and point . Thus, is the region of the plane that contains the points closest to node than to any other node in . In Fig. 1 the nodes are numbered from "1" to "5" and the Voronoï cell related to node 5 is depicted. The construction of that Voronoï cell is quite simple. We start identifying the natural neighbors [1] of node 5. Now, we consider the segments joining node 5 to each neighbor node, and a straight line, normal to each one of these segments, traced at the central point of those segments as depicted in Fig. 1. These straight lines defined the Voronoï cell. Fig. 2 illustrates the construction of the Voronoï cells related to the other nodes. For additional details the interested reader can refers to [1] and the references therein.

B. Sibson Shape Functions
Once the Voronoï diagram is constructed, we can calculate the value of the NEM shape function , related to the node at the position . For this purpose we consider, in the scenario depicted in Fig. 3, the insertion of a point and the computation of the shape function related to node 5 at that position, that is, . We consider two possible placements of such . After introduction of point , the Voronoï diagram is modified because the consideration of the cell related to the point just introduced that in our example consists of the polygon "abcd". The cell just inserted, the one related to point , intersects the former cells related to the original nodes. In Fig. 3 we emphasize in grey color the intersection between this cell and the one initially related to node 5. Now, the shape function related to node 5 at position is computed from the ratio between the intersected area (the grey area in Fig. 3) and the one related to new cell (polygon "abcd"), that is It is easy to prove [1] that the support of the shape function related to a node is given by the union of the circles passing through this node and any other two neighbor nodes of such that no other node falls within the circle. Fig. 4(a) depicts those circles and the corresponding support for node 5. By computing , by using the procedure just described at different positions , the shape function can be represented. Fig. 4(b) depicts the shape function related to node 5. One of the reasons of the robustness of the Sibson interpolation against the distortion of the background mesh defining the Delaunay triangles and the associated Voronoï diagram, is that even if a triangle becomes very distorted, the support of the shape functions is not restricted to this triangle (as it is the case in the finite element method). Nowadays there is not a mathematical proof of the higher robustness of this kind of interpolation; however, there are numerous examples showing the low sensibility of NEM to the mesh distortion [7]. Several simulations in solid and fluid mechanics [2], [3] have been successfully performed using the same cloud of nodes throughout the whole simulation.

C. Nodal Sibson Shape Function Properties
The NEM shape functions are except at the positions of nodes where they are and on the boundary of the Delaunay circles where they are . The shape functions have a compact support delimited approximately by the neighbour nodes, leading to a sparse system. The band width of the resulting linear systems is a little bit higher than the bandwith associated to the FEM. That is precisely one of the reason of the where is the position of node and the Kronecker's delta. The so-called delta property (the first one in (10)) enables to enforce essential boundary conditions as in FEM. The other two properties in (10) prove the linear completeness, that is, the Sibson interpolation can represent exactly linear functions. Moreover, the partition of unity, second relation in (10) allows functional enrichments in the context of the partition of unity paradigm.
The computing time related to the construction of the Sibson shape functions is higher than the one associated with the definition of such functions in the finite element framework. However, nowadays, there are fast algorithms coming from the computational geometry community that are robust, fast even in 3-D, and some of them freely distributed. If we take the CPU time as comparison criterion and consider simulations involving massively remeshning operations (that define the natural context of application of meshless methods), finite elements simulations are nowadays faster than natural element method, because the existence of efficient and robust remeshing procedures in 2-D. If one consider convergence and accuracy criteria, the convergence rate for linear FEM and NEM is the same (first order in the energy norm) but as we proved in some of our former works in mechanics [8] (and we will prove also in magnetostatics later in this paper) the NEM accuracy is higher than the one related to the finite elements. Moreover, in models requiring field transfer after each remeshing process the superiority of the NEM is more relevant as argued previously.

IV. APPLICATIONS
First, two magnets face to face are considered (Fig. 5). Ten meshes with an increasing number of nodes were used. Fig. 5 gives the Voronoï diagram associated to the clouds consisting of 107, 800 and 2200 nodes. For the FEM simulation, linear interpolation on the Delaunay triangles was used.
For each mesh, both potential formulations have been solved using both the NEM and the FEM. Fig. 6 depicts the evolution of energy with the number of nodes. We can notice that the energies calculated with both formulations bound the exact energy. However, the energies computed from NEM seem to be more accurate that the ones coming from FEM calculations.  By representing the evolution of the numerical error computed from (6) versus the number of nodes in a log-log scale (see Fig. 7) for both the finite element and the natural element solutions, we notice that both evolutions are linear and that the convergence rate, slope of such evolutions, is the same and approximately equal to one. We can also notice the higher accuracy of the NEM solution.
A more complex model in now considered. It concerns the variable reluctance machine (VRM) depicted in Fig. 8. The modelling of VRM is a very trick problem because of regions with sharp angles where the concentration of the magnetic flux density is very dramatic. Some simple error indicators are available [9] allowing the nodal refinement. The relative permeability of the rotor and stator yokes is equal to 1000. Different nodal distributions have been tested (the coarsest one containing 6000 nodes and the finest one more than one million) and both potential formulations solved by using the NEM and the FEM. The errors and have been computed from (5). The error ratio evolves from a value close to 2 related to the coarsest nodal distribution, to 2.5 for the finest one.

V. CONCLUSION
In this paper, the accuracy of meshless natural element method has been compared to the FEM one. For the comparison purposes we proposed the use of an error estimator