Anomalous Magnetization Spikes in the Locally Corrected Nyström Discretization of Static Volume Integral Equation using Tetrahedral Cells

John C. Young1, Robert J. Adams1, and Stephen D. Gedney2

1Department of Electrical & Computer Engineering, University of Kentucky, Lexington, Kentucky 40506, USA
john.c.young@uky.edu, robert.adams@uky.edu

2Department of Electrical Engineering, University of Colorado – Denver, Denver, Colorado 80204, USA

Submitted On: July 31, 2023; Accepted On: September 26, 2023

ABSTRACT

A locally corrected Nyström (LCN) discretization of a magnetostatic volume integral equation is presented. Anomalous magnetization spikes can occur when the underlying mesh uses tetrahedral cells regardless of discretization order. The mechanism for the anomalous magnetization spikes is discussed, and mitigation of the spikes through use of an LCN-to-Moment Method conversion is investigated. Results are presented validating that the LCN-to-Moment Method suppresses the anomalous spikes.

Index Terms: integral equation methods, locally corrected Nyström method, moment method.

I. INTRODUCTION

The locally corrected Nyström (LCN) method [13] is one of the primary methods for discretizing integral equations. Advantages of the LCN method over the Moment Method (MoM) include not having to explicitly enforce continuity of physical quantities across mesh cell boundaries, less strict mesh conformality requirements, more efficient system matrix fill, and ease of implementing higher-order codes. Divergence-conforming formulations require that normal components of quantities such as fields and currents be continuous (unless physically discontinuous) across mesh cell boundaries. Although the Nyström method does not explicitly enforce this normal continuity, the underlying Nyström degrees-of-freedom should permit the proper space and continuity properties required by the formulation to be achieved.

In this paper, a magnetostatic volume integral equation is presented that is discretized by the locally corrected Nyström method [4]. Results are provided in which anomalous (non-physical) magnetization [5] are observed for complex geometries when tetrahedral mesh elements are used in the discretization. The mechanism for the anomalous magnetization spikes is discussed by noting that the typical Nyström representation does not span the same mixed-order divergence-conforming space that commonly used Moment Method bases span. Further, use of an LCN-to-MoM conversion [7, 8] of the discretized LCN system is observed to suppress the anomalous magnetization spikes since the LCN-to-MoM conversion provides a representation with the appropriate degrees-of-freedom to model a mixed-order divergence-conforming space as well as the appropriate normalcontinuity.

II. LCN FORMULATION FOR TETRAHEDRA

Consider the magnetostatic volume integral equation (VIE)

χ-1M(r)=Hexcitation(r)+Hm(r),rV, (1)

defined over a material volume V where M is the magnetization, χ is the magnetic susceptibility tensor, and Hm is the demagnetizing field

Hm(r)=V[14π|r-r|M(r)]dv. (2)

The VIE in (1) is discretized using the locally corrected Nyström (LCN) method. For simple geometries, good results and convergence have been obtained for both hexahedral and tetrahedral cells. For sphere and spherical shell geometries with isotropic, homogeneous magnetic susceptibility, high-order convergence is obtained for higher basis orders when higher-order mesh representations are used. For some complex geometries, however, the magnetization is observed to exhibit anomalous spikes for tetrahedral meshes. Similar spikes in the magnetization are not observed for hexahedral meshes of the same geometry when using a mixed-order LCN formulation [10].

images

Figure 1: Primary tetrahedral cell Tc adjacent to four secondary tetrahedral cells Tk. Shared face Sk is shared by Tc and Tk.

For basis order p = 0, the phenomenon seems to be at least partially due to the inability of the underlying function space of the Nyström representation of the magnetization to sufficiently model the continuity of the normal component of magnetization across tetrahedral cell boundaries. More generally, for order p0 the standard Nyström representation on tetrahedral cells (and similarly on triangular cells) does not span the same space as typical divergence-conforming MoM representations (e.g,, [9]). The ability to maintain continuity across all cell boundaries is discussed for the basis order p = 0 and p = 1 representations. For orders p > 1, analogous results hold. Further, the number of degrees-of-freedom of the standard Nyström representation and mixed-order divergence-conforming MoM representations are not the same, indicating different spaces are being spanned.

