Journal cover Journal topic
Advances in Radio Science An open-access journal of the U.R.S.I. Landesausschuss in der Bundesrepublik Deutschland e.V.
Journal topic
https://doi.org/10.5194/ars-17-59-2019
https://doi.org/10.5194/ars-17-59-2019

19 Sep 2019

19 Sep 2019

# Towards a Spectral Method of Moments using Computer Aided Design

Towards a Spectral Method of Moments using Computer Aided Design
Stefan Kurz, Sebastian Schöps, and Felix Wolf Stefan Kurz et al.
• Technische Universität Darmstadt, Institute for Accelerator Science and Electromagnetic Fields and Centre for Computational Engineering, Darmstadt, Germany

Abstract

We present first numerical examples of how the framework of isogeometric boundary element methods, in the context of electromagnetism also known as method of moments, can be used to achieve higher accuracies by elevation of the degree of basis functions. Our numerical examples demonstrate the computation of the electric field in the exterior domain.

1 Introduction

Spectral methods, cf. , are classes of methods for solving differential equations, closely related to classical element based methods. They share the same idea: The approximation of the solution through a series of basis functions. While classical methods choose to refine a mesh, and with this to further localise the support of each individual basis function, spectral methods employ global basis functions to approximate the solution of the problem. Since such a global basis is not always readily available, quite often any finite and boundary element method which relies solely on p-refinement, i.e., the increase of the degree of the local without mesh (h-) refinement, is referred to as a spectral element method.

In engineering applications, spectral element methods are rarely considered; the reason simply being that to fully enjoy their convergence properties, meshes with curved elements of increasing orders must be generated. This poses challenges to mesh generation and pre-processing. In contrast, classical h-refinement based mesh generators are well understood. However, with the introduction of isogeometric analysis by the problem of efficient mesh generation can be avoided by the use of exact geometry mappings, allowing computations directly on CAD-generated objects.

In this document, we present numerical experiments which showcase how the isogeometric framework can be used to obtain an implementation of a spectral boundary element method. We demonstrate the implementation by the solution of electromagnetic scattering problems through p-refinement. For this, we employ a solution strategy via the electric field integral equation , which is also referred to as method of moments (MOM), and an isogeometric boundary element framework . We essentially employ p-refinement to a patchwise polynomial basis, in our case based on Bernstein polynomials. Other variants of spectral MOM have already been studied, eg. by or .

The organisation of the paper is straight forward. We first introduce basic notions of the electric field integral equation, and our p-refinement based discretization scheme, built on top of the framework of isogeometric analysis. Afterwards, we comment on the matrix assembly, followed by a discussion of our numerical examples.

Figure 1Visualisations of the geometries. The maximum diameter of the boat is smaller than 4.5, the sphere is of radius 1.

Figure 2Visualisations of the real part of the unknown surface current.

2 The Electric Field Integral Equation

We consider the scattering of an electromagnetic wave under the assumption of constant material coefficients μ and ϵ in Ωc, i.e., the surroundings of a scatterer Ω, PEC boundary condition on $\mathrm{\Gamma }:=\partial \mathrm{\Omega }$ and the Silver-Müller radiation condition. Prescribing an incident wave g we arrive at

with wave number $\mathit{\kappa }:=\mathit{\omega }\sqrt{\mathit{ϵ}\mathit{\mu }}.$ Herein, n denotes the outward unit normal of Γ. Under these conditions, there exists a surface current j such that the scattered field can be represented by the electric field integral equation (EFIE), given by

with

$\begin{array}{ll}\stackrel{\mathrm{̃}}{{V}}\left(\mathbit{j}\right)\left(\mathbit{x}\right)& :=\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right)\mathbit{j}\left(\mathbit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{y}+\frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}{\mathbf{grad}}_{x}\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right)\\ \text{(3)}& & \cdot {\mathrm{div}}_{\mathrm{\Gamma }}\left(\mathbit{j}\left(\mathbit{y}\right)\right)\mathrm{d}{\mathrm{\Gamma }}_{y},\end{array}$

for all $\mathbit{x}\notin \mathrm{\Gamma }.$ Herein, ${G}_{\mathit{\kappa }}\left(\cdot ,\cdot \right)$ denotes the Green's function given by ${G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right):=\frac{{e}^{i\mathit{\kappa }|\mathbit{x}-\mathbit{y}|}}{\mathrm{4}\mathit{\pi }|\mathbit{x}-\mathbit{y}|},$ cf. for details.

In a continuous setting, this can be recast as the variational problem of finding an unknown surface current j in the trace of the space H(curl,Ω), such that

