The finite element method (FEM) is one of the most powerful numerical methods to solve practical engineering problems due to its strong capability of handling complex and multiscale structures  and inhomogeneous media . In the computational electromagnetics (CEM), the FEM is widely used to solve various electromagnetic problems, like the electrostatic currents , magneticstatic problems , electromagnetic scattering , integrated circuits modeling , electrical-thermal co-simulations .
However, investigations found that the traditional FEM is quite sensitive to the quality of background meshes . More specifically, equilateral triangles and regular tetrahedrons are much preferred to guarantee the accuracy. However, when complex or multiscale structures are discretized into many nonoverlapping elements, a large number of irregular elements may exist. Then, the accuracy of the FEM can be severely degenerated and even totally unacceptable. However, generation of high good quality triangles and tetrahedrons is quite challenging when multiscale and complex structures are involved . Therefore, the results cannot be always trusted when the mesh quality cannot be guaranteed. In addition, when low order basis functions are used, only a low level of accuracy of results can be obtained. That implies that a large number of unknowns are required if we need accurate results .
In this paper, we address these two issues in the CEM by introducing the edge-based smoothed FEM (ES-FEM) to accurately solve the electrostatic problems in the two dimensional cylindrical system and three dimensional cartesian system. The ES-FEM was firstly introduced to solve non-electromagnetic problems by properly combining the traditional FEM and mesh-free methods, termed as smoothed point interpolation methods(S-PIMs). Various versions of smoothed FEM (SFEM) are proposed based on how to construct smoothing domains, like node-based FEM (NS-FEM) , face-based FEM (FS-FEM) , edge-based FEM (ES-FEM) . They all show better performance in terms of accuracy and numerical stability for thermal analysis , acoustic analysis  and computational mechanics . In most of applications, the ES-FEM shows best performance among them. In , the FS-PIM is introduced for the electromagnetic analysis. It is found that the FS-PIM shows improved accuracy. We introduced the ES-FEM to model the high-speed interconnects  and the electrostatic lens . In this paper, we introduced the ES-FEM for electromagnetic analysis and comprehensively investigate its numerical properties, which has significant performance improvement compared with the FEM.
This paper is organized as follows. In Section II, detailed formulations for the ES-FEM are presented. In Section III, its accuracy and numerical properties are comprehensively investigated through several numerical experiments. At last, we draw some conclusions in Section IV.
Ii-a Problem Configurations
There are various axis-symmetrical systems, like electrostatic lens 
, which can be simplified as two dimensional problems in cylindrical system. The boundary-value problems can be defined by the following partial differential equation (PDE) as
subject to the boundary condition
where denotes the potential in the computational domain, is the excitation, and are constants for the third kind of boundary conditions which is the first kind of boundary condition with and , and the second kind of boundary condition with and .
For those fully three dimensional structures, the PDE for the three dimensional cartesian system can be defined as
subject to the boundary condition
where , , and are constant material parameters, denotes the potential, is the excitation, and are constants.
Ii-B Formulations of the FEM
In the FEM, we formulate (1) and (3) into a set of linear equations with the expansion of through linear basis functions in each element through the Ritz method . In matrix form, those linear equations can be rewritten into
where is assembled from its elemental entity in cylindrical system given by
In three dimensional cartesian system, is expressed as
is a linear shape function in the FEM.
Ii-C Creation of Smoothing Domain
Assuming that the solution domain has been divided into nonoverlapping elements with nodes and edges. For two dimensional domain, by connecting the centroid of each element with its neighbor nodes, as shown in Fig. 1(a), the smoothing domain associated with the th edge is created. Therefore, the whole computational domain is divided into nonoverlapping smoothing domains. Each smoothing domain may consist of 3 or 4 segments and nodes depending on the location of the edge .
For three dimensional domain, one of the sub-domain of the smoothing domain for th edge located in the th tetrahedron is created by connecting the centroid of the tetrahedron element and two centroids of the triangle surfaces, which are shared by edge and vertexes of edge , as shown in Fig. 1(b). By assembling all sub-domains of th edge, the th smoothing domain is created.
In traditional FEM, each node has only interactions with those belonging to the same element. However, in the ES-FEM, each node interacts with those contributing to the same smoothing domain. That means the interactions are performed on the smoothing domain rather than triangle or tetrahedron elements, which makes the ES-FEM more accurate and more stable to irregular meshes, especially, for multiscale structures without increasing the overall count of unknowns compared with the traditional FEM.
Ii-D Detailed Formulations of the ES-FEM
The difference between the ES-FEM and the FEM lies in how we construct the stiff matrix which contains the gradient of the shape function . We first introduce the generalized gradient smoothing technique to replace the gradient component by its smoothed counterpart . In the ES-FEM formulation, the smoothed stiff matrix can be expressed as
The generalized gradient smoothing technique  is applied in each smoothing domain by averaging the gradient shape function
where denotes the smoothing function defined as
where is area or volume of the th smoothing domain in two or three dimensional domains, respectively. Through the divergence theorem, the gradient smoothed can be expressed as
where is the boundary of smoothing domain ,
is the unit normal vector pointing to outside of the smoothing domain on.
where in two dimensional cylindrical domain and in three dimensional domain, is the segment count of , and are abscissas and weight for Gaussian numerical integration, respectively, is the th component of unit outward vector, is the length of the th segment. Therefore, we can denote the stiff matrix as
where is the number of nodes which contribute to the th smoothing domain and is the smoothed gradient matrix which can expressed as
in two and three dimensional systems, respectively, is the number of elements that share the edge , is the area or volume of the th triangle or tetrahedron element.
By applying the boundary conditions and excitation exactly the same as those in the traditional FEM, we can solve the electromagnetic problems through the ES-FEM.
Iii Numerical Results and Discussion
Iii-a A Cubical Box
To analytically investigate the numerical properties of the ES-FEM, we first consider a cubical box as shown in Fig. 2.
The top surface of cubical box is applied to a boundary condition
and other surfaces are set to 0 . Therefore, the potential inside the box is given by
Four background meshes with different mesh sizes are considered in our numerical analysis. In addition, to investigate the numerical stability of the ES-FEM, four irregular meshes, which are generated from a small random offset enforced at each node except boundary nodes in the regular meshes, are used. The averaged edge length of all elements is denoted as . The relative error of of all the nodes in the computational domain is defined as
where is the numerical result and is the reference potential obtained by (17).
Fig. 3 quantitively illustrates the percentage of different ratio of elements in the four regular and irregular meshes used in the FEM and the ES-FEM. The ratio is defined as the maximum to minimum edge length in each element, which denotes mesh quality. Obviously, the ratio should be as near as one for good quality meshes. As shown in Fig. 3, the ratio of only around 20% of elements is larger than 2. However, in irregular meshes, the ratio of over 40% of elements is larger than 2, which implies that much more distorted elements exist in the irregular meshes.
Fig. 4 shows the relative error of the traditional FEM and the ES-FEM. It is easy to find that the ES-FEM outperforms the traditional FEM for both regular and irregular meshes. For regular meshes, the ES-FEM is more accurate than the FEM. When it comes to irregular meshes, it becomes more obvious. As shown in Fig.4 , the accuracy of the FEM is severely degenerated up to almost one order. However, the ES-FEM shows great numerical stability and its accuracy is not affected by the irregular meshes as shown in Fig. 4, which shows great numerical stability.
Iii-B A Multiscale Electronic Lens
A multiscale electronic lens is considered which consists of an emitter, a suppression electrode, an absorbing electrode and an anode, as is shown in Fig. 5(a). Detailed geometry configurations are demonstrated in Fig. 5(b). The size of the whole structure is 12 millimeters, but the size of the emitter tip is only 0.0003 millimeters. The ratio of the size of the electronic lens to the size of the emitter tip reaches , which is a typical multiscale structure.
As we can see, an electronic lens is an axis-symmetric structure which can be simplified into a two dimensional model. The potential of the emitter, the absorbing electrode and the anodes are set to 0 V, 7,000 V and 70,000 V, respectively. Both the cylindrical ES-FEM and fully three dimensional ES-FEM are applied in this case. Since the potential at the axis is quite important for the design, we extract the potential at the central axis to verify the accuracy of FEM and ES-FEM.
Fig. 6 shows that the ratio of elements in the four meshes used in our simulations. As shown in Fig. 6, the ratio of almost 40% elements is larger than 2. Further investigations find most irregular elements exist near the tip of the emitter. Since the lens is a multiscale structure, the quality of meshes cannot be improved no matter how we refine the meshes. As we mentioned in the introduction section, it is challenging to generate good quality meshes for multiscale structures. Since the FEM is sensitive to mesh quality, the accuracy cannot be guaranteed for multiscale structures.
The whole solution domain is first divided into 290,953 nonoverlapping tetrahedrons with 64,765 nodes and 387,700 edges. The overall count of unknowns in the FEM and the ES-FEM are 64,765. The potential along the central axis obtained from the ES-FEM with the linear shape function, the traditional FEM with the linear shape function. The reference solution is obtained from the Comsol with 278,488 unknowns and the curved second order shape function.
Potential and relative errors along the axis are shown in Fig. 7. Both the 2D and 3D ES-FEM can obtain much more accurate results, especially in the mesh distortion region, near the tip of the emitter, compared with the FEM. The relative errors of the traditional FEM reach 2, which is almost 20 times of that of the ES-FEM results. The ES-FEM can nearly immune to the distortion mesh due to the smoothing gradient technique. Furthermore, numerical results show that the ES-FEM in two dimensional system performs better than its three dimensional counterpart due to the geometry symmetry which reduces the error caused by mesh discretization in multiscale structures. The relative errors of the 3D ES-FEM with 4 different meshes is shown in Fig. 8. It is easy to find that the ES-FEM is much more accurate than the traditional FEM.
Iii-C A Micromotor
We consider a electrostatic micromotor, which is a fully three dimensional structure . Due to its symmetrical property, only a quarter of the structure is considered, which has 3 poles at stator and 2 poles at the rotor. The radius of stator is 100 m with 20 m in height. The inner and outer radius of the rotor are 20 m and 50 m with 18 m in height, respectively. The potential of the central stator tooth is set to 100 V and other teeth are set to 0 V. Other boundaries are set as Neumann boundary conditions.
The whole solution domain is first discretized into 45,519 nonoverlapping tetrahedrons with 8,839 nodes and 56,815 edges. The potential distribution is computed at m plane as shown in Fig. 10. The potential is also evaluated at m in the m plane. The reference solution is the numerical solution computed by the FEM with a fine mesh with 509,451 unknowns. As shown in Fig. 10, the ES-FEM with much less unknowns can obtain the exactly the same pattern as the reference. Furthermore, as shown in Fig. 11, the potential obtained from the ES-FEM agree well with the reference solution. It implies that the ES-FEM can solve complex electromagnetic structures.
In this paper, we introduce the ES-FEM to accurately solve two dimensional cylindrical and three dimensional electromagnetic problems. By applying the generalized gradient smoothing technique to smoothing the gradient component in the FEM and constructing the smoothing domains based on all edges in the background meshes, a highly accurate and numerically stable ES-FEM is constructed. As numerical results show, the ES-FEM can indeed obtain much more accurate results and is stable to irregular meshes, especially when a large number of irregular meshes may exist in the multiscale structures. Therefore, it shows great potential to solve practical electromagnetic problems.
The authors would like to thank financial supports from Beijing Natural Science Foundation through Grant 4194082, National Natural Science Foundation of China through Grant 61801010, and Fundamental Research Funds for the Central Universities.
-  J. M. Jin, The finite element method in electromagnetics. John Wiley&Sons, 2014.
-  S. H. Lee, K. Mao, and J. M. Jin, “A Complete Finite-Element Analysis of Multilayer Anisotropic Transmission Lines From DC to Terahertz Frequencies,” IEEE Trans. Adv. Packag., vol. 31, pp. 326-338, Jun. 2008.
-  T. Wang, R. F. Harrington, and J. R. Mautz, “Quasistatic analysis of a microstrip via through a hole in a ground plane,” IEEE Trans. Microw. Theory Tech., vol. 36, no. 6, pp. 1008-1013, 1988.
-  P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers. 2nd ed., Cambridge: Cambridge University Press, 1990.
-  A. F. Peterson and S. P. Castillo, “A frequency-domain differential equation formulation for electromagnetic scattering from inhomogeneous cylinders,” IEEE Trans. Antennas Propagat., vol. 37, pp. 601–607, Jun. 1989.
-  S. H. Lee, K. Mao, and J. M. Jin, “A complete finite-element analysis of multilayer anisotropic transmission lines from DC to terahertz frequencies,” IEEE Trans. Adv. Packag., vol. 31, pp. 326–338, Jun. 2008.
-  T. Lu, J. M. Jin, “Electrical-Thermal Co-Simulation for Analysis of High-Power RF/Microwave Components,” IEEE Trans. Electromagn. Compat., vol. 59, no. 1, pp. 93-102, Sep. 2016.
-  G. Liu, Meshfree methods: moving beyond the finite element method. 2nd edn. CRC Press, 2009.
-  G. Liu, H. Nguyen-Xuan, and T. Nguyen-Thoi, “A theoretical study on the smoothed FEM (S-FEM) models: Properties, accuracy and convergence rates,” Int. J. Numer. Meth. Engng, vol. 84, pp. 1222-1256, Dec. 2010.
-  G. Liu, N. Trung, Smoothed finite element methods. CRC press, 2016.
-  G. Liu, T. Nguyen-Thoi, and H. Nguyen-Xuan, “A node-based smoothed finite element method (NS-FEM) for upper bound solutions to solid mechanics problems,” Comput. Meth. Appl. Mech. Eng., vol. 87, no. 2, pp. 14-26, Jan. 2009.
-  T. Nguyen-Thoi, G. Liu, and K. Lam, “A facebased smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements,” Int. J. Numer. Methods Eng., vol. 78, no. 3, pp. 324–353, Apr. 2009.
-  Z. He, G. Liu, and Z. Zhong, “An edge-based smoothed finite element method (ES-FEM) for analyzing three-dimensional acoustic problems,” Comput. Meth. Appl. Mech. Eng., vol. 199, pp. 20–33, Dec. 2009.
-  E. Li, G. Liu, “Simulation of hyperthermia treatment using the edge-based smoothed finite-element method,” Numer. Heat Tranf. A-Appl., vol. 57(11), pp. 822–847, Jun. 2010.
-  Z. He, A. Cheng, and Z. Zhong, “Dispersion error reduction for acoustic problems using the edge-based smoothed finite element method (ESFEM)” Int. J. Numer. Methods Eng., vol. 86, no. 11, pp. 1322–1338, Jun. 2011.
-  T. Nguyen-Thoi, G. Liu, and H. Nguyen-Xuan, “Additional properties of the node-based smoothed finite element method (NS-FEM) for solid mechanics problems,” Int. J. Numer. Methods Eng., vol. 6, no. 4, pp. 633-666, 2009.
-  Lima, Naísses Z., and Renato C. Mesquita, “Face-based gradient smoothing point interpolation method applied to 3-D electromagnetics,” IEEE Trans. Magn., vol. 50, no. 2, pp. 537-540, Feb. 2014.
-  Y. Zhang, S. Yang, and D. Su, “An Edge-based Smoothed FEM for Accurate High-Speed Interconnect Modeling,” in Conf. Numerical Electromagnetic and Multiphysics Modeling and Optimization, Boston, U.S.A., May. 2019.
-  Y. Zhang, P. Wang, W. Li and S. Yang, “An Edge-Based Smoothed FEM for Multiscale Electrostatic Lens Modeling,” in Conf. ICEAA-IEEE APWC, Granada, Spain, Sep. 2019.
-  F. Guimaraes, R. Saldanha, R. Mesquita, D. Lowther, and J. Ramirez, “A meshless method for electromagnetic field computation based on the multiquadric technique”, IEEE Trans. Magn., vol. 43, no. 4, pp. 1281–1284, May. 2007.