Consider the set of tetrahedral cells depicted in Fig. 1. A primary cell Tc is adjacent to four secondary cells Tk for k = 0,1,2,3. The boundary face Sk is shared by cells Tc and Tk. Local coordinates are defined such that uix is the coordinate along the ith unitary axis for i=1,2,3. A dependent local coordinate u0x=1-u1x-u2x-u3x is also defined. The local coordinates are such that ukx=0 on face Sk in cell Tx. For i=1,2,3, uix(u1x,u2x,u3x) is the ith unitary vector and uix(u1x,u2x,u3x)is the ith reciprocal unitary vector in cell x. Furthermore, u0x=-(u1x+u2x+u3x). The cell vertices are ordered such that the outward normal to boundary Sk is -ukc in cell Tc and is -ukk in cell Tk, and, so, ukc=-ukk on Sk.

Let the Nyström degrees-of-freedom be cast onto a set of bases within tetrahedron x with degrees-of-freedom αix,j for the jth degree-of-freedom associated with the ith unitary direction. The pth-order polynomial-complete basis representation on tetrahedron x is

Mxp(u1x,u2x,u3x)=1gxi=13fixp(u1x,u2x,u3x)uix, (3)

where gx(u1x,u2x,u3x) is the cell Jacobian and fixp is a polynomial of degree p. Note that on the boundary face Sk, gc=gk at each point on the face. Hence, in the following development, the cell Jacobians will cancel for constraints on Sk. For both the p = 0 and p = 1 discussion, the degrees-of-freedom in cell Tk are considered fixed, and the degrees-of-freedom in cell Tc will be constrained (if possible) to achieve continuity of normal magnetization across all faces.

For p = 0, there are three Nyström degrees-of-freedom and three bases, and

fix0(u1x,u2x,u3x)=αix,00. (4)

Enforcing continuity at each face Sk for k=0,1,2,3 gives the constraints

[ukc1gci=13αic,00uic]ukc=0
=[ukc1gki=13αik,00uik]ukk=0. (5)

Since the set of unitary and reciprocal unitary vectors are orthonormal, the constraints in (2) reduce to

αkc,00=-αkk,00,k=1,2,3, (6)

and

k=13αkc,00=-k=13αk0,00,k=0. (7)

Note there are only three degrees-of-freedom in cell Tc but four constraints that must be satisfied. The three degrees of freedom within the cell permit the normal component of the magnetization to be matched continuously at three of the tetrahedral cell faces. However, there are insufficient degrees-of-freedom within the cell to match the normal component of magnetization at all four of the faces of the cell as required in a divergence-conforming formulation.

For p = 1, there are twelve Nyström degrees-of-freedom and twelve bases, and [3]

fix1(u1x,u2x,u3x)=αix,01+αix,11u1x+αix,21u2x+αix,31u3x. (8)

Enforcing continuity at each face Sk for k=0,1,2,3 gives the constraints

[ukcgci=13(αic,01+αic,11u1c+αic,21u2c+αic,31u3c)uic]ukc=0=[ukcgki=13(αik,01+αik,11u1k+αik,21u2k+αik,31u3k)uik]ukk=0. (9)

The constraints in (9) then reduce to

(αic,01+αic,11u1c+αic,21u2c+αic,31u3c)|ukc=0=-(αik,01+αik,11u1k+αik,21u2k+αik,31u3k)|ukk=0, (10)

for k=1,2,3. For each k=1,2,3, three degrees-of-freedom on Tc are specified leaving 3 remaining degrees-of-freedom to match continuity across S0. For S0, the constraint reduces to

i=13(αic,01+αic,11u1c+αic,21u2c+αic,31u3c)u0cuic=i=13(αi0,01+αi0,11u10+αi0,21u20+αi0,31u30)u0cui0. (11)

Given the normal

u0c=-(u1c+u2c+u3c)=-u00, (12)

at S0 and u1x+u2x+u3x=1 on the shared face, there are three constraints and maintaining normal continuity is possible. Note that [5] mistakenly indicates normal continuity is not possible for p1 using the polynomial-complete representation.