$\begin{array}{ll}& \underset{\mathrm{\Gamma }}{\int }\left[\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right)\mathbit{j}\left(\mathbit{x}\right)\cdot \mathit{\xi }\left(\mathbit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{y}+\frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right)\right\\ & \phantom{\rule{1em}{0ex}}\left({\mathrm{div}}_{\mathrm{\Gamma }}\circ \phantom{\rule{0.125em}{0ex}}\mathbit{j}\right)\left(\mathbit{x}\right)\left({\mathrm{div}}_{\mathrm{\Gamma }}\circ \phantom{\rule{0.125em}{0ex}}\mathbit{\xi }\right)\left(\mathbit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{y}]\mathrm{d}{\mathrm{\Gamma }}_{x}=\\ \text{(4)}& & -\underset{\mathrm{\Gamma }}{\int }\left({\mathbit{e}}_{i}×\mathbit{n}\right)\cdot \mathbit{\xi }\mathrm{d}\mathrm{\Gamma }\end{array}$

holds for all test functions ξ in the same space. The corresponding trace space is often denoted by ${H}_{×}^{-\mathrm{1}/\mathrm{2}}\left({\mathrm{div}}_{\mathrm{\Gamma }},\mathrm{\Gamma }\right)$ and requires a divergence-conforming discretisation.

Figure 3Numerical examples on the unit sphere. Wave number κ=1. The error refers to the maximum error obtained by comparing the values of the analytical solution with the field values through evaluation of Eq. (2) on a selection of 10 000 points on a sphere of radius 2 around the origin. Computed on a desktop PC with 16 GB RAM and an Intel i7 8700k processor.

Figure 4Numerical examples for the ship geometry as presented in Fig. 1a. Highest p generates a discretisation with 5600 complex-valued degrees of freedom. Wave number κ=1. The error refers to the maximum error obtained by comparing the values of the analytical solution with the field values through evaluation of Eq. (2) on a selection of 10 000 points on a sphere of radius 6 around the origin. Computed on a desktop PC with 16 GB RAM and an Intel i7 8700k processor.

3 Discretisation

To solve the electric wave equation via a boundary element approach, the unknown is reduced to a vector field on Γ, often discretised by divergence-conforming elements. Implementations of such and related numerical schemes are, among others, given by , , or . We utilise an approach based on the framework of isogeometric analysis, dealing with the special case where the B-spline basis reduces to the set of Bernstein polynomials, for order p≥0 given by

${b}_{\mathit{\alpha },p}=\left(\begin{array}{c}n\\ \mathit{\alpha }\end{array}\right){x}^{\mathit{\alpha }}\left(\mathrm{1}-x{\right)}^{p-\mathit{\alpha }},\phantom{\rule{1em}{0ex}}\mathrm{for}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{0}\le \mathit{\alpha }\le p.$

To apply p-refinement, one needs sufficiently smooth (patchwise) geometry mappings, since otherwise they limit the order of ansatz functions which can be applied effectively. Thus, we choose geometry mappings given by mappings from the unit square to parts of the geometry Γj parametrized as tensor product mappings of rational Bernstein polynomials. Such representations can be easily extracted from NURBS (non uniform rational B-Spline) parametrisation through Bézier extraction .

Following , a suitable discretisation of the trace spaces on so-called multipatch domains $\mathrm{\Gamma }={\bigcup }_{\mathrm{0}\le j, i.e., elements in the context of spectral element methods, is provided in . By tensor product construction one constructs discrete spaces that are conforming w.r.t. the differential operators. As an example, if Sp denotes the Bernstein polynomials of order p≥1 on (0,1), one finds that

${S}^{p}\otimes {S}^{p}\stackrel{\mathbf{curl}}{⟶}\left(\begin{array}{c}{S}^{p}\otimes {S}^{p-\mathrm{1}}\\ {S}^{p-\mathrm{1}}\otimes {S}^{p}\end{array}\right)\stackrel{\mathrm{div}}{⟶}{S}^{p-\mathrm{1}}\otimes {S}^{p-\mathrm{1}}$

holds. This way, using the NURBS geometry mappings induced by the CAD geometries to seamlessly map between (0,1)2 and patches Γj⊆Γ, one can define a conforming discretisation of the entire de Rham complex

$\begin{array}{}\text{(5)}& {H}^{\mathrm{1}}\left(\mathrm{\Omega }\right)\stackrel{\mathbf{grad}}{⟶}H\left(\mathbf{curl},\mathrm{\Omega }\right)\stackrel{\mathbf{curl}}{⟶}H\left(\mathrm{div},\mathrm{\Omega }\right)\stackrel{\mathrm{div}}{⟶}{L}^{\mathrm{2}}\left(\mathrm{\Omega }\right),\end{array}$

as well as its trace spaces

$\begin{array}{}\text{(6)}& {H}^{\mathrm{1}/\mathrm{2}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathbf{curl}}_{\mathrm{\Gamma }}}{⟶}{H}_{×}^{-\mathrm{1}/\mathrm{2}}\left({\mathrm{div}}_{\mathrm{\Gamma }},\mathrm{\Gamma }\right)\stackrel{{\mathrm{div}}_{\mathrm{\Gamma }}}{⟶}{H}^{-\mathrm{1}/\mathrm{2}}\left(\mathrm{\Gamma }\right),\end{array}$