In the p = 0 case, the Nyström degrees-of-freedom are insufficient by one DoF to be able to enforce normal continuity at all four faces simultaneously. The p = 1 case allows a linear tangential, linear normal representation instead of the typical linear tangential, constant normal representation. However, the linear tangential, linear normal representation includes degrees-of-freedom associated with the null space of the divergence operator, which is undesirable.

In general, the order p divergence-conforming interpolatory vector basis set for tetrahedral cells in [9] give the number of boundary face bases and internal cell bases as

2(p+1)(p+2), (13)

and

p(p+1)(p+2)/2, (14)

respectively. Table 1 lists the number of degrees-of-freedom for the Nyström representation and a divergence-conforming interpolatory [9] representation on a tetrahedral cell for basis orders p = 0 through p = 2. It is noted that there are fewer Nyström bases than divergence-conforming bases. What is missing from the Nyström function space are mixed-order p+1 terms that ensure that the divergence of the basis function space is complete to order p. Without these terms, the divergence of the Nyström basis space is only of order p-1. In addition, for the p = 0 basis, the divergence is zero.

Table 1: Total degrees-of-freedom (DoF) in Nystrom and an interpolatory divergence-conforming representations versus basis order p. (Number of boundary face and internal cell DoF in parenthesis)

Basis Order p Nyström DoF Interpolatory Divergence Conforming DoF
0 3 4 (4,0)
1 12 15 (12, 3)
2 30 36 (24, 12)

For p > 0, while the Nyström basis is sufficient to represent the magnetization to order p and provide normal continuity, it has insufficient DoF to represent the charge to order p. Therefore, the Nyström basis has an insufficient number of DoF to represent the mixed-order divergence-conforming space to order p. Furthermore, it over specifies the p-1 function space. Hence, the standard Nyström representation on tetrahedra (and triangles) cannot be truly divergence-conforming although the method seems to provide very good solutions across a wide range of problems.

For higher-order LCN representations, a variety of quadrature rules with differing numbers of points and properties exist thus complicating the discussion. For Nyström discretizations with p > 0, it is desirable to choose a quadrature rule of degree q = 2p, but, for p > 1, rules of degree q=2p usually lead to non-square local correction matrices. Non-square local correction matrices may compromise the stability of the solution, so choosing quadrature rules of degree q(p+1) such that the number of points equals the number k=1p+1p(p+1)(p+2)/6 of Nyström bases for each vector component at order p is recommended. The quadrature rules need not be symmetric even though symmetric rules are preferred when available.

The magnetization spike phenomenon is not observed for hexahedral meshes. For example, for the p = 0 Nyström representation on hexahedral cells, a mixed-order representation [10] has six degrees-of-freedom which are sufficient to match continuity of the normal magnetization at each of the six faces of the cell. If the Nyström degrees-of-freedom are cast onto a set of bases within the cell, the basis representation is

Mx0=1gxi=13(αix,00+αix,10uix)uix, (15)

Along each unitary direction there are two degrees-of-freedom allowing for a constant plus linear representation of the magnetization so that the magnetization at one face can vary sufficiently to its opposite face (unlike the p = 0 tetrahedral representation). Furthermore, the hexahedral Nyström degrees-of-freedom span the same space that a typical divergence-conforming Moment Method basis set spans. A polynomial-complete representation on hexahedral cells, however, does not span the proper divergence-conforming space and may suffer from other spurious effects [10].

III. LCN-TO-MOM FORMULATION

Nyström methods are desirable since the system matrix fill avoids the costly double integrations that arise in Method of Moment (MoM) discretizations. An LCN-to-MoM conversion [7, 8] allows an LCN system to be easily converted to a MoM system without sacrificing too many of the LCN method’s advantages. Further, for the magnetostatic VIE formulation presented, it is not necessary to strictly enforce the continuity of MoM bases across tetrahedral cell faces. Hence, independent half-MoM bases can be assigned to each face in each tetrahedral cell, which further simplifies the LCN-to-MoM conversion. For the magnetostatic formulation, independent half-MoM bases must be assigned to shared faces between two cells that have different susceptibilities.

In the LCN-to-MoM conversion on tetrahedral cells, an order p MoM representation requires the LCN system to be filled at order (p+1) to accommodate the polynomial space requirements of divergence-conforming bases [9]. In view of Table 1, the (p+1) Nyström representation has too many degrees-of-freedom for the analogous divergence-conforming representation. Hence, the LCN-to-MoM conversion matrices can be viewed as the appropriate constraints to remove the extraneous degrees-of-freedom in a Nyström representation and achieve a divergence-conforming representation.