which is required for the analysis of boundary element methods. In the case of the divergence-conforming space, one must require additional normal continuity of the vector field across patch interfaces. This can be achieved through identification of the corresponding degrees of freedom (DOFs) with one another. We denote the discretisation on the boundary defined by

$\begin{array}{}\text{(7)}& {S}_{p}^{\mathrm{0}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathbf{curl}}_{\mathrm{\Gamma }}}{⟶}{S}_{p}^{\mathrm{1}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathrm{div}}_{\mathrm{\Gamma }}}{⟶}{S}_{p}^{\mathrm{2}}\left(\mathrm{\Gamma }\right).\end{array}$

For these spaces, it is possible to show existence, uniqueness, and quasi-optimality of the solution . The discrete problem to Eq. (4) is that of finding a discrete surface current ${\mathbit{j}}_{h}\in {S}_{p}^{\mathrm{1}}\left(\mathrm{\Gamma }\right)$ such that

$\begin{array}{ll}& \underset{\mathrm{\Gamma }}{\int }\left[\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right){\mathbit{j}}_{h}\left(\mathbit{x}\right)\cdot {\mathbit{\xi }}_{h}\left(\mathbit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathbit{y}}+\frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}\left(\mathbit{x},\mathbit{y}\right)\left({\mathbf{div}}_{\mathrm{\Gamma }}\circ \phantom{\rule{0.125em}{0ex}}{\mathbit{j}}_{h}\right)\right\\ \text{(8)}& & \phantom{\rule{1em}{0ex}}\left(\mathbit{x}\right)\left({\mathbf{div}}_{\mathrm{\Gamma }}\circ \phantom{\rule{0.125em}{0ex}}{\mathbit{\xi }}_{h}\right)\left(\mathbit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathbit{y}}]\mathrm{d}{\mathrm{\Gamma }}_{\mathbit{x}}=-\underset{\mathrm{\Gamma }}{\int }\left({\mathbit{e}}_{i}×\mathbit{n}\right)\cdot {\mathbit{\xi }}_{h}\mathrm{d}\mathrm{\Gamma }\end{array}$

holds for all test functions ${\mathbit{\xi }}_{h}\in {S}_{p}^{\mathrm{1}}\left(\mathrm{\Gamma }\right).$ The coefficients to represent ${\mathbit{j}}_{h}=\sum {c}_{j}{\mathbit{\xi }}_{h,j}$, where the ξh,j denote the basis functions of ${S}_{p}^{\mathrm{1}}\left(\mathrm{\Gamma }\right),$ can be obtained in form of the vector c given by the solution to the linear system Ac=r, where A and r can be assembled in direct analogy to Eq. (8).

4 Matrix Assembly

We apply a modified version of the superspace approach as used by e.g. to construct the dense system matrix A resulting from Eq. (4) via the representation