images

Figure 2: Average relative field error in Nyström solution of (1) for a magnetic spherical shell vs. maximum mesh edge length for various basis orders p and mesh orders o. The relative mesh discretization error for the spherical shell surface area is also plotted.

images

Figure 3: Subsection of circular-cylindrical shell with extending circular frustrum shell overlaid with tetrahedral mesh.

images

Figure 4: Magnitude of magnetization in circular-cylindrical shell section with circular-cylindrical frustrum shell extension for (a) tetrahedral LCN simulation (b) tetrahedral LCN-to-MoM simulation, and (c) hexahedral LCN simulation.

IV. DISCUSSION

First, a convergence analysis was performed for a locally-corrected Nyström (LCN) discretization of (1) for a magnetic spherical shell. The inner radius is 0.9 m, the outer radius 1.0 m, and the relative permeability is 50. The shell was meshed with a sequence of three tetrahedral meshes with 1514, 3320, and 6426 cells, respectively for mesh orders o = 1 (linear) and o = 2 (quadratic). The convergence analysis was performed for LCN basis orders p = 0, 1, and 2. The magnetic field was computed at various points outside the shell, and the average relative error was calculated using the analytic solution [4] as a reference. The results of the convergence analysis are plotted in Fig. 2. Also plotted are the relative mesh discretization error of the total surface area of the shell. The average relative error is observed to be limited by the mesh discretization error. However, the convergence rate increases with basis order until limited by the mesh discretization error. Hence, when the geometry is fairly smooth and uncomplicated, the LCN method usually produces good results.

As a second example, a subsection of a circular-cylindrical shell with a hollow circular-cylindrical frustrum extending outward, depicted in Fig. 3, was analyzed. The cylinder has a height of 6 m, an outer radius of 5 m, a wall thickness of 2 cm and is aligned on the z axis and centered at the origin. The frustrum has a height of 3 m, a wall thickness of 2 cm, and extends out from the cylinder to a final outer radius of 5.5 m. Only the subsection of the structure for x>3 m is retained. The circular-cylindrical shell has a relative permeability of 150, and the frustrum extension has a relative permeability of 100. The excitation Hexcitation=-z^A/m, and the geometry was meshed with a 28728 cell, linear tetrahedralmesh.

The LCN simulation was performed at p = 0, and the LCN-to-MoM conversion used the Schaubert-Wilton-Glisson (SWG) bases [6] with four degrees-of-freedom per tetrahedral cell. Independent half-MoM bases were applied to each shared face. Plots of the magnetization for the LCN solution and LCN-to-MoM solution are provided in Fig. 4 (a) and Fig. 4 (b), respectively. In Fig. 4 (c) is shown the magnetization for a p = 0 LCN solution using a hexahedral mesh. The LCN-to-MoM solution and hexahedral LCN solutions are observed to be regular while the tetrahedral LCN solution exhibits anomalous spikes in magnetization across the mesh. The maximum magnetization magnitude for the tetrahedral LCN-MoM solution is approximately 156 A/m while the tetrahedral LCN solution spikes to almost 1275 A/m. Furthermore, some of the spikes in LCN magnetization occur at smooth parts of the mesh and not near edges or other complex features.

images

Figure 5: Scattered magnetic field Hz along a line centered at x = 4.9 m, y = 0 m for the cylindrical shell section with frustrum extension example. Included are the fields due to the tetrahedral LCN solution, the tetrahedral LCN-MoM solution, and the hexahedral LCN solution.

In Fig. 5 is plotted the z-component of the scattered magnetic field vs z along a line centered at (x,y)=(4.9,0)m for the tetrahedral LCN, tetrahedral LCN-to-MoM, and hexahedral LCN solutions. The line passes close to some of the magnetization spikes observed in Fig. 4 (a) for the tetrahedral LCN solution. While the scattered field of the LCN-to-MoM and hexahedral LCN are visually identical, the field of the tetrahedral LCN solution is seen to be corrupted by the spuriousmagnetization.

V. CONCLUSION

The locally-corrected Nyström (LCN) method was observed to potentially produce anomalous spikes in magnetization when used to solve magnetostatic volume integral equations with tetrahedral mesh cells to model geometry. The failure of the LCN method for tetrahedra was discussed in terms of the inability of the Nyström representation to appropriately model a typical divergence-conforming space. Further, it was also observed that the p = 0 LCN discretization for tetrahedra does not allow continuity of normal magnetization between mesh cells to be appropriately modeled. For p > 0, the LCN basis does allow continuity, but fails to represent the charge to order p, and does not fully represent the divergence-conforming basis to order p.

An LCN-to-MoM discretization, however, does appropriately model a mixed-order divergence-conforming space since the MoM bases are constructed to be divergence-conforming. Further, proper continuity of the normal magnetization across cell boundaries is modeled at all orders even though it is not required to explicitly enforce the MoM bases to be continuous across cell boundaries. Hence, many of the advantages of an LCN method are maintained in an LCN-to-MoMdiscretization.

Convergence results for a spherical magnetic shell were presented that showed the LCN method usually performs well for simple structures. Higher-order convergence was achieved until stagnated by mesh discretization error. However, for a more complicated model, the LCN method produced anomalous spikes while the LCN-to-MoM did not.

In conclusion, it is remarkable that the typical LCN representations used with tetrahedral meshes perform as well as they do as the representation do not span the appropriate space. Still, over a wide range of problems, no issues are observed, and good convergence characteristics are achieved. It is extremely difficult to predict a priori for a specific geometry whether the LCN method will produce spurious results. For large complex geometries with millions of cells, the probability of anomalous spikes occurring greatly increases. Although mesh quality can affect both the LCN and LCN-to-MoM solution, mesh quality seems to be more important to obtain robust LCN results when using tetrahedral cells. Refining the mesh may help suppress the spikes in some cases, but the increase in system size, as well as not knowing a priori whether the refinement is sufficient, is prohibitive for very large problems. Further, while the spikes seem to be local in nature and do not seem to corrupt the whole solution, fields in the vicinity of the spikes may have significant errors. Hence, the use of the LCN-to-MoM is advisable when working with tetrahedral meshes.

ACKNOWLEDGMENT

This work was supported in part by Office of Naval Research Grants N00014-16-1-3066 and N00014-23-1-2138.

REFERENCES

[1] L. F. Canino, J. J. Ottusch, M. A. Stalzer, J. L. Visher, and S. M. Wandzura, “Numerical solution of the Helmholtz equation in 2D and 3D using a high-order Nyström discretization,” (in English), Journal of Computational Physics, vol. 146, no. 2, pp. 627-663, Nov. 1, 1998.

[2] S. D. Gedney and J. C. Young, “The locally corrected Nyström method for electromagnetics,” in Computational Electromagnetics: Recent Advances and Engineering Applications, R. Mittra, Ed. New York: Springer, 2014, pp. 149-198.

[3] M. S. Tong, Z.-G. Qian, and W. C. Chew, “Nyström method solution of volume integral equations for electromagnetic scattering by 3D penetrable objects,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 5, pp. 1645-1652, May 2010.

[4] J. C. Young and S. D. Gedney, “A Locally Corrected Nyström formulation for the magnetostatic volume integral equation,” IEEE Transactions on Magnetics, vol. 47, no. 9, pp. 2163-2170, Sep. 2011.

[5] J. C. Young, R. J. Adams, and S. D. Gedney, “Anomalous current spikes in the kocally corrected Nyström discretization of volume integral equations,” 2023 International Applied Computational Electromagnetics Society Symposium, Denver, CO, March 26-March 30, 2023.

[6] D. Schaubert, D. Wilton, and A. Glisson, “A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE Transactions on Antennas and Propagation, vol. 32, no. 1, pp. 77-85, 1984.

[7] R. A. Pfeiffer, J. C. Young, R. J. Adams, and S. D. Gedney, “Locally corrected Nyström to moment method conversion for volume integral equations,” IEEE Transactions on Magnetics, vol. 55, no. 4, pp. 1-7, 2019.