$\mathbf{A}={\mathbf{P}}^{\mathit{⊺}}{\mathbf{A}}^{\mathit{†}}\mathbf{P}={\mathbf{P}}^{\mathit{⊺}}\left(\sum _{i,j=\mathrm{0}}^{\mathit{#}\text{elements}}{\mathbf{A}}_{j,i}^{\mathit{†}}\right)\mathbf{P},$

where A is the dense system matrix w.r.t an elementwise polynomial, globally discontinuous basis, and P is the superspace matrix assembling the divergence-conforming basis. Herein, each matrix ${\mathbf{A}}_{j,i}^{\mathit{†}}$ is sparse, including only the interaction of the local basis on element i with that of element j. This way, the small dense system A is assembled without the need to store A as a whole. In the special case of only p-refinement, the matrix P incorporates only linear combinations necessary to reduce the degree in one tensor product direction for the construction of ${S}_{p}^{\mathrm{1}}\left(\mathrm{\Gamma }\right)$, and to achieve normal continuity across patch interfaces. For the solution of the arising linear system, we utilise a partially pivoted LU decomposition of the Eigen linear algebra library .

5 Numerical Examples

To showcase the possibility of obtaining an increased accuracy through (mainly) p-refinement, we present a simple numerical example in analogy to the ones for h-refinement presented by .

We excite our model setup by a Hertz-Dipole as defined by Jackson (1998, p. 411) as

$\begin{array}{ll}{\mathbit{E}}_{\mathrm{DP}}\left(\mathbit{x}\right)& :={e}^{i\mathit{\kappa }r}\left(\frac{{\mathit{\kappa }}^{\mathrm{2}}}{r}\left(\mathbit{n}×{\mathbit{p}}_{\mathrm{0}}\right)×\mathbit{n}+\left(\frac{\mathrm{1}}{{r}^{\mathrm{3}}}-\frac{i\mathit{\kappa }}{{r}^{\mathrm{2}}}\right)\\ \text{(9)}& & \left(\mathrm{3}\mathbit{n}\left(\mathbit{n}\cdot {\mathbit{p}}_{\mathrm{0}}\right)-{\mathbit{p}}_{\mathrm{0}}\right)\right),\end{array}$

with $r=‖\mathbit{x}-{\mathbit{x}}_{\mathrm{0}}‖$ and $\mathbit{n}=\left(\mathbit{x}-{\mathbit{x}}_{\mathrm{0}}\right)/r$. The dipole's singularity is placed inside Ω. Away from its singularity, specifically within the exterior domain, the dipole fulfils Eq. (1). Thus, by existence and uniqueness of the solution of the exterior problem Eq. (4), cf. , we know that Eq. (8) will approximate the surface current required to represent the field ${{\mathbit{E}}_{\mathrm{DP}}}_{|{\mathrm{\Omega }}^{\mathrm{c}}},$ cf. . This construction using a prescribed solution is known as a method of manufactured solutions, cf. Oberkampf and Roy (2010, Chap. 6.3).

Summarised, we solve for j in Eq. (8) with an excitation given by Eq. (9). We then evaluate Eq. (2) numerically with the approximated surface current jh whose induced eh(x) yields the numerical solution. By existence and uniqueness eh is an approximation to EDP in the exterior. We compute the error ${max}_{\mathbit{x}\in V}\parallel {\mathbit{E}}_{\mathrm{DP}}\left(\mathbit{x}\right)-{\mathbit{e}}_{h}\left(\mathbit{x}\right){\parallel }_{{C}^{\mathrm{3}}}$ at points x in a set V containing points placed on a sphere enclosing the geometry, see the captions of Figs. 3 or 4 for details.

As a first example, we choose the example of a unit sphere given by 6 patches and define a dipole with ${\mathbit{x}}_{\mathrm{0}}={\mathbit{p}}_{\mathrm{0}}=\left(\mathrm{0},\mathrm{0.1},\mathrm{0.1}\right)$ as an excitation. Although a boundary element framework has been utilised to solve the problem, no compression was applied to the systems due to their small system size, cf. Fig. 3c.

Thus the condition of the system matters little compared to the case in which compression must be applied.

On the sphere, a stable exponential rate of convergence w.r.t. p is observed. The application of p-refinement reduces the time required for matrix assembly as well as the overall system size, cf. Fig. 3, significantly.

We present another example. As geometry, we choose the boat depicted in Fig. 1 consisting of 28 Bézier patches. Placing the dipole inside with ${\mathbit{x}}_{\mathrm{0}}=\left(\mathrm{1},\mathrm{0},\mathrm{0}\right)$ and ${\mathbit{p}}_{\mathrm{0}}=\left(\mathrm{0},\mathrm{0},\mathrm{0.1}\right)$, i.e., under the “bridge”, we evaluate the electric field around the geometry. For this non-smooth and non-convex geometry, the rate of convergence, as seen in Fig. 4, is not as pronounced as before. Still, one can see an exponential convergence behaviour paired with excellent times to solution.

6 Conclusions

The adaptation of isogeometric to spectral element methods is straight forward. Through the use of Bézier extraction, one extracts piecewise smooth parametrisations of geometries. Using known constructions from isogeometric analysis, one can define a global, divergence-conforming basis that consists of tensor product Bernstein polynomials, with supports consisting of one patch each, or in the case of functions identified to achieve normal continuity, multiple patches. Application of p-refinement to this basis yields a spectral method of moments, which achieves high accuracies without the need for system compression.

These results are exceptionally promising since the reduced system size makes the application of direct solvers possible, thus circumventing the need for preconditioners, which are still a challenging topic for boundary element methods for electromagnetic problems.

Code availability
Code availability.

The code basis, as well as the geometries for these computations are available at http://www.bembel.eu/ (last access: 21 May 2019). A fork of the repository ready to recreate the computations can be found at https://github.com/flx-wlf/bembel/tree/ars2019aa/ (last access: 21 May 2019). See also (https://doi.org/10.5281/zenodo.2671596).

Author contributions
Author contributions.

All authors have jointly carried out research and worked together on the manuscript. The numerical tests have been conducted by the last author. All authors read and approved the final manuscript.

Competing interests
Competing interests.

Stefan Kurz is also affiliated as Chief Expert with Robert Bosch GmbH.

Special issue statement
Special issue statement.

Acknowledgements
Acknowledgements.

The work of Felix Wolf is supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at TU Darmstadt.

Financial support
Financial support.

This research has been supported by the DFG (grant no. SCHO1562/3-1) and the DFG (grant no. KU1553/4-1).

Review statement
Review statement.

This paper was edited by Thomas Eibert and reviewed by three anonymous referees.

References

Benoit, C., Royer E., and Poussigue, G.: The spectral moments method, J. Phys. Condens. Matter., 4, 3125–3152, 1992. a

Borden, J., Scott, M. A., Evans, J. A., and Hughes T. J. R.: Isogeometric finite element data structures based on Bézier extraction of NURBS, Int. J. Numer. Meth. Eng., 87, 15–47, 2011. a

Buffa, A. and Hiptmair, R.: Galerkin boundary element methods for electromagnetic scattering, Lect. Notes Comp. Sci., 31, 83–124, 2003. a, b, c

Buffa, A., Rivas, J., Sangalli, G., and Vázquez, R.: Isogeometric discrete differential forms in three dimensions, SIAM J. Numer. Anal., 49, 818–844, 2011. a

Buffa, A., Dölz, J., Kurz, S., Schöps, S., Vázquez, R., and Wolf, F.: Multipatch approximation of the de Rham sequence and its traces in isogeometric analysis, submitted, preprint available: arXiv:1806.01062, 2018. a, b

Di Ruscio, D., Burghignoli, P., Baccarelli, P., Comite, D., and Galli, A.: Spectral Method of Moments for Planar Structures With Azimuthal Symmetry, IEEE T. Antenn. Propag., 62, 2317–2322, https://doi.org/10.1109/TAP.2014.2302831, 2014. a

Dölz, J., Kurz, S., Schöps, S., and Wolf, F.: A Numerical Comparison of an Isogeometric and a Classical Higher-Order Approach to the Electric Field Integral Equation, submitted, preprint available: arXiv:1807.03628, 2018a. a, b

Dölz, J., Kurz, S., Schöps, S., and Wolf, F.: Isogeometric Boundary Elements in Electromagnetism: Rigorous Analysis, Fast Methods, and Examples, submitted, preprint available: arXiv:1807.03097, 2018b. a, b, c

Dölz, J., Harbrecht, H., Kurz, S., Multerer, M., Schöps, S., and Wolf, F.: Bembel v0.9, https://doi.org/10.5281/zenodo.2671596, 2019. a

Guennebaud, G. and Jacob, B.: Eigen v3, available at: http://eigen.tuxfamily.org (last access: 21 May 2019), 2010. a

Hughes, T. J. R., Cottrell, J. A., and Bazilevs Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Meth. Appl. Mech. Eng., 194, 4135–4195, 2005. a

Jackson, J. D.: Classical Electrodynamics, Wiley and Sons, New York, 3rd edition, 1998. a

Oberkampf, W. L. and Roy, C. J.: Verification and Validation in Scientific Computing, Cambridge University Press, Caimbridge, 2010. a

Trefethen, L. N.: Spectral Methods in MATLAB, SIAM, Philadelphia, 2000. a

Tzoulis, A. and Eibert, T. F.: A Hybrid FEBI-MLFMM-UTD Method for Numerical Solutions of Electromagnetic Problems Including Arbitrarily Shaped and Electrically Large Objects, IEEE T. Antenn. Propag., 53, 3358–3366, 2005. a

Hiptmair, R. and Kielhorn, L.: BETL – A generic boundary element template library, Seminar for Applied Mathematics, ETH Zürich, Rep. no. 36, 2012. a

Weggler, L.: High Order Boundary Element Methods, Dissertation, Universität des Saarlandes, Saarbrücken, 2011. a