[8] M. Shafieipour, I. Jeffrey, J. Aronsson, and V. I. Okhmatovski, “On the equivalence of RWG method of moments and the locally corrected Nyström method for solving the electric field integral equation,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 2, pp. 772-782, Feb. 2014.

[9] R. D. Graglia and A. F. Peterson, Higher-Order Techniques in Computational Electromagnetics, Sci-Tech Publishing, 2016.

[10] S. D. Gedney, A. Zhu, and C.-C. Lu, “Study of mixed-order basis functions for the locally corrected Nyström method,” IEEE Transactions on Antennas and Propagation, vol. 52, no. 11, pp. 2996-3004, Nov. 2004.

BIOGRAPHIES

images

John C. Young received the B.E.E. degree in electrical engineering from Auburn University in 1997, the M.S. degree in electrical engineering from Clemson University in 2000, and the Ph.D. degree in electrical engineering also from Clemson University in 2002.

From January 2003 to April 2003, he served as a post-doctoral researcher at Clemson University, and from 2003 to 2005, he served as a post-doctoral researcher at Tokyo Institute of Technology, Tokyo, Japan. From 2005 to 2008, he worked at Japan Radio Co. Since 2008, he has been with the Department of Electrical and Computer Engineering at the University of Kentucky, Lexington, KY where he is currently an associate professor.

Dr. Young’s research interests include integral equation methods, finite element methods, electromagnetic theory, waveguides, array antennas, and magnetic signature modeling of hysteretic materials. He is a member of IEEE, the Applied Electromagnetics Society (ACES), and URSI Commission B. He currently serves as an Associate Editor for the IEEE Transactions on Antennas and Propagation and on the Education Committee of the Antennas and Propagation Society. He also served (2020-2023) on the Board of Directors of ACES where he is currently Secretary.

images

Robert J. Adams received the B.S. degree from Michigan Technological University, Houghton, MI, USA, in 1993, and the M.S. and Ph.D. degrees in electrical engineering from Virginia Polytechnic Institute and State University (Virginia Tech), Blacksburg, VA, USA, in 1995 and 1998, respectively.

From 1999 to 2000, he was a Research Assistant Professor with Virginia Tech. Dr. Adams joined the University of Kentucky in 2001, where he is currently a Professor with the Department of Electrical and Computer Engineering.

He has made novel contributions to mesh and frequency stable integral equation formulations of electromagnetic problems, constraint-based methods for high-order MOM discretizations, spectral splitting methods for implementing shadowing effects in integral equations at high frequencies, and sparse direct solution methods for low-to-moderate frequency electromagnetics applications. Dr. Adams is a senior member of the IEEE.

images

Stephen D. Gedney received the B.Eng.-Honors degree from McGill University, Montreal, P.Q., in 1985, and the M.S. and Ph.D. degrees in Electrical Engineering from the University of Illinois, Urbana-Champaign, IL, in 1987 and 1991, respectively.

He is currently the Don and Karen White Professor of the Department of Electrical Engineering at the University of Colorado Denver (CUD). Previously he was a Professor of Electrical Engineering at the University of Kentucky from 1991 – 2014. He worked for the U.S. Army Corps of Engineers, Champaign, IL (’85-’87). He was a visiting Professor at the Jet Propulsion Laboratory, (92’, 93’), HRL laboratories (’96-’97) and Alpha Omega Electromagnetics (’04-’05). He received the Tau Beta Pi Outstanding Teacher Award in 1995 and 2013. From 2002 – 2014, he was the Reese Terry Professor of Electrical and Computer Engineering at the University of Kentucky. He was titled as a Distinguished Professor of the University of Colorado in 2022. He is a past Associate Editor of the IEEE Transactions on Antennas and Propagation (1997 – 2004), a member of the IEEE Antennas and Propagation Society ADCOM (2000 – 2003), and served as the chair of the IEEE Antennas and Propagation Society Membership Committee (1995 – 2002). He is a Fellow of the IEEE and member of Tau Beta Pi.

ABSTRACT

I. INTRODUCTION

II. LCN FORMULATION FOR TETRAHEDRA

images

III. LCN-TO-MOM FORMULATION

images

images

images

IV. DISCUSSION

images

V. CONCLUSION

ACKNOWLEDGMENT

REFERENCES

BIOGRAPHIES