MANIFOLD HARMONICS FOR INTEGRAL EQUATION BASED ISOGEOMETRIC ANALYSIS OF ACOUSTIC AND ELECTROMAGNETIC PHENOMENA AND SHAPE OPTIMIZATION By A. M. A. ALSNAYYAN A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering โ€“ Doctor of Philosophy 2023 ABSTRACT Loop subdivision has become a standard tool in computer graphics for generating (second order) smooth surfaces of arbitrary topology. A few of its most attractive features include that it is inherently multiresolution, scalable, allows for local adaptive refinement, compression, and can readily handle smooth deformations of the surface. These features make Loop subdivision an attractive candidate for use in an isogeometric analysis (IGA) framework, in which geometry and physical quantities (e.g., surface currents, pressures, charges, etc.) are expressed in the same basis. The goal of this thesis is to develop a collection of efficient and accurate methods that inherit these features for both analysis and design in acoustics and electromagnetics problems. This thesis addresses some of the fundamental problems that arise in IGA of acoustic and electromagnetic phenomena. First, we develop representation of scalar and vector quantities on both simply connected and multiply connected Loop subdivision surfaces. This leverages the high order nature of Loop subdivision to construct a complete Helmholtz decomposition that avoids issues that plague standard methods. We use this complete Helmholtz decomposition to discretize well-conditioned integral equations for scattering analysis. Next, we address the challenge of computational cost associated with evaluating the many additional inner products that arise in these integral equation formulations; this bottleneck is overcome by using a wideband multilevel fast multipole algorithm (MLFMA) and mixed-potential field representations to accelerate computations in the iterative solver. Then, we introduce the notion of manifold harmonics as a means to (i) compress geometry and integral equation matrices, (ii) enrich physical quantities hierarchically for acoustic and electromagnetic scattering, and (iii) utilize local enrichment to better handle representation of physical quantities on structures with challenging features. Finally, we utilize the manifold harmonics basis on a Loop subdivision surface for the shape reconstruction problem; specifically, they allow for multi-resolution representation (and editing) of complex surfaces and functions defined thereupon. Taken together, these contributions address several open problems and yield an integrated design-analysis approach that is highly accurate, flexible, and fast with applicability to a range of engineering problems. Copyright by A. M. A. ALSNAYYAN 2023 This thesis is dedicated to Mom and Dad. Thank you for always believing in me. v ACKNOWLEDGMENTS I would like to express my sincere gratitude to my advisor Dr. Shanker Balasubramaniam for his unwavering guidance and encouragement, both in and outside the academic world. Working with such a talented advisor and researcher has been inspiring. Without him, truly none of this work would be possible. I would also like to express my deepest appreciation to my committee which included Dr. Alex Diaz, Dr. Leo Kempel, Dr. Jeffery Nanzer, and Dr. Edward Rothwell. I would also like to show appreciation to all the friends and colleagues that I have met at Michigan State University, both in and out the Electromagnetics Research Group. Namely, I would like to thank Dr. Scott Oโ€™Connor, Dr. Zane Crawford, Dr. Steve Hughey, Dr. Jie Li, Dr. Connor Glosser, Dr. Stavros Vakalis, Omkar Ramachandran, Luke Baumann, Doga Dikbayir, Jacob Hawkins, Anton Schlegel, Daniel Chen, Sean DePalma, and Jason Merlo. Finally, I would like to thank my family. My siblings Gharam, Ghalih, and Essam for supporting me. My parents, Mahmoud and Aysheh, to whom words fail to do justice in expressing my ever lasting appreciation for all they have done for me. vi TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 ISOGEOEMTRIC LOOP SUBDIVISION . . . . . . . . . . . . . . . 6 CHAPTER 3 MANIFOLD HARMONICS . . . . . . . . . . . . . . . . . . . . . . . 20 CHAPTER 4 ACOUSTIC THEORY AND INTEGRAL EQUATION FORMULATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 CHAPTER 5 AN EFFICENT ELECTROMAGNETIC ISO-GEOMETRIC INTEGRAL EQUATION SOLVER . . . . . . . . . . . . . . . . . . . . 61 CHAPTER 6 SHAPE OPTIMIZATION . . . . . . . . . . . . . . . . . . . . . . . . . 90 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 vii CHAPTER 1 INTRODUCTION 1.1 Introduction Computational methods have seen a significant investment in their development across numerous disciplines over the past few decades. This period has witnessed development of both higher order basis functions that admit accurate discretizations for numerical solution spaces, higher order representations of geometry for accurate representation of geometries (at least locally), and lastly the development of acceleration methodologies to reduce the computational cost of solving integral equations for large scattering bodies, amongst many other equally important advances. In this period, the first two advancement have largely been disconnected leading to a host of complications. The standard paradigm for computational analysis proceeds as follows: (i) taking a CAD model (higher order differentiable representation of the geometry), (ii) generate an analysis ready mesh, and (iii) post-process. This process introduces a host of complications including, but not limited to, geometric approximation issues and large computational overhead. Moreover, it is not suitable for design-through-analysis problems, i.e. design optimization. The reason for these complications are due to the fundamental distinction between the geometry processing and analysis based on this geometry. Isogeomtric analysis bridges the gap between the two technologies by using the same basis to represent both the geometric as well as the physics on the geometry. 1.2 Challenges arising in an Isogeometric Framework In standard frameworks, the process begins by reducing the exact geometry natively defined using CAD tools to a surface mesh, or an approximation of the exact geometry. This approximation can lead to a host of complications including geometric error and lack of mesh adaptivity leading to analytically incorrect answers and high-precision results are impossible. The advant of (IGA) based methods available to at least partly solve the 1 problem [79]. Considering that the choice of physical function representation is usually strongly related to the choice of geometric representation. A fundamental step is to use a geometric model that can be utilized directly as an analysis model. This deviates from the standard framework by using the CAD tools for analysis; choosing a good geometry technology is necessary to reduce the numerical error as it embeds the native geometry information in the basis design for physics; this implies that IGA is a strong candidate for physical simulations. Several aspects, at least, of IGA should be developed further for future industry applications. In this body of work we address the following: (1) the problem of efficient numerical integration, (2) representing vector fields on multiply connected surfaces, (3) further compressed representation of the physics and the geometry and its application to (4) acoustic and electromagnetic scattering and (5) shape reconstruction. 1.2.1 Representation of physical quantities on manifolds The accurate representation of equivalent current and pressure densities on scattering boundaries are critical for analysis. Hence, in this body of work, we develop a representation of a complete Helmholtz decomposition on a multiply connected Loop subdivision surface. On simply connected surfaces, this permits the definition of basis sets that are divergence- free and curl-free. What is missing from this suite is a basis set that is both divergence-free and curl-free, a necessary ingredient for a complete Helmholtz decomposition of currents on multiply-connected structures. Hence, we effect this missing ingredient numerically using random polynomial vector fields. We show that this basis set is analytically divergence-free and curl-free. Furthermore, we show that these basis recovers curl-free, divergence-free, and curl-free and divergence-free fields. Finally, we use this basis set to discretize a well- conditioned integral equation for analyzing perfectly conducting objects and demonstrate excellent agreement with other methods. 1.2.2 Computational cost of inner products In constructing a Loop subdivision IGA boundary integral formulation our basis function is a 4๐‘ก โ„Ž order polynomial (second order smooth) with a large support domain. 2 When we discretize this integral equation using a Galerkin prescription we find that we must compute an integration, or inner product, over a very large support domain. To help visualize challenges in evaluating the necessary inner products, it should be noted that unlike Lagrangian discretization, the domain of a subdivision basis function is significantly larger; it includes the union of incident triangles. Since we use the same basis function for physics as well, it implies that one needs to evaluate products of higher order polynomials convolved with a Greenโ€™s function on higher order surfaces, for all pairwise interactions as is dictated by the Surface Integral Equation (SIE) formulation. This limits the applicability to larger scatters and as such the need for an acceleration technique is imminent. Therefore, in this thesis we develop such a technique that specifically handles this bottleneck that arises in both acoustic and electromagnetic integral equations. In doing so, we shall propose and demonstrate an efficient method based on a wideband fast multipole method to evaluate the required inner products and matrix vector products. 1.2.3 Manifold Harmonics The manifold harmonic basis (MHB), a smooth, orthogonal global basis set, give us a framework wherein we can apply many signal processing operators, including filtering, compression, and so on, to both the forward problem as well as the shape optimization problem. Moreover, given that we define them on the Loop subdivision surface, it discrete representation inherits properties that are of interest when representing physical quantities on a surface: high-order continuity and multiresolution. That is, one can compose a MHB on the limit surface via a composition of those computed on a series of refined control meshes form the hierarchical manifold harmonic basis (H-MHB). While this does provide an excellent representation on objects that are sufficiently smooth, the number of basis sets can be high if on desires higher fidelity near singular regions. But it is possible to incorporate localised enrichment based on known local physics or heuristics, in tandem with the standard basis set known as localized manifold harmonic basis (L-MHB). Here we explore development of both these types of basis sets and their application to problems in 3 acoustics and electromagnetics. 1.2.4 Shape Reconstruction The last avenue we explore is the recovery of the shape of scattering obstacles from phaseless far-field data, otherwise referred to as the shape reconstruction problem. The main challenge we aim to address is minimizing the search space, or the number of degrees of freedom in this problem. The solution proceeds by parameterizing some initial guess and then finding the optimal paramteric coefficents that minimize a cost functional that measures the discrepancy between the prescribed far-field data and the far-field pattern corresponding to the current approximation of the paramterized scatterer. As is apparent, the geometry parameterisation and its interplay with the analysis plays a crucial role in shape optimisation/reconstruction. When a geometric representation, distinct from the computer-aided design (CAD) model is used to represent the shape, additional errors are incurred as well as a large parameterisation space is needed to represent high fidelity shapes, drastically increasing the difficulty of the optimization scheme. In addition, it leads to non-physical oscillations in the optimised geometry [75] as well as severely distorted mesh requiring auxiliary mesh smoothing, further complicating the problem. To remedy these difficulties, the shape optimisation/reconstruction can be constructed in a high-order isogeometric framework. 1.3 Scope and Outline of the Thesis The thesis is structured as follows. The second chapter of the thesis details the implementation of discretization of physical quantities on simply connected and multiply connected Loop subdivision surfaces. This includes expanding on the multiresolution features of Loop subdivision, as well as numerical examples to show the effectiveness of the decomposition. The third chapter employs the MHB. MHB allows us to associate a frequency spectrum to a function on a manifold, analogous to the Fourier decomposition. This insight can be used to build a framework for analysis. The purpose of this chapter is to 4 review and illustrate such possibilities for computational acoustics and electromagnetics as well as chart a potential path forward. To this end, we introduce three features of MHB: (a) enrichment for analysis of multiply-connected domains, (b) local enrichment (L-MHB) and (c) hierarchical representation (H-MHB) for reuse of data from coarser to fine geometry discretizations. Several results highlighting the efficacy of these methods are presented. The fourth chapter details the implementation of an IGA acoustics integral equation that address the challenge of integration associated with the higher order nature of our basis set by using an adaptive quadrature rule that is accelerated via MLFMA. We provide several numerical examples demonstrating the validty of the approach. In the fifth chapter we implement a fast and well-conditioned boundary integral equation based isogeometric framework that governs electromagnetic scattering phenomena from multiply-connected Loop subdivision surfaces. This new framework provides us better accuracy, more stability and high flexibility. Finally, the last chapter lies in utilizing manifold harmonics on a Loop subdivision surface for the shape reconstruction problem. These MHs provide the means to smoothly interpolate data on a manifold and are highly effective, specifically as it relates to geometry representation and editing; MHB form a natural basis for multi-resolution representation (and editing) of complex surfaces and functions defined therein. 5 CHAPTER 2 ISOGEOEMTRIC LOOP SUBDIVISION 2.1 Introduction Computational methods have become critical for guiding scientific investigation in numerous disciplines. In the past few decades, considerable effort have been invested into integral equation and differential equation based methods. The efficacy of the solution to these equations depends on the discrete representation of (a) the geometry and (b) the numerical solution space. The development of both higher order basis functions [135, 83] and higher order representations of geometry (at least locally) [71, 37] have been among the many other important advances. However, despite advances in these areas, the technologies of geometric design and numerical simulation have been fundamentally disjoint. Traditional analysis proceeds by selecting a discrete representation of the โ€œexact" CAD geometry [79] via a piecewise continuous tessellations and a different discrete representation for the numerical simulation. This geometric approximation yields an awkward communication with the CAD software for refining and remeshing which greatly diminishing the accuracy of the model. This is especially highlighted when considering the geometrical representation of the computational domain and for the evaluation of the geometric quantities; lack of higher order continuity in geometry can introduce errors in the normal and the curvatures, as well as tangential differential operators. Therefore, particular care must be devoted to the choice of a suitable geometric representation and numerical space discretization methods, as the error introduced by an approximation of the geometry can affect the accuracy of the numerical solutions; this could also potentially lead to the adoption of inefficient computational meshes, requiring large amounts of Degrees Of Freedom (DOFs), in order to control the geometrical error. One remedy is Isogeometric Analysis (IGA), a means of unification of the the disjoint 6 technologies of geometric design and numerical simulation. This is achieved by using the same representation for the geometry of the computational domain as well as for the numerical solution [119]. The advantage of the isogeometric framework is the elimination of both the geometry approximation error and the computational cost of changing the representation between design and analysis. Ideally, this means one circumvents the costly operations, in traditional analysis, of building an approximation of the CAD model that is suitable for numerical analysis. Furthermore, avoiding geometry representation reduction explicitly can also improve the efficiency of the design-through-analysis process including optimization problems, uncertainty quantification, verification and validation, analysis of deformable scatters, etc. The traditional workhorse representation tool of IGA methods has been non-uniform rational B-splines (NURBS). NURBS use bi/tri-variate spline based patches/solids such as those based on Bezier, B-splines, and are the the most ubiquitous CAD tools. NURBSโ€™ properties such as non-negativity and the fact that it provides a partition of unity make it an excellent candidate for defining function spaces. The choice of representation of the physics is usually strongly related to the choice of geometric representation which necessary to mitigate overall errors. Finite element methods based on NURBS basis functions that exhibit h- and p-adaptivity have been demonstrated [79]. Unfortunately, NURBS faces difficulty when the object being meshed is topologically complex or is multi-scale [55, 24], which can lead to surfaces that are not watertight. Some other popular IGA basis candidates that are ideal for handling shapes that are complex are T-splines, THB-splines, and subdivision surface. T-splines greatly reduce the number of the control points in the control mesh. Analysis-ready T-splines are a good candidate for constructing isogeometric analysis. More detailed work on T-splines and its application in IGA can be found in [24] and references therein. THB-splines [31], properly satisfy the requirements needed for developing not only effective geometric modeling tools, but also efficient computational schemes that provides local refinement possibilities 7 [31]. The varying strengths and weaknesses for these methods can be found in [63]. In contrast to the rest, subdivision surfaces have played a significant role in the computer animation industry. Among its many advantages are the ease with which one can represent complex topologies, scalability, inherent multiresolution features, efficiency and ease of implementation. Furthermore, it converges to a smooth limit surface that is ๐ถ 2 almost everywhere except at isolated points where it is ๐ถ 1 continuous. In the next section, we focus on subdivision schemes. 2.2 Subdivision Surfaces Subdivision surfaces are a representation of a smooth limit surface by a control mesh. This control mesh can be refined by various schemes, and in the limit of infinitely many refinements it approaches the limit surface. In practice these refinements are repeated until the surface is sufficiently fine. This is the approach used in applications such as computer graphics and visualization, where subdivision surfaces are most commonly employed. The properties of the limit surface depend on the type of control mesh used and the choice of refinement scheme. Subdivision schemes come in several different flavors; among the most popular are Loopโ€™s scheme on triangular meshes [107], and Catmull-Clark [38] and Doo-Sabin [57] schemes on quadrilateral meshes. Generally speaking, all of the three of these schemes alluded to above can be used to construct IGA methods. The subdivision approach has been used in thin-shell deformation problems [44] and a boundary element method for electrodynamics [148], some work on IGA based on Catmullโ€“Clark can be found in [22], and where IGA is used to solve PDEs defined on a surface [85, 121, 101]. In the remaining portions of this work we focus on Loop subdivision basis functions have been extensively studied, most interestingly regarding their smoothness [131, 132], (local) linear independence [123], approximation power [15], and the robust evaluation around so called extraordinary vertices [149]. 8 2.2.1 Loopโ€™s Subdivision Scheme Let ๐‘‡ ๐‘™ denote a ๐‘™-th refined control mesh, with vertices ๐‘‰ ๐‘™ := {v๐‘™๐‘– , ๐‘– = 1, . . . , ๐‘๐‘ฃ๐‘™ } and triangular faces ๐‘ƒ ๐‘™ := {p๐‘™๐‘– , ๐‘– = 1, . . . , ๐‘ ๐‘™๐‘“ }. In short, we can represent a ๐ถ 2 (almost everywhere) smooth limit surface ฮ“, through an infinite number of iterative refinements of the control mesh ๐‘‡ 0 , following the loop subdivision scheme [170]. Each Loop subdivision step involves dividing each triangle into four subtriangles by inserting one vertex per edge, and adjusting vertex locations. The location of the newly inserted vertex is an affine combination of the end vertices of the edge, each with weight 3/8, and the other two vertices of the two incident triangles, each with weight 1/8. The location of the existing vertex on the refined mesh is an affine combination of its one-ring vertices with weight ๐›ผ ๐‘› = (5/8โˆ’(3/8+ ๐‘๐‘œ๐‘ (2๐œ‹/๐‘›)/4)โˆ—2)/๐‘›, and itself with weight 1โˆ’ ๐‘›๐›ผ ๐‘› , where ๐‘› is the valence; both the old and new vertex positions are required to update during each refinement [106, 38]. In the limit of an infinite number of refinements of the control mesh, we approach the limit surface with ๐ถ 2 smoothness almost everywhere, except at the extraordinary vertices, where it has ๐ถ 1 smoothness. In practice, this prescription is not followed. There exists closed form expressions for computing the limit surface ฮ“ for a given control mesh ๐‘‡ ๐‘™ in terms of quantities defined on the given control mesh [149]. Assume that a subdivision surface admits a natural parameterization of the surface ฮ“ in terms of the barycentric coordinates defined on each face ๐œ– โˆˆ ๐‘ƒ ๐‘™ , for some ๐‘™. We begin by considering any patch ๐œ– โˆˆ ๐‘ƒ ๐‘™ for some ๐‘˜, as depicted in Fig. 2.1. We define the 0-ring of a patch (triangle) as the vertices that belong to the patch, and the 1-ring as the set of all vertices, ๐‘›๐‘ฃ , that can be reached by traversing no more than two edges, as shown in Fig. 2.1. We define the regularity of the triangle by the characterization of its verticesโ€™ valence (0-ring); the valence of a given vertex is the number of edges incident on itself. A vertex is considered regular if its valence is equal to 6, otherwise, it is called an irregular or extraordinary vertex. A triangle is regular if its vertices are all regular, and irregular otherwise. Using these definition, we can define the 9 mapping from the barycentric coordinates on a given patch, ๐œ–, to the limit surface by a weighted average of the effective basis functions associated to its 1-ring [149]. As a result, we can define the limit surface as ๐‘๐‘ฃ๐‘™ ร• ฮ“(r) = c๐‘›๐‘™ ๐œ‰๐‘›๐‘™ (r), (2.1) ๐‘›=1 where c๐‘›๐‘™ are vertex locations of the ๐‘๐‘ฃ๐‘™ control points, and ๐œ‰๐‘›๐‘™ is the effective basis function that is associated with quantities associated with c๐‘›๐‘™ and has a support ฮ“๐‘›๐‘™ ; note, โˆช๐‘› ฮ“๐‘›๐‘™ = ฮ“๐‘™ . We note that we have only provided a brief overview of Loop subdivision as an isogeometric tool; information provided is purely for completeness and omits details that can be found in [44, 149, 106, 97, 6, 121, 85] and references therein. The basis functions {๐œ‰๐‘›๐‘™ } span a IGA finite dimensional space ฮจ that is the subspace of the Sobolev space ๐ป 2 (ฮ“) [121, 85]. Its properties are (a) positivity, (b) compact support, (c) forming a partition of unity, and (d) ๐ถ 2 continuity almost everywhere. Furthermore, the hierarchy of control meshes underlying a subdivision surface lends itself naturally to multiresolution analysis [172, 108]. Figure 2.1: Regular triangular patch defined by its 1-ring vertices. 10 2.2.2 Multiresolution subdivision surfaces Subdivision surfaces represent a limit surface associated with a nested hierarchy of control meshes of increasing resolution. This leads to a natural representation at various levels of detail which is particularly valuable for various applications, in both computer graphics and compression of functions defined on surfaces, numerical solution of integral and differential equations, because the refinement rules can be applied adaptively. This multiresolution analysis can be applied by any uniformly convergent subdivision scheme [108]. Each application of a subdivision rule applied to a surface generates a refined surface of increased detail. A sequence of refinements produces a hierarchy that enables the selection of a representation of a surface at a particular level of detail. The basic idea behind multiresolution analysis is to decompose a complicated function into a โ€œsimplerโ€ low resolution part and a โ€œdetailโ€ high resolution part necessary to recover the original function. In order to exploit the hierarchical refinement, we must define these prolongation and restriction operators that allow us to map data defined on one mesh to another mesh. It then becomes natural to consider the relation between functions on multiple levels. The resulting algorithms allow the concurrent editing of coarse and fine levels, in multiresolution analysis of surfaces the refinement and coarsening matrices [๐‘†] and [๐‘…], respectively, are used to construct a decomposition of the geometry; we define these operators in the next section. 2.2.2.1 Prolongation and Restriction Operators The Loop subdivision refinement process needed to generated a family of meshes can be represented using matrix operations because the subdivision rules are linear combinations of the mesh vertices; note, the subdivision matrix is uniquely defined by the subdivision level, the number of control maps and its connectivity map [60]. As such, we encode the weights of a loop subdivision scheme for the vertices at level ๐‘™ as as sparse matrix ๐‘™+1 ร—๐‘ ๐‘™ ๐‘† ๐‘™ โˆˆ R ๐‘๐‘ฃ ๐‘ฃ , [144, 53, 165]. Therefore, ๐‘‰ ๐‘™+1 = ๐‘† ๐‘™ ๐‘‰ 0 and, where ๐‘† ๐‘™ = ๐‘† ๐‘™โˆ’1 ร— ๐‘† ๐‘™โˆ’2 ร— ยท ยท ยท ร— ๐‘†0 represents the weighting rule of the subdivision mesh ๐‘‡ ๐‘™ generated from the initial control 11 mesh ๐‘‡ 0 , obtained by sparse matrix multiplication. As ๐‘™ โ†’ โˆž, we approach the position of any vertex on the limit surface. The prolongation operator is defined as [๐‘‰ ๐‘™+1 ] = [๐‘† ๐‘™ ][๐‘‰ ๐‘™ ], (2.2) and the coarsening operator is defined as [๐‘‰ ๐‘™ ] = [๐‘… ๐‘™ ][๐‘‰ ๐‘™+1 ]. (2.3) [๐‘… ๐‘™ ] enables the computation of the coarse representation [๐‘‰ ๐‘™ ] corresponding to a given fine representation [๐‘‰ ๐‘™+1 ]. This matrix can be determined by the least squares fitting, which leads [๐‘… ๐‘™ ] = ([๐‘† ๐‘™ ]๐‘‡ [๐‘†])โˆ’1 ๐‘†๐‘‡ . Notice that [๐‘… ๐‘™ ] is the pseudoinverse of [๐‘† ๐‘™ ] such that [๐‘… ๐‘™ ][๐‘† ๐‘™ ] = [๐ผ] yields the identity matrix; we note, different choices for the coarsening matrix are possible. 2.3 Iso-geometric Basis Sets 2.3.1 Loop subdivision Basis Functions To discretize and solve integro-differential equations on closed subdivision surfaces we consider isogeometric function spaces defined via Loopโ€™s subdivision. To define this isogeometric basis set, we assume that, again, there exists a net of control function values, coincident with the location of the control net. Thus, any scalar function ( ๐‘“ (r)) can then be expressed in terms of the Loop subdivision basis set via, ๐‘™ ๐‘๐‘ฃ ร• ๐‘“ (r) = ๐‘Ž ๐‘›๐‘™ ๐œ‰๐‘›๐‘™ (r), (2.4) ๐‘–=๐‘› where ๐‘๐‘ฃ๐‘™ and ๐œ‰๐‘›๐‘™ (r) retain the same definition as those prescribed above. The properties of this representation follow from those for subdivision. Henceforth, the functions ๐œ‰๐‘›๐‘™ (r) will be referred to as Loop basis. The aforementioned properties are critical to the development of isogeometric analysis. 12 2.3.1.1 Complete Helmholtz decomposition on Multiply connected Loop Subdivison Surfaces This Loop basis can also be used to define vector quantities on the surface as well. To do so, we begin by representing currents on any closed surface ฮ“, via the Helmholtz decomposition as J(r) = โˆ‡ฮ“ ๐œ™(r) + โˆ‡ฮ“ ร— (nฬ‚๐œ“(r)) + ๐œ”(r), ยฏ (2.5) where ๐œ”(r) ยฏ is the harmonic field and has a trivial solution, (๐œ”(r) ยฏ = 0), for simply connected structures, โˆ‡ฮ“ is the surface gradient, and ๐œ“(r) and ๐œ™(r) are scalar potentials. In what follows, we construct currents in terms of the scalar potentials using the loop subdivision basis sets. Using (2.4) we can define the scalar potentials ๐œ™(r) and ๐œ“(r) on the limit surface as ๐‘๐‘ฃ ๐‘™ ร• ๐‘™ หœ = ๐œ™(r) โ‰ˆ ๐œ™(r) ๐‘Ž1,๐‘› ๐œ‰๐‘›๐‘™ (r), ๐‘›=1 ๐‘™ (2.6) ๐‘๐‘ฃ ร• ๐‘™ หœ = ๐œ“(r) โ‰ˆ ๐œ“(r) ๐‘Ž2,๐‘› ๐œ‰๐‘›๐‘™ (r). ๐‘›=1 ๐‘™ Here ๐‘๐‘ฃ๐‘™ , ๐œ‰๐‘›๐‘™ (r) are defined earlier and {๐‘Ž ๐‘˜,๐‘› } for ๐‘˜ โˆˆ {1, 2} are the unknown coefficients. The curl-free and divergence-free components of the current can be represented using this approximation of the potentials, viz., ๐‘™ ๐‘๐‘ฃ h ร• i J(r) โ‰ˆ J๐‘™๐‘ (r) = ๐‘™ ๐‘Ž 1,๐‘› ๐‘™ J1,๐‘› (r) + ๐‘™ ๐‘Ž 2,๐‘› ๐‘™ J2,๐‘› (r) ยฏ ๐‘™ (r), +๐Ž (2.7a) ๐‘›=1 ๐‘™ J1,๐‘› (r) = โˆ‡ฮ“ ๐œ‰๐‘›๐‘™ (r), ๐‘™ J2,๐‘› (r) = nฬ‚(r) ร— โˆ‡ฮ“ ๐œ‰๐‘›๐‘™ (r), (2.7b) ร•๐‘” h i ๐Žยฏ ๐‘™ (r) = ๐‘™ ๐‘Ž 3,๐‘› ๐‘™ J3,๐‘› (r) + ๐‘Ž4,๐‘› ๐‘™ ๐‘™ J4,๐‘› (r) , ๐‘›=1 ๐‘™ ๐‘™ wherein we define J3,๐‘› (r) and J4,๐‘› (r) in later sections. Note that since the representation is constructed using conditions on currents that rely on derivatives of the potentials ๐œ™(r) หœ 13 หœ and ๐œ“(r), this leads to the existence of nontrivial solutions to (2.7) and therefore we must enforce uniqueness. In order to ensure uniqueness, we impose a zero-mean constraint for both divergence-free and curl-free components. 2.3.1.2 Representation of harmonic components Next, we detail the representation of the harmonic component. We begin by highlighting particular features of the Loop subdivision representation of the current thus far; more details can be found in [98, 97]. They are as follows: 1. The representation of currents is ๐ถ 1 continuous on ฮ“(r) almost everywhere. 2. The normal is ๐ถ 1 almost everywhere. 3. On ฮ“๐‘›๐‘™ โˆฉ ฮ“๐‘š ๐‘™ , hJ๐‘™ (r), J๐‘™ (r)i = โˆ’hJ ๐‘™ (r), J ๐‘™ (r)i. 1,๐‘› 2,๐‘š 1,๐‘š 2,๐‘› ๐‘™ ๐‘™ 4. The inner product, hJ1,๐‘š (r), J2,๐‘› (r)i, is identically zero. This can be proven using Greenโ€™s theorems, ๐œ‰๐‘š ๐‘™ (r) = 0 for r โˆˆ ๐œ•ฮ“๐‘™ , and properties stated earlier. To wit, ๐‘š โˆซ ๐‘™ ๐‘™ ๐‘™ ๐‘™ hJ1,๐‘š (r), J2,๐‘› (r)i = โˆ‡ฮ“ ๐œ‰๐‘š (r) ยท J2,๐‘› (r)๐‘‘r ๐‘™ ฮ“๐‘š โˆฉฮ“๐‘›๐‘™ โˆซ ๐‘™ ๐‘™ = ๐œ‰๐‘š (r)uฬ‚(r) ยท J2,๐‘› (r)๐‘‘r (2.8) ๐‘™ ๐œ•(ฮ“๐‘š โˆฉฮ“๐‘›๐‘™ ) โˆซ ๐‘™ ๐‘™ โˆ’ ๐œ‰๐‘š (r)โˆ‡ฮ“ ยท J2,๐‘› (r)๐‘‘r = 0, ฮ“๐‘›๐‘™ โˆฉฮ“๐‘š ๐‘™ where ๐œ•(ฮ“๐‘š ๐‘™ โˆฉ ฮ“๐‘™ ) is the boundary of ฮ“๐‘™ โˆฉ ฮ“ ๐‘™ , and uฬ‚(r) is a outward pointing normal ๐‘› ๐‘š ๐‘› tangential to the boundary. Given the properties stated earlier, the integrals are all zero. 5. As a result, the basis functions defined earlier are exaclty orthogonal. ๐‘™ ๐‘™ ๐‘™ ๐‘™ ๐‘™ ๐›ฟ where, 6. Finally, we denote hJ1,๐‘š (r), J1,๐‘› (r)i = hJ2,๐‘š (r), J2,๐‘› (r)i = ๐›พ๐‘š๐‘› ๐‘™๐‘˜ โˆซ ๐‘™ ๐›พ๐‘š๐‘› = โˆ‡๐‘  ๐œ‰๐‘›๐‘™ (r) ยท โˆ‡๐‘  ๐œ‰๐‘š ๐‘™ (r)๐‘‘r. (2.9) ฮ“๐‘›๐‘™ โˆฉฮ“๐‘›๐‘™ 14 These properties can be exploited readily in constructing a purely numerical basis for the harmonic components. We build on the ideas presented in [4]. Prior to proceeding, we note that the space of harmonic fields is 2๐‘”-dimensional, where ๐‘” is the genus of the object. Consider a function F(r) = โˆ’nฬ‚(r) ร— nฬ‚(r) ร— V(r) where V(r) โˆˆ R3 and whose components are a low degree random polynomial. Specifically, let ฮ“(r) โˆˆ (0, ๐ฟ ๐‘ฅ ) ร— (0, ๐ฟ ๐‘ฆ ) ร— (0, ๐ฟ ๐‘ง ) .  Then each component of V(r) can be defined as a tensor product of Legendre polynomials via ๐‘ƒ๐›ผ (๐‘ฅ/2๐ฟ ๐‘ฅ )๐‘ƒ๐›ฝ (๐‘ฆ/2๐ฟ ๐‘ฆ )๐‘ƒ๐œ (๐‘ง/2๐ฟ ๐‘ง ) where ๐‘ƒ๐œ (ยท) is the standard Legendre polynomials of degree ๐œ defined on [โˆ’1, 1]. Moreover, the maximum degree of Legendre polynomials should be at least d(2๐‘”)1/3 e so that dimension of the space of polynomials is greater than the dimension of the space of the vector fields. For the examples in this paper, we choose the degree to be ๐›ผ + ๐›ฝ + ๐œ = d(2๐‘”)1/3 e + 5. To obtain ๐‘” harmonic basis functions that project to ๐Ž(r), ยฏ we start with the following: consider a collection of ๐‘” random vector fields formed using Legendre polynomials defined earlier and denoted using F๐‘– (r) for ๐‘– โˆˆ {1, ยท ยท ยท , ๐‘”}. For each such field F๐‘– (r), we can obtain ร•h i ๐‘™ ๐‘™ ๐‘™ ๐‘™ ๐‘™ J3,๐‘– (r) = F๐‘– (r) โˆ’ ๐‘Ž 1,๐‘›,๐‘– J1,๐‘› (r) + ๐‘Ž 2,๐‘›,๐‘– J2,๐‘› (r) (2.10) ๐‘› h i๐‘‡ h i๐‘‡ ๐‘™ ๐‘™ ๐‘™ ๐‘™ ๐‘™ ๐‘™ where the coefficients are arranged as โ„1,๐‘– = ๐‘Ž 1,1,๐‘– ,ยทยทยท , ๐‘Ž1,๐‘›,๐‘– and โ„2,๐‘– = ๐‘Ž 2,1,๐‘– ,ยทยทยท , ๐‘Ž2,๐‘›,๐‘– , and obtained by solving ๏ฃฎ 11 ๏ฃฎ ๐‘™ ๏ฃน ๏ฃฎ ๐‘™ ๏ฃน ๏ฃฏ๐’ข 0 ๏ฃบ ๏ฃฏโ„1,๐‘– ๏ฃบ ๏ฃฏ๐’ฑ1,๐‘– ๏ฃบ ๏ฃน ๏ฃฏ ๐‘™ ๏ฃบ ๏ฃฏ ๐‘™ ๏ฃบ. (2.11a) ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ=๏ฃฏ ๏ฃบ ๏ฃฏ 0 ๐’ข 22 ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏโ„2,๐‘– ๏ฃบ ๏ฃฏ๐’ฑ2,๐‘– ๏ฃบ ๏ฃฐ ๏ฃป ๏ฃฐ ๏ฃป ๏ฃฐ ๏ฃป Here, ๐’ข๐‘›๐‘š11 = ๐’ข 22 = ๐›พ ๐‘›๐‘š ๐‘›๐‘š , and โˆซ ๐‘™ ๐’ฑ๐‘˜,๐‘›,๐‘– = J๐‘™๐‘˜,๐‘› (r) ยท F๐‘– (r)๐‘‘r. (2.11b) ฮ“๐‘›๐‘™ ๐‘™ ๐‘™ ๐‘™ Furthermore, note that if J3,๐‘– (r) is a harmonic vector field, then J4,๐‘– (r) = nฬ‚(r) ร— J3,๐‘– (r) is also a ๐‘™ harmonic vector field that is independent of J3,๐‘– (r). These basis span the space of harmonic fields (with high probability that the set of random vector fields are linearly independent) [4]. 15 There are several properties that make this decomposition viable for analysis and as such we examine these properties below. Specifically, we note the following: ๐‘™ 1. The inner products of a field F(r) with basis J1,๐‘› ๐‘™ (r) and J2,๐‘› (r), for ๐‘› โˆˆ {1, . . . , ๐‘๐‘ฃ๐‘™ }, effectively measures the surface divergence and curl of the field. 2. It it trivial to show that both the surface divergence and the surface curl of J3,๐‘› ๐‘™ (r) and ๐‘™ J4,๐‘› (r), for ๐‘› โˆˆ {1, . . . , ๐‘”} as defined in (2.10), are identically zero via (2.11). 3. Along the same lines, it can be shown that hJ๐‘™๐‘˜,๐‘› (r), J๐‘™๐‘–,๐‘› (r)i = 0 for ๐‘˜ โˆˆ {1, 2} and ๐‘– โˆˆ {3, 4}. ๐‘™ ๐‘™ 4. The support of J3,๐‘– (r) and J4,๐‘– (r) is global and as such we need additional infrastructure to evaluate the necessary matrix elements. A direct consequence of these properties and those in section 2.3.1.1 is that the prescribed basis J๐‘™๐‘˜,๐‘› (r) for ๐‘› โˆˆ {1, . . . , ๐‘๐‘ฃ๐‘™ } and ๐‘˜ โˆˆ {1, 2}, and J๐‘™๐‘ž,๐‘– (r) for ๐‘ž โˆˆ {3, 4} and ๐‘– โˆˆ {1, ยท ยท ยท , ๐‘”}, provide a complete Helmholtz decomposition of vector quantities defined on the manifold. 2.3.1.3 Helmholtz Decomposition Example In this example, we demonstrate the effectiveness of our prescribed basis set to affect a Helmholtz decomposition of a vector field that is a linear combination of analytically prescribed curl-free, divergence-free, and harmonic vector fields, defined on a torus. We choose a torus, see Fig. 2.2, with a major axis of 3 m and a minor axis of 1 m, modeled using an initial control mesh comprising of 2048 vertices and 4096 faces. We refine the control mesh twice such that surface area and the total curvature of the limit surface ฮ“ agree within 99% to the analytical torus. We aim to report the ๐ฟ2 norm of each component of the Helmholtz decomposition, for 16 Figure 2.2: The two basis harmonic vector fields for a torus. the following vector field Xฬ„(r). In particular, we have X๐‘‘ (r) = โˆ’nฬ‚(r) ร— nฬ‚(r) ร— r, X๐‘ (r) = nฬ‚ ร— X๐‘‘ (r), t (2.12) X โ„Ž (r) = โˆ’nฬ‚(r) ร— nฬ‚(r) ร— , ||t|| 2 Xฬ„(r) = X๐‘‘ (r) + X๐‘ (r) + X โ„Ž (r) + nฬ‚(r) ร— X โ„Ž (r) where t = ๐‘ฆ ๐‘ฅห† โˆ’ ๐‘ฅ ๐‘ฆ, ห† and X๐‘ is a divergence-free field, X๐‘‘ is a curl-free field, and X โ„Ž and nฬ‚ ร— X โ„Ž are purely harmonic fields. Next, we choose to approximate these fields using the basis sets defined earlier; specifically, Xฬ„(r) โ‰ˆ J1 (r) + J2 (r) + J3 (r) + J4 (r) ร• ร• = ๐‘Ž 1,๐‘› J1,๐‘› (r) + ๐‘Ž 2,๐‘› J2,๐‘› (r) ๐‘› ๐‘› (2.13) ร•๐‘” ร• ๐‘” + ๐‘Ž3,๐‘› J3,๐‘› (r) + ๐‘Ž 4,๐‘› J4,๐‘› (r), ๐‘›=1 ๐‘›=1 where, J1 (r) โ‰ˆ X๐‘‘ (r), J2 (r) โ‰ˆ X๐‘ (r), J3 (r) โ‰ˆ X โ„Ž (r), and J4 (r) โ‰ˆ nฬ‚(r) ร— X โ„Ž (r). Note, we drop the superscript that denotes the level of refinement for brevityโ€™s sake. We compute the unknown coefficients using the orthogonality properties of the basis set as presented in sections 2.3.1.1 and 2.3.1.2. Table. 2.1 lists the ๐ฟ2 -norm of each component of the decomposition of the vector fields, respectively. Furthermore, in Fig. 2.2 we plot the two resulting harmonic fields. Ideally, the four components are recovered in the corresponding discrete spaces. As expected, the dominant components in each decomposition are reflected in the correct subspaces. The 17 X๐‘‘ X๐‘ Xโ„Ž nฬ‚ ร— X โ„Ž k J1 k2 5.05E2 2.28E-17 3.02E-10 3.31E-19 k J2 k2 2.28E-17 5.05E2 3.31E-19 3.02E-10 k J3 k2 1.01E-26 1.67E-26 1.36E2 5.26E-18 k J4 k2 1.98E-29 6.87E-24 5.02E-17 1.36E2 Table 2.1: ๐ฟ2 -norms of the components of the decomposed vector fields. projection on to other subspaces are significantly smaller in magnitude, but do not vanish. This is largely due the fact that while (2.12) is analytically exact on a canonical torus, it is not so on a geometric approximation. 2.4 Conclusion In this chapter, we explore Loop subdivision as a tool for IGA, a method that uses the same basis sets for both geometry representation as well as representation of physics on this geometry, which inherits the rather significant advantages of modern CAD toolsโ€“smoothness of geometry representation, morphing, dynamic meshes. Furthermore, this work develops the numerical infrastructure onto existing Loop subdivision IGA methods to enable analysis of vector fields on multiply connected structures. As a result, we are able to prescribe a complete Helmholtz decomposition on the surface wherein orthogonality of different components is preserved. Numerical examples are presented for validation of the efficacy and excellent convergence performance as well. This chapter, ยฉ2022 IEEE, in part, is reprinted, with permission from A. M. A. Alsnayyan and B. Shanker, "A complete Helmholtz decomposition on multiply connected subdivision surfaces and its application to integral equations," IEEE Transactions on Antennas and Propagation on, August 2022, with slight modifications to fit in this dissertation. This chapter, ยฉ2022 IEEE, in part, is reprinted, with permission from A. M. A. Al- snayyan and B. Shanker, "Manifold Harmonics and their application to computational 18 electromagnetics," IEEE Journal on Multiscale and Multiphysics Computational Techniques on, August 2022, with slight modifications to fit in this dissertation. 19 CHAPTER 3 MANIFOLD HARMONICS 3.1 Introduction The Laplace-Beltrami operator (LBO) is a second-order differential operator which is widely used in the fields of computer graphics [144, 96, 113], acoustics [9], electromagnetics [6, 7, 5, 4, 69], and machine learning [26, 163]. In particular, the use of the spectrum and eigenfunctions of the LBO, otherwise referred to as manifold harmonic basis (MHB), are of interest in this body of work. On a manifold, the set of eigenfunctions of the LBO form a complete and orthonormal basis on the space of square-integrable functions. The eigenbasis of the manifold bears significant similarity with the classical Fourier series in Euclidean space, and similarly, are leveraged to define the manifold harmonic transform (i.e. the generalized Fourier transform on the manifold). As such, the MHB provides a unique, compact multi-resolution basis for spectral shape processing that is independent of the actual shape representation and parameterization. Several successful applications have been proposed that take advantage of these desirable properties, such as spectral geometry filtering, compression, and surface deformation [157, 168]. While the manifold harmonic transform has been commonplace in the compute graphics literature for an array of applications, to our knowledge it has not been utilized in analysis of physics on manifolds. Typically, obtaining MHB proceeds by finding the eigenfunctions of the graph Laplacian [168]. But, in this body of work, we define the basis on Loop subdivision surfaces, as it inherits the attractive properties that are of interest when representing physical quantities on a surface: high-order continuity and multiresolution. For instance, given that a hierarchy of subdivision meshes point to same limit surface, it follows that if a function is associated with control points, then this will follow the refinement rules of subdivision. That is, one can compose a MHB on the limit surface via a composition of those computed on a series 20 of refined control meshes. Of course, one needs the prolongation and restriction operators [144, 19] for communication of data defined on different levels of refinement. This opens the door to applications such as shape and topology optimization, machine learning, target recognition, and data fusion for inverse problems where one can utilize a sequence of meshes. Furthermore, the MHB provides a smooth, orthogonal global basis set. While this does provide an excellent representation on objects that are sufficiently smooth, the number of basis sets can be high if on desires higher fidelity near sharp features (say tips). In response, it is possible to incorporate localised enrichment based on known local physics or heuristics, in tandem with the standard basis set. The notion of MHB is powerful as is evident from its numerous applications in computer graphics and machine learning. It has been used extensively in these areas, including (a) to effect mesh compression [86], (b) shape analysis, as it is robust to shape deformation [136, 134], (c) can be integrated with a localized harmonic basis for better representation [113], and (d) can be used to map data across meshes of different resolution [144]. Given the diverse advantages, the question we ask in this chapter is whether some of these advantages can be used for analysis. There has been an attempt to do just that for compression and shape optimization [9, 8]; the former for acoustics analysis and the latter for electromagnetic analysis. In this chapter, we explore that avenue several steps further. Specifically, we extend the notion of MHB to (a) multiply connected domains for electromagnetic analysis, (b) to local manifold harmonic bases (L-MHB) to better represent local rapidly varying fields, and (c) hierarchical manifold harmonic bases (H-MHB) to re-use data from a series of geometries of different resolutions. 3.2 Laplace-Beltrami Operator and the Manifold Harmonic Basis The LBO, ฮ”ฮ“ , is a self-adjoint and semi-positive definite operator [137] defined on a compact 2D Riemannian manifold ฮ“ embedded in R3 as follows, ฮ”ฮ“ ๐œ’(r) := โˆ‡ ยท (โˆ‡๐œ’(r)), (3.1) where ๐œ’(r) is a real-valued function defined on ฮ“. The LBO admits a complete and 21 countable sequence of eigenfunctions which form an orthonormal basis in ๐ฟ2 (ฮ“) [122], denoted by {๐ป๐‘š } such that โˆ’ ฮ”ฮ“ ๐ป๐‘š = ๐œ†๐‘š ๐ป๐‘š . (3.2) Here, {๐œ†๐‘š }โˆž๐‘š=0 is the corresponding spectrum. These eigenvalues constitute a countable set of non-negative real eigenvalues 0 = ๐œ†1 โ‰ค ๐œ†2 โ‰ค ยท ยท ยท [134]. Discrete computation: In order to numerically compute MHBs, we employ the Loop Subdivision FEM Galerkin method. This is akin to similar efforts using Lagrangian surface descriptions [133, 134] that have shown both โ„Žโˆ’ and ๐‘โˆ’ convergence [133, 67, 134]. The numerics necessary for computing eigenfuctions of the LBO relies on casting the Laplacian eigenvalue problem in a variational setting. The solution of this variational problem is approximated using the finite element Galerkin technique on the surface. We begin by evaluating an inner product of (3.2) with some test function ๐‘ฃ(r) โˆˆ {๐œ‰๐‘– (r)} and then use Greenโ€™s theorems to arrive to the following: hโˆ‡๐‘  ๐‘ฃ(r), โˆ‡๐‘  ๐ป๐‘š (r)i ฮ“ = โˆ’๐œ†๐‘š h๐‘ฃ(r), ๐ป๐‘š (r)i ฮ“ . (3.3) โˆซ where ๐‘“ (r), ๐‘”(r) ฮ“ = ฮ“ ๐‘“ (r) ยท ๐‘”(r)๐‘‘r follows the standard inner product definition. The MH ๐ป๐‘š (r) is represented in the same fashion as (2.4) leading to ร• ๐‘๐‘ฃ ๐‘™ ๐‘™ ๐ป๐‘š โ‰ˆ ๐ป g ๐‘š (r) = โ„Ž ๐‘–,๐‘š ๐œ‰๐‘– (r), (3.4) ๐‘– ๐‘™ for โ„Ž ๐‘–,๐‘š โˆˆ R. This leads to a generalized eigenvalue problem [๐ด ๐‘™ ][๐ป ๐‘™ ] = โˆ’[ฮ›][๐ต ๐‘™ ][๐ป ๐‘™ ], (3.5) where, โˆซ ๐‘™ ๐ด โˆ‡๐‘  ๐œ‰๐‘–๐‘™ (r) ยท โˆ‡๐‘  ๐œ‰ ๐‘™๐‘— (r)๐‘‘r,   ๐‘–๐‘— = (3.6a) ฮ“๐‘™๐‘– โˆซ  ๐‘™ ๐ต ๐‘–๐‘— = ๐œ‰๐‘–๐‘™ (r)๐œ‰ ๐‘™๐‘— (r)๐‘‘r. (3.6b) ฮ“๐‘™๐‘– 22 For this generalized symmetric eigenvalue problem [๐ด ๐‘™ ] โˆˆ R๐‘๐‘ฃ ร—๐‘๐‘ฃ is positive semi-definite, [๐ต ๐‘™ ] โˆˆ R๐‘๐‘ฃ ร—๐‘๐‘ฃ is positive definite, [ฮ›๐‘™ ] โˆˆ R๐‘๐‘ฃ ร—๐‘๐‘ฃ contains ๐‘๐‘ฃ eigenvalues along its diagonal, and [๐ป ๐‘™ ] โˆˆ R๐‘๐‘ฃ ร—๐‘๐‘ฃ contains the solution vectors, i.e. the coefficients of each eigenvector defined in (3.4), in its column space. For this symmetric generalized eigenvalue problem we have [๐ป ๐‘™ ]๐‘‡ [๐ด ๐‘™ ][๐ป ๐‘™ ] = [ฮ›๐‘™ ] and [๐ป ๐‘™ ]๐‘‡ [๐ต ๐‘™ ][๐ป ๐‘™ ] = [๐ผ ๐‘™ ], where [๐ผ ๐‘™ ] is the identity matrix. From the previous relations, it follows that the eigenfunctions are orthogonal with respect to the [๐ต ๐‘™ ]-based scalar product (i.e., hH๐‘™๐‘– , H๐‘™๐‘— i[๐ต ๐‘™ ] = H๐‘™,๐‘‡ ๐‘– [๐ต ๐‘™ ]H๐‘™๐‘— ). The eigenvectors with corresponding eigenvalues can be calculated with a direct eigensolver or by using the efficient band-by-band computation method presented in [157]. There is a extensive body of literature on efficient computation of these functions, largely applied to computer graphics [168]. 3.2.0.1 Compressed Representation of the manifold Let ฮ“ be a 2D manifold embedded in R3 . Considering [๐ป ๐‘™ ], we seek a representation ๐‘๐‘ฃ๐‘™ ร• ฮ“(r) = ๐›ฝ ๐‘™๐‘– ๐ป๐‘–๐‘™ (r). (3.7) ๐‘–=1 It is apparent from the orthogonal property of MHB that the coefficients can be defined as ๐›ฝ ๐‘™๐‘– = hฮ“(r), ๐ป๐‘–๐‘™ (r)i๐ฟ2 (ฮ“) = ฮ“(r)๐ป๐‘–๐‘™ (r)๐‘‘r; this, in effect, is dubbed a manifold harmonic โˆซ ฮ“ transform (MHT) of the surface [157]. Furthermore, it is important to remember that each eigenfunction ๐ป๐‘–๐‘™ (r) corresponds to different spatial โ€œfrequencyโ€ on the manifold. Thus, the magnitude of ๐›ฝ ๐‘™๐‘– , dictate the relative importance of a given eigenfunction at a given spatial โ€œfrequencyโ€. This transformation makes it possible to define equivalent signal processing operations: this includes under/over sampling, filtering, multi-resolution, windowing, and so on. To wit, this representation makes it possible to enrich shapes (and functions defined on the shape) systematically. To illustrate, the effectiveness of this approach, consider the 23 manifold ฮ“(r) again. Equation (3.7) can be rewritten as ๐‘๐‘ฃ๐‘™ ร•๐‘€๐‘™ ร• ฮ“(r) = ๐›ฝ ๐‘™๐‘– ๐ป๐‘–๐‘™ (r) + ๐›ฝ ๐‘™๐‘– ๐ป๐‘–๐‘™ (r) . (3.8) ๐‘–=1 ๐‘–=๐‘€+1 | {z } | {z } Low Frequency High Frequency In this expression, we have designated some cutoff coefficient ๐›ฝ ๐‘™๐‘€ , or correspondingly, q ๐œ† ๐‘™๐‘€ to be a maximum cutoff โ€œfrequencyโ€. This effectively is a versatile tool for many applications including EM anaylysis and shape reconstruction that enables us to enrich data at different levels of fidelity. A demonstration of MHT is shown in Figure 3.1, which depicts the reconstruction of a geometry as we increase the number of MHs. As is evident from these figures, it does not take too many MHs to capture the general shape of the object. As we increase the number of MHs, we rapidly approach the true shape, i.e. we capture more localized features. Note, the total number of MHs in the original system is ๐‘๐‘ฃ = 5002. (a) Original model. (b) With 100 MHs. (c) With 250 MHs. (d) With 750 MHs. Figure 3.1: The statue of the girl in (3.1a) is reconstructed with increasing number of MHs in (3.1b)-(3.1d). 24 (a) ๐ป1 (b) ๐ป9 (c) ๐ป500 (d) ๐ป1500 Figure 3.2: A select few MHs of the jet airliner. (a) ๐ป1 . (b) ๐ป9 . (c) ๐ป500 . (d) ๐ป1500 . 3.2.0.2 Analysis of Scalar Fields using Manifold Harmonics In this section we aim to demonstrate the utility of the MHT in analysis of physics on manifolds. In particular, we take advantage of its ability to rigorously compress the system. In what follows, we illustrate some of the features of this approach. The MHB defined earlier can be utilized as an isogeomtric basis set enabling us to not only represent geometry, but as a basis for representing physical quantities defined on the geometry; we note that the optimality of MHs for signal approximation has been addressed [3]. The salient features of this representation are as follows: (a) subdivision-based MHs are ๐ถ 2 smooth allowing for an excellent basis for representing physics on the geometry; (b) the orthogonal property allows for hierarchical signal representation, de-noising, and compression; (c) low computational cost and storage overhead feasible for applications in analysis; and (d) subdivision-based MHs span the same space as ฮจ. Considering the aforementioned properties, we aim to construct a spectral repre- sentation of scalar functions defined on the surface using MHs. Following (3.7), we have ๐‘™ ๐‘๐‘ฃ ร• ฮฆ(r) = ๐›ผ ๐‘™๐‘– ๐ป๐‘–๐‘™ (r), (3.9) ๐‘–=1 wherein it is possible to choose an error threshold, ๐œ–, for representation of desired functions on the manifold and thereby, truncate the number of MHs used. Numerical Example 1: To assess the efficiency and accuracy of our MH spectral represen- 25 (a) ๐ป1 (b) ๐ป9 (c) ๐ป500 (d) ๐ป1500 Figure 3.3: A select few MHs of the bumpy cube. (a) ๐ป1 . (b) ๐ป9 . (c) ๐ป500 . (d) ๐ป1500 . tation method for both physics and geometry we consider a 3.21๐œ† ร— 3.21๐œ† ร— 3.04๐œ† bumpy cube in Fig. 6.5c represented using 2562 subdivision basis functions as our candidate object; the physics on this surface is represented using an increasing number of MHs, and solution obtained is compared against an isogeometric solution for a plane wave incident along the ๐‘งห† direction. We define the error between the full isogeometric solution to those using MHs as ฮฆฬƒ โˆ’ ฮฆ ๐œ– ๐ฟ2 = , (3.10) kฮฆk where kยทk denotes the ๐ฟ2 norm, ฮฆ denotes the scattered far field obtained using purely the subdivision basis for physics and ฮฆฬƒ the solution obtained using MHs. As is evident from the first row in Table. 3.1, we observe convergence with increasing MHs. Numerical Example 2: Next, we study the reconstruction of both the geometry and the physics using MHs in tandem. In the second and third row in Table. 3.1, we present the geometric errors with respect to the number of MHs. The metrics we will use here and k๐‘† ๐‘€ โˆ’๐‘†๐บ k throughout the chapter are as follows: (a) Surface area error ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘€ ) = k๐‘†๐บ k , where ๐‘† ๐‘€ is the reconstructed surface area for ๐‘€ MHs and ๐‘†๐บ is the exact surface area. And, (b) ๐ป (ฮ“๐บ , ฮ“๐‘€ ), which denotes the maximum Hausdorff distance [43], between the exact geometry and the reconstructed geometry for ๐‘€ MHs. In addition, in Figs. 3.4a-3.4d both the reconstructed surface as well as the reconstructed and exact surfaces are prejected onto the ๐‘ฅ-๐‘ง, ๐‘ฆ-๐‘ง, and ๐‘ฅ-๐‘ฆ planes for direct comparisons, indicated as black meshed and gray solid regions, respectively. Finally, Fig. 3.4e demonstrates convergence as one increases the number of MHs for 26 ๐‘€ 500 1000 1500 2000 2562 Rel. ๐œ– ๐ฟ2 error ๐œ– 6.19E-4 3.50E-5 1.58E-5 9.55E-6 1.25E-13 ๐‘€ 10 50 100 200 ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ ๐œ– 2.46E-1 1.89E-2 4.13E-3 2.32E-3 ๐‘€ 10 50 100 200 ๐ป (ฮ“๐บ , ฮ“๐‘› ) ๐œ– 3.10E-1 4.42E-2 3.97E-2 3.30E-17 Table 3.1: Relative error metric with respect to number of manifold harmonics M. both geometry and physics. Each of the plots corresponds to fixed MHs for geometry with increasing MHs for physics. The error is measured against those obtained via an isogeometric solve. As is evident, it is possible to represent both the geometry and scattered field using significantly fewer MHs. It should also be noted that the errors plateau soon after a threshold error is reached in geometry reconstruction. Parenthetically, we note that the number of MHs used for shape reconstruction is typically not the same as those used to construct the physics on the surface. A good illustrative example is a sphere. To reconstruct the sphere, one needs only the lowest order, but reconstruction of the fields on the surface requires far more MHs proportional to the frequency of these fields. As a result, given ๐‘๐‘ฃ MHs and a desired fidelity, one would have different number of MHs to represent the surface and physics on the surface. 3.2.0.3 Analysis of Vector Fields using Manifold Harmonics Furthermore, we aim to construct a spectral representation of vector field quantities on the surface using MHB. The MHB forms the building block for a complete system of eigenfunctions of the vector Laplaceโ€“Beltrami operator (or Hodge-Laplace operator) ฮ” ยฎฮ“ = โˆ‡ฮ“ divฮ“ โˆ’ curlฮ“ curlฮ“ . Indeed, the system {โˆ‡ฮ“ ๐ป๐‘š , curlฮ“ ๐ป๐‘š } forms a system of orthogonal nontrival eigenvectors for ฮ” ยฎ ฮ“ with the same eigenvalues ๐œ†๐‘š , โˆ’ฮ” ยฎ ฮ“ โˆ‡ฮ“ ๐ป๐‘› = ๐œ†๐‘› โˆ‡ฮ“ ๐ป๐‘› , (3.11) โˆ’ฮ”ยฎ ฮ“ curlฮ“ ๐ป๐‘› = ๐œ†๐‘› curlฮ“ ๐ป๐‘› . (3.12) 27 (a) 10 Geometric MHs. (b) 50 Geometric MHs. (c) 100 Geometric MHs. (d) 200 Geometric MHs. (e) ๐œ– ๐ฟ2 relative error. Figure 3.4: Convergence in ๐œ– ๐ฟ2 relative error of reconstructed far field as we increase MHs used to represent the surface and fields on the surface. 28 (a) | Jฬƒ500 | (b) | Jฬƒ1000 | (c) | Jฬƒ2000 | (d) | Jฬƒ12130 | (e) |J| Figure 3.5: Reconstruction of the target current J obtained with an increasing number of MHs for the jet airliner. Given the representation of each of the eigenfunction (3.4), it follows that currents defined on the manifold can be written in terms of these eigenfunctions. Specifically, for the Helmholtz decomposition (2.5) prescribed earlier, we seek the gradient and rotation as follows: ๐‘๐‘ฃ๐‘™ ร• โˆ‡ฮ“ ๐ป๐‘š (r) โ‰ˆ โˆ‡ฮ“ ๐ปe๐‘š๐‘™ (r) = ๐‘™ โ„Ž ๐‘š,๐‘› โˆ‡ฮ“ ๐œ‰๐‘–๐‘™ (r), (3.13a) ๐‘›=1 ๐‘๐‘ฃ๐‘™ ร• e๐‘š๐‘™ (r) = curlฮ“ ๐ป๐‘š (r) โ‰ˆ curlฮ“ ๐ป ๐‘™ โ„Ž ๐‘š,๐‘› nฬ‚ ร— โˆ‡ฮ“ ๐œ‰๐‘›๐‘™ (r). (3.13b) ๐‘›=1 Using these expressions, the currents, J(r), may alternatively be written in in terms of this basis as (a) Figure 3.6: Example of a Helmholtz decomposition of a vector field on a double torus into curl-free, div-free and harmonic field. 29 (a) | Jฬƒ200 | (b) | Jฬƒ1000 | (c) | Jฬƒ2000 | (d) | Jฬƒ5122 | (e) |J| Figure 3.7: Reconstruction of the target current J obtained with an increasing number of MHs for the bumpy cube. J๐‘™๐‘€ (r) = Jฬƒ1,๐‘€ ๐‘™ ๐‘™ (r) + Jฬƒ2,๐‘€ ๐‘™ (r) + Jฬƒ3,๐‘€ ๐‘™ (r) + Jฬƒ4,๐‘€ (r) ๐‘€ h ร• i (3.14a) ๐‘™ ๐‘™ ๐‘™ ๐‘™ = ๐‘ฃ๐‘š Jฬƒ1,๐‘š (r) + ๐‘ค๐‘š Jฬƒ2,๐‘š (r) +๐Ž ยฏ ๐‘€ (r), ๐‘š=1 โˆ‡ฮ“ ๐ปe๐‘š๐‘™ (r) ๐‘™ Jฬƒ1,๐‘š (r) = , ๐‘™ p ๐œ†๐‘š curlฮ“ ๐ป e๐‘š๐‘™ (r) ๐‘™ Jฬƒ2,๐‘š (r) = (3.14b) ๐‘™ p ๐œ†๐‘š ร•๐‘” h i ๐Žยฏ ๐‘™๐‘€ (r) = ๐‘™ ๐‘Ž3,๐‘š ๐‘™ Jฬƒ3,๐‘š (r) + ๐‘Ž4,๐‘š๐‘™ ๐‘™ Jฬƒ4,๐‘š (r) . ๐‘š=1 such that {Jฬƒ1,๐‘š๐‘™ ๐‘™ , Jฬƒ2,๐‘š } is an othronormal basis [30]. And, ๐Ž ยฏ ๐‘™๐‘€ is the harmonic component represented using the MHB following (2.10). Numerical Example 1: In this example we consider representation of currents on two different simply connected objects: a bumpy cube and a jet airliner. Our goal is to examine the convergence of the representation of the current to a bandwidth of ๐‘€ harmonics. In both instances, we reconstruct a surface current generated by a 1 GHz plane wave incident in the โˆ’zฬ‚, respectively. In Fig. 3.3, we visualize the manifold harmonic representation of the current for a bumpy cube, with 5124 DoFs, and in Fig. 3.2 for a jet airliner, with 12132 DoFs. As can be seen in both figures, the first J๐‘š functions capture the coarse features of the current and the next, high frequency ones, correspond to the details. Table. 3.2 demonstrates the precision of the inverse MHT (3.14) w.r.t original current J๐‘ as we increase the number of 1/2 MHs. Our metric for validation is the reconstruction error ๐œ– = J๐‘ (x๐‘– ) โˆ’ Jฬƒ(x๐‘– ) [๐ต] . Note, 30 ๐‘€ 200 1000 2000 5122 Bumpy Cube ๐œ– 9.87E-4 4.02E-4 4.16E-5 3.65E-17 ๐‘€ 500 1000 2000 12130 Jet airliner ๐œ– 4.21E-4 9.93E-5 4.14E-5 9.03E-17 Table 3.2: Relative ๐œ– error, with respect to number of manifold harmonics M, in the reconstructed surface currents density. J๐‘ (x) is the current on the surface as approximated by the Loop-subdivision basis set. In both candidate objects, we find that as expected, ๐œ– decreases as the number of MH ๐‘€ increases, eventually approaching machine precision. Numerical Example 2: In what follows, we illustrate the prowess of MHB for vector field representation on a multiply connected surface. To wit, we consider representation of vector fields on the double torus (Fig. 3.6a) with 6140 DoF. Our goal is to examine the convergence of the representation of the vector field to a bandwidth of ๐‘€ harmonics. In this instance, we reconstruct the tangential projection of the vector field, L ร— (x โˆ’ x0) B(x) = , (3.15) |x โˆ’ x0 | 3 on the surface; note, this field is centered at the origin and oriented in the direction (0.1, 0.2, 2.1). In general, this field will have non-zero projections on to all three components in the Hodge decomposition. M k ๐œ–1 k 2 k ๐œ–2 k 2 k ๐œ–3 k 2 k ๐œ–4 k 2 500 3.37E-3 2.97E-3 2.68E-2 1.30E-1 1000 5.41E-5 5.62E-5 5.67E-3 1.32E-2 2000 1.69E-5 5.56E-7 2.04E-3 4.02E-3 6138 6.94E-27 1.25E-26 5.54E-24 1.31E-23 Table 3.3: ๐ฟ2 -error of the reconstructed components of a vector field on the double torus (Fig. 3.6a), using M MHs. 31 Table. 3.3 demonstrates the precision of the inverse MHT (3.14) w.r.t original vector field J๐‘ on the surface, as approximated by the Loop-subdivision basis set, as we increase the number of MHs. Our metric for validation is the ๐ฟ2 error metric as follows: ๐œ– ๐œˆ = q k Jฬƒ๐œˆ,๐‘€ โˆ’ J๐œˆ,๐‘ k 2 /k J๐œˆ,๐‘ k 2 , where ๐œˆ โˆˆ {1, 2, 3, 4}. We find that as expected, ๐œ– decreases as the number of MH ๐‘€ increases. Where the ๐ฟ2 norm of each reference component is as follows: k J1,๐‘ k 1 = 6.02E+2, k J2,๐‘ k 1 = 5.67E+2, k J3,๐‘ k 2 = 2.18E-1, and k J4,๐‘ k 2 = 2.67E-5. In Fig. 3.6a, we visualize the Helmholtz decomposition of the Loop subdivision representation of the current. 3.3 Localized Manifold Harmonics 3.3.1 Formulation In this section, we follow a framework for spectral analysis that is designed to be at the same time local and compatible with the existing manifold harmonic spectral constructions [113]. In practice, this approach boils down to the computation of the eigenfunctions of a new operator, which is realized as a simple update to the Laplace-beltrami operator (3.2) thus fully retaining the computational efficiency and theoretical guarantees of the resulting optimization process. As alluded to earlier, the principal challenge with MHB is the number of basis functions necessary to capture fine features, be it geometry or physics. A potential remedy is the use of L-MHB [113]. L-MHB is a local and compatible basis with the existing MHB. In practice, computing L-MHB boils down to the computation of the eigenfunctions of a new operator, which is realized as a simple update to the LBO (3.2) thus fully retaining the computational efficiency and theoretical guarantees of the resulting optimization process. Considering a region ฮ“๐›ผ โŠ‚ ฮ“, we seek a set of L-MHB {๐œ1 , . . . , ๐œ ๐‘€ 0 } as the solution to the following optimization problem: ๐‘€0 ร• ๐‘š๐‘–๐‘› โ„ฐ๐‘† (๐œ ๐‘— ) + ๐œ‡๐‘… โ„ฐ ๐‘… (๐œ ๐‘— ), (3.16a) ๐œ1 ,...,๐œ ๐‘€0 ๐‘—=1 32 s.t. h๐œ ๐‘– , ๐œ ๐‘— i๐ฟ2 (ฮ“) = ๐›ฟ ๐‘–๐‘— ๐‘–, ๐‘— = 1, . . . , ๐‘€ 0 , (3.16b) h๐œ ๐‘– , ๐ป ๐‘— i๐ฟ2 (ฮ“) = 0, ๐‘– = 1, . . . , ๐‘€ 0; ๐‘— = 1, . . . , ๐‘€, (3.16c) where the first term, โ„ฐ ๐‘  (ยท), is as defined earlier and โˆซ 2 โ„ฐ๐‘… ( ๐‘“ ) = ๐‘“ (r)(1 โˆ’ ๐‘ข(r) ๐‘‘r,  (3.17) ฮ“ is a quadratic penalty promoting the localization of the basis functions on the given region ฮ“๐›ผ โŠ† ฮ“, via the local penalty term ๐œ‡๐‘… . Here, we choose ๐‘ข(r) such that ๐‘ข(r) = 1 for r โˆˆ ฮ“๐›ผ and ๐‘ข(r) = 0 otherwise. As we will see in what follows, (3.16) allows constructing an incremental set of functions that are smooth, orthogonal to a given set of MHB, and localized to ฮ“๐›ผ . 3.3.2 Implementation Details Next, to obtain the solution to (3.16) we use a penalty method with the following penalty term ๐œ‡โŠฅ [113]. With this in mind, let the limit surface ฮ“ be parametrized by the ๐‘™-th refined control mesh ๐‘‡ ๐‘™ and [๐œ ๐‘™ ] โˆˆ R๐‘›ร—๐‘€ be a matrix containing the ๐‘€ 0 L-MHB in its columns. 0 Equation (3.16) can then be written in discrete form as, [๐‘„ ๐‘™ ][๐œ ๐‘™ ] = [๐ต ๐‘™ ][๐œ ๐‘™ ][ฮฃ๐‘™ ], (3.18a) [๐‘„ ๐‘™ ] = [๐ด ๐‘™ ] + ๐œ‡๐‘… [๐ต ๐‘™ ]diag(v) + ๐œ‡โŠฅ [๐ต ๐‘™ ][๐ป ๐‘™ ][๐ป ๐‘™ ]๐‘‡ [๐ต ๐‘™ ], (3.18b) where [๐‘„ ๐‘™ ] is symmetric, positive semi-definite, and v denotes the discrete version of v(r) = (1 โˆ’ ๐‘ข(r))2 that serves as a diagonal update to enforce locality in (3.17). By casting the optimization problem as a generalized eigenvalue problem, we obtain a natural ordering of the basis functions; for ๐œ‡๐‘… > 0 and ๐œ‡โŠฅ > 0, we obtain a new set of ๐‘€ 0 functions { ๐œหœ 1 , . . . , ๐œหœ ๐‘€ 0 } localized to a given region ฮ“๐›ผ โŠ‚ ฮ“. These functions 33 (a) Figure 3.8: Localized manifold harmonics for the bumpy cube. We show the first few localized manifold harmonics. effectively extend the MHB, in the sense that the new set { ๐ปหœ 1 , . . . , ๐ปหœ ๐‘€ , ๐œหœ 1 , . . . , ๐œหœ ๐‘€ 0 } forms an orthogonal basis for a ๐‘€ + ๐‘€ 0-dimensional subspace of ๐ฟ2 (ฮ“). Note, we are guilty of a slight abuse of notation in the interest of brevity; we have directly used the discrete or quantities denoted using โ€˜tilde.โ€™ The localized harmonics bring additional (higher frequency) information to the global basis formed by the first ๐‘€ MHB. Compared to only using only the MHB, the new representation provides a more compressed model since fewer localized harmonics are needed to capture the high-frequency content within ฮ“๐›ผ , than the number of harmonics that would be needed when considering only the MHB. Lastly, Fig. 3.8a provides an illustration of the first few LMH on the a bumpy cube. 3.4 Hierhical Manifold Harmonic Basis on Loop subdivised surfaces 3.4.1 Formulation In this section, we demonstrate the construction of H-MHB via a series of control meshes and apply it to the analysis of electromagnetic scattering. A salient feature of subdivision is that all control meshes point to the same limit surface. It follows, therefore, that it is possible to relate MHB defined on a mesh ๐‘‡ ๐‘™ to mesh ๐‘‡ ๐‘™+1 . In other words, we can define a hierarchical basis set that is defined band-by-band over a series of Loop subdivision control meshes. This approach fits well with the common numerical computation of eigenvectors and eigenvalues, which is done per band, leading to an efficient basis computation. Furthermore, this allows reduction in complexity when using this basis for electromagnetic analysis, as will be demonstrated below. 34 (a) (b) Figure 3.9: A select few low manifold harmonics and the corresponding prolonged manifold harmonics from a level below on a tetrahedron. Fortunately, the mapping between different control meshes are well defined within the rubric of subdivision; specifically, prolongation and restriction operators allow us to map data defined on one mesh onto another mesh. We use this framework highlighted in the last chapter to construct prolonged manifold harmonics [144]. Specifically, the prolonged MHB at level ๐‘™ + 1 from level ๐‘™ are given by: [๐ปห† ๐‘™+1 ] = [๐‘ƒ ๐‘™ ][๐ป ๐‘™ ], (3.19) and, the coarsening operator is defined as [๐ป ๐‘™ ] = ๐‘… ๐‘™ [๐ปห† ๐‘™+1 ], (3.20) where, the prolongation operator is the Loop subdivision matrix (2.2) and restriction operator is defined as [๐‘… ๐‘™ ] = ([๐‘† ๐‘™ ]๐‘‡ [๐‘† ๐‘™ ])โˆ’1 [๐‘†]๐‘‡ , see (2.3). 3.4.2 Implementation Details We now introduce the notations for hierarchically partitioning the MHB, as similarly proposed in [144]. The computation of the H-MHB basis requires splitting the eigenspace into ๐‘“ + 1 non-overlapping bands defined over their respective ๐‘™ ๐‘ก โ„Ž refined mesh ๐‘‡ ๐‘™ , given ร๐‘“ the required band widths {๐‘˜ ๐‘™ โ‰ฅ 0, ๐‘™ 0, . . . , ๐‘“ }, such that ๐‘™=0 ๐‘˜ ๐‘™ = ๐‘˜หœ ๐‘“ . The hierarchical basis is given by [๐ปหœ 0 ] = [๐ป 0 ], [ฮ›ฬƒ0 ] = [ฮ›0 ], and [๐ปหœ ๐‘™+1 ] = [๐‘† ๐‘™ ๐ปหœ ๐‘™ , ๐ป ๐‘™+1 ], and [ฮ›ฬƒ๐‘™+1 ] = [ฮ›ฬƒ๐‘™ , ฮ›๐‘™+1 ] where [๐ป 0 ], [ฮ›0 ] are the first ๐‘˜ 0 MHB and respective eigenvalues and [๐ปหœ ๐‘™+1 ], [ฮ›ฬƒ๐‘™+1 ] are the first ๐‘˜ ๐‘™+1 H-MHB and eigenvalues in the band starting after the largest eigenvalue 35 Algorithm 1: Compute H-MHB. input : A series of control meshes ๐‘‡ ๐‘™ , ๐‘™ โˆˆ {0, . . . , ๐‘“ } ร๐‘“ Total dimension of MHB ๐‘˜หœ ๐‘“ = ๐‘™=0 ๐‘˜ ๐‘™ h - size of MHB band to compute Control error ๐œ– output : Set of H-MHB [๐ป] หœ and eigenvalues [ฮ›ฬƒ] 1 // Compute [๐ป 0 ] and its eigenvalues [ฮ›0 ] 2 ๐œŽ = max(ฮ›0 ) + 0.4(max(ฮ›0 ) โˆ’ min(ฮ›0 )) 3 // Store [๐ป 0 ] 4 ๐ป โˆ— โ† [๐ป 0 ] 5 for ๐‘™ := 1 to ๐‘“ do Mesh loop 6 ๐‘˜โ†0 7 while ๐‘˜ < ๐‘˜ ๐‘™ do MH loop 8 Compute โ„Ž MHB ๐ปยฏ and eigenvalues ฮ›ฬ„ near ๐œŽ 9 Remove ๐ปยฏ for which h๐ปห† โˆ— , ๐ปi ยฏ ฮ“>๐œ– 10 //Concatenate remaining ๐‘˜ 0 ๐ปยฏ and ฮ›ฬ„ 11 หœ = [๐ป โˆ— , ๐ป] [๐ป] ยฏ 12 [ฮ›ฬƒ] = [ฮ›โˆ— , ฮ›ฬ„] 13 ๐œŽ = max(ฮ›ฬ„) + 0.4(max(ฮ›ฬ„) โˆ’ min(ฮ›ฬ„)) ๐‘˜ = ๐‘˜ + ๐‘˜ 0 14 end 15 //Store ๐ปยฏ 16 ๐ป โˆ— โ† ๐ปยฏ 17 end of [ฮ›ฬƒ๐‘™ ] and the operator [ยท, ยท] denotes either column concatenation, or diagonal matrix concatenation according to the context. For each band we employ an eigenvalue solver, that computes eigenvalues and eigen- vectors near a given ๐œŽ that is close to an eigenvalue [144, 134]. We iteratively compute the bands as shown in Alg. 1; this tactic is a good fit with the common numerical computation of eigenvectors and eigenvalues, which is done per band [153], leading to an efficient basis computation. We note that, as the eigenvalues are not exactly linear, we overcompensate in the computation of the bands, by computing โ„Ž eigenpairs more than what is required. Then, we leverage the approximate orthogonality property, to remove MHs that are already well represented in the basis ๐ปหœ ๐‘™ , where we discard eigenvectors with a maximal projection norm larger than ๐œ–. In all our experiments we choose โ„Ž = 50 and ๐œ– = 0.9. In this scheme, it 36 (a) (b) (c) Figure 3.10: Two successive Loop subdivison refinements, from coarse to refined, of control mesh (a) 540 vertices, (b) 2050 vertices and (c) 8194 vertices. difficult to guarantee completeness of the basis set, i.e., none of the information is โ€œlost", as otherwise, this reduces the basisโ€™ representation power; a formal guaranteed of its completeness doesnโ€™t exist [144]. The method in which we split the eigenspectrum is heuristic, as such, issues like repeating eigenvalues may be problematic. However, in our experiments, we have not experienced problems due to this limitation. Lastly, Fig. 3.9 shows the visual similarity between the MHB for [๐ป] and the prolonged MHB [๐ป] ห† for a few low eigenfunctions. 3.5 Conclusion In this chapter, we have explored the benefits of MHB defined on Loop subdivision surfaces as applied to both the manifold and physics on the manifold. In particular we have exploited the compression and multi-resolution framework that this basis provides to address every facet of the analysis process; from representing and compressing/filtering geometry to scalar and vector fields on Loop subdivision surfaces; including representing currents on multiply connected surfaces. In addition, we review a framework for local spectral analysis and processing that naturally extends the MHB spectrum, allowing one to perform operations which are localized to a given region of interest on the surface. Lastly, we focused on explotiing the multi-resolution nature of both the MHB subdivision surfaces to implement the H-MHB as a means of using a composition of data from a hierarchy of refined geometries to compute results from the finest geometry. The application of these bases is explored in later chapters including problems in inverse design, acoustic 37 scattering, and electromagnetic scattering. Moreover, MHB opens the door to far-reaching applications like machine learning, target recognition, data fusion for inverse problems, shape and topology optimization. This chapter, ยฉ2022 JASA, in part, is reprinted, with permission, from A. M. A. Alsnayyan and B. Shanker, "Laplace-Beltrami based multi-resolution shape reconstruction on subdivision surfaces," The Journal of the Acoustical Society of America on, March 2022, with slight modifications to fit in this dissertation. This chapter, ยฉ2022 IEEE, in part, is reprinted, with permission from A. M. A. Al- snayyan and B. Shanker, "Manifold Harmonics and their application to computational electromagnetics," IEEE Journal on Multiscale and Multiphysics Computational Techniques on, August 2022, with slight modifications to fit in this dissertation. 38 CHAPTER 4 ACOUSTIC THEORY AND INTEGRAL EQUATION FORMULATIONS 4.1 Introduction Surface integral equation (SIE) posit excellent benefits such as (a) the Greenโ€™s function embeds the necessary radiation condition; (b) reduction of a volumetric problem into a surface embedded in a volume, and therefore, the required degrees of freedom; and (c) readily handle arbitrary geometries. As such, this method has become commonplace for computing scattering from piecewise homogeneous scatterers [89, 156, 91]. The traditional downsides of these methods are (a) non-uniqueness at irregular frequencies [36], (b) computational complexity that scales as ๐’ช(๐‘๐‘2 ), where ๐‘๐‘ denotes the number of degrees of freedom used in the analysis [11], and (c) the integration of singular functions. By computational complexity, we refer to the scaling of memory required as well as execution time of each matrix vector product required in any iterative solver. We note that both downsides have been overcome [68]. Methods have been developed to obtain a unique solution at all frequencies and to reduce the cost to ๐’ช(๐‘๐‘ log ๐‘๐‘ )[65, 74]. Given these advances, SIE-based computational methods are both reliable and robust. However, despite all appearances, a number of interesting challenges remain and these are not specific to acoustics. Indeed, as will be evident in this paragraph, some of these comments apply to all SIEs. A traditional SIE solution proceeds as follows: the surface of a candidate object is represented using a collection of discrete elements, the physical quantity of interest is then represented using basis functions described on these discrete elements, and a variational technique is used to reduce the continuous equations to a set of discrete equations. Typically, the object is represented using a collection of triangles, and low order ansatz basis functions are used to represent the physical quantity of interest on the scatterer. The resulting discrete set inherits (and depending on the operator, exacerbates) errors due to geometric approximation. It follows, then, that these challenges can be alleviated by 39 resorting to higher order geometric representation and the corresponding basis functions [126]. One would imagine that resorting to higher order geometric representation would result in fewer degrees of freedom, simply due to the fact that one does not need as many piece-wise flat triangles to represent the object. But one should note that in most higher order schemes the discretization is still ๐ถ 0 and functions are higher order within one triangle. Our approach to higher order analysis is different. In the past, we developed a method for discretizing SIEs that was highly flexible and adaptive called the Generalized Method of Moments (GMM) [116, 51]. GMM is effective and efficient, and has been used for different problems in both acoustics and electromagnetics. It permits refinement in the order of geometric representation and the order of basis functions used for physics, and spatial refinement. In GMM, geometry is constructed using a set of overlapping patches, with physical basis functions defined on each patch. Since each patch is constructed using a least-squares approach, a persistent drawback of GMM is that these parametrizations are not watertight. Robust higher order watertight surface representation is indeed vogue in the computer graphics literature. Here, the surface representation is watertight and approximately ๐ถ 2 almost everywhere (provided one uses subdivision [106, 149]). The higher order continuity of the surface representation offers a number of tantalizing opportunities from both mathematical and engineering perspectives such as analysis of convergence, development of Calderรณn type preconditioners, Laplace-Beltrami based compression, higher order basis for physics that has the same continuity properties as geometric basis sets, shape and topology optimization methods that rely on compressed systems, and so on. In this chapter, we set the stage for this body of work by asking a more fundamental questionโ€“what should we have in our arsenal so as to realize computationally efficient higher order (in geometry and physics) analysis tools. Our approach to meeting the goals alluded to earlier is to use isogeometric analysis. Here, the basis functions used to represent the physics are the same as basis functions 40 used to represent the geometry. As a result, all the geometric modeling features that one desires, such as adaptivity, refinement, etc., directly carry over to modeling the physics. In a manner of speaking, in the long term such an analysis tool would enable seamless integration of computer aided design (CAD) and computer aided engineering (CAE). The advent of isogeometric analysis has its genesis in early papers by Hughes [79], and has since been adapted for a number of different types of problems [49, 77] in structural mechanics and electromagnetics [98, 97]. While our earlier body of work on extending isogeometric analysis based on subdi- vision surfaces to electromagnetics [98, 67] demonstrated its viability, several challenges related to computational cost, both evaluation of matrix elements as well as matrix vector products, were left unaddressed. It is important to note that the Greenโ€™s kernels used for electromagnetics were all due to single layer potentials, unlike those for acoustics [146]. The latter makes evaluation significantly more challenging. To help visualize challenges in evaluating the necessary inner products, it should be noted that unlike Lagrangian discretization, the domain of a subdivision basis function is significantly larger; it includes all triangles that share a node (or control point). The basis function used for geometry is a product of box-splines and therefore, fourth order. Since we use the same basis function for physics as well, it implies that one needs to evaluate products of higher order polynomials convolved with a Greenโ€™s function on higher order surfaces, for all pairwise interactions as is dictated by an SIE formulation. Assuming that ๐‘๐‘ is the number of degrees of freedom, and ๐‘๐‘ž is the number of quadrature points, the cost of evaluating the matrix elements scales as ๐’ช(๐‘๐‘ž2 ๐‘๐‘2 ). We note that ๐‘๐‘ž can be very high as both the surface and the basis functions are fourth order polynomials. Our objective is to develop an isogeometric analysis technique for acoustic scattering that is computationally efficient. Specifically, we seek to reduce the cost of evaluation of both the matrix elements as well as a matrix vector product. In order to do so, we exploit a version of the Multi-level Fast Multipole Method (MLFMA) that is stable at both low and high frequencies [112, 80] and has been used extensively 41 for wideband electromagnetic analysis. To this end, the principal contributions of are as follows: we present โ€ข isogeometric analysis (IGA) on smoothly represented structures of arbitrary topology, โ€ข a study of convergence properties of IGA based SIE solvers, โ€ข methods to linearize both the cost of evaluation of the necessary matrices and matrix-vector products, โ€ข a demonstration of convergence of IGA-SIE solver for conical geometries as well as complex targets with multiscale features, โ€ข results validating the scaling and accuracy of wideband MLFMA scheme, and The rest of the chapter is organized as follows; in Section 4.2.1, we present the formulation, describe subdivision basis for geometry, and discuss the resulting linear system. In Section 4.3, we present a methodology to efficiently evaluate all the requisite inner products. This is followed by asymptotic complexity estimates in Section 4.3.3, and a slew of numerical results in Section 6.4. Finally, in Section 6.5, we summarize our contributions. 4.2 Acoustics: Helmholtzโ€™s Equation and Integral Equation Formula- tions 4.2.1 Problem Formulation Consider a scatter, in a homogeneous medium, with a hard boundary ฮฉ โˆˆ R3 and a uniquely defined normal ๐‘›(r),ห† where r โˆˆ ฮฉ. It is assumed that the rest of space is homogeneous, and that the object is at rest. Consider a velocity field incident on this scatterer denoted by v๐‘– (r). This generates a scattered velocity field given by v๐‘  (r) and we define the total velocity field as v๐‘ก (r) = v๐‘– (r) + v๐‘  (r). These fields can be represented by equivalent potentials such that v๐œ‰ (r) = โˆ‡ฮฆ๐œ‰ (r), for ๐œ‰ โˆˆ {๐‘ก, ๐‘ , ๐‘–}. A unique solution in ฮฆ๐‘ก (r) is 42 guaranteed through the Burton-Miller formulation, written in terms of these potentials[36]: โ„’ ๐ต๐‘€ [ฮฆ, ฮฉ] (r) = ๐‘‰ ๐‘– (r),   ฮฆ โ„’ ๐ต๐‘€ [ฮฆ, ฮฉ] (r)  (1 โˆ’ ๐›ผ) โˆ’ ๐’Ÿ[ฮฆ, ฮฉ](r) โˆ’ ๐›ผ๐›ฝ๐’ฉ[ฮฆ, ฮฉ](r), (4.1) 2 ๐‘‰ ๐‘– (r)  (1 โˆ’ ๐›ผ)ฮฆ๐‘– (r) + ๐›ผ๐›ฝ ๐‘›ห† ยท โˆ‡ฮฆ๐‘– (r). where 0 โ‰ค ๐›ผ โ‰ค 1, the constant weighting factor and ๐›ฝ, the coupling factor, make the solution unique at all frequencies; the optimal choice of alpha and beta has been the subject of many studies[12, 110, 169]. In the remainder of the chapter we follow the standard Burton-Miller formulation[36]. Note, setting ๐›ผ = 0 or ๐›ผ = 1 introduces a non-trivial null space at frequencies that correspond to the interior resonance of the structure. The reader is referred to [36, 104] for a theoretical explanation. Insight into the null spaces of these operators for spherical objects can be obtained using methods presented in [99]. The operators ๐’Ÿ and ๐’ฉ are defined as โˆซ ๐’Ÿ[ฮฆ, ฮฉ](r)  ฮฆ(r0)nฬ‚(r0) ยท โˆ‡0 ๐บ(r, r0)๐‘‘r0 , (4.2a) ฮฉ โˆซ  ๐’ฉ[ฮฆ, ฮฉ](r)  nฬ‚(r) ยท โˆ‡ ฮฆ(r0)nฬ‚(r0) ยท โˆ‡0 ๐บ(r, r0)๐‘‘r0 . (4.2b) ฮฉ We define the kernel as the free-space Helmholtz kernel, ๐บ(r, r0) = ๐‘’ โˆ’๐‘– ๐‘˜0 |rโˆ’r | /4๐œ‹|r โˆ’ r0 |, ๐‘˜0 is 0 the wavenumber, ๐›ฝ = ๐‘–/๐‘˜ 0 . An ๐‘’ ๐‘–๐œ”๐‘ก dependence is assumed and suppressed. The solution to these equations is typically effected in a discrete setting as follows; one discretizes ฮฉ ร ๐‘๐‘ using a collection of patches such that ฮฉ = ๐‘™ ฮฉ๐‘™ , where ฮฉ๐‘™ denotes a patch. Then each patch ฮฉ๐‘™ supports a unique, locally indexed, set of basis functions ๐œ“ ๐‘› (r), ๐‘› โˆˆ {1, ยท ยท ยท , ๐‘๐‘ฃ }. In other words, each basis function, ๐œ“ ๐‘› , is defined over a unique collection of patches ฮ“๐‘› = {ฮฉ๐›ผ , ยท ยท ยท , ฮฉ๐‘™ , ยท ยท ยท , ฮฉ๐›ฝ }. The exact solution of (4.1) is generally not available. The numerical solution of the integral equations is effected in a discrete setting using an isogeometric method based on Loop subdivision. 43 4.2.2 Isogeometric basis sets and discrete systems In chapter 2, we have defined effective basis functions that act upon control vertices so as to produce a limit surface. In what follows, we postulate that one can define quantities at these control nodes that represent the physics on the limit surface. To this end, one can write the potential over any given patch as ร•๐‘๐‘ฃ ฮฆ(r(๐‘ข, ๐‘ฃ)) = ๐‘Ž ๐‘– ๐œ“ ๐‘– (r(๐‘ข, ๐‘ฃ)). (4.3) ๐‘–=1 Where now we define ๐‘Ž ๐‘– as weights assigned to the locally indexed ๐‘– ๐‘ก โ„Ž vertex; ๐‘๐‘ฃ , ๐œ“ ๐‘– (r(๐‘ข, ๐‘ฃ)), and (๐‘ข, ๐‘ฃ) retain the same definition as those prescribed above. Employing Galerkin testing with these prescribed basis sets, one obtains ๐’ตโ„ = ๐’ฑ, (4.4a) โˆซ ๐‘ ๐‘š๐‘› = ๐œ“ ๐‘š (r)โ„’ ๐ต๐‘€ ๐œ“ ๐‘› , ฮ“๐‘› (r)๐‘‘r,   (4.4b) ฮ“๐‘š where โ„ = [๐‘Ž 1 , ๐‘Ž2 , ยท ยท ยท , ๐‘Ž ๐‘ ]๐‘‡ , ๐’ฑ = [๐‘‰1 , ๐‘‰2 , ยท ยท ยท , ๐‘‰๐‘ ]๐‘‡ and โˆซ ๐‘‰๐‘š = ๐œ“ ๐‘š (r)๐‘‰ ๐‘– (r)๐‘‘r. (4.4c) ฮ“๐‘š While simply stated, evaluation of these integrals is challenging from mulitple perspectives: the hypersingular nature of these integrals needs to be carefully taken into account when the support ฮ“๐‘š โˆฉ ฮ“๐‘› โˆ‰ โˆ…. Likewise, higher order quadrature poses significant bottlenecks. These issues are addressed in the next section. 4.3 Wideband MLFMA for higher order isogeometric bases As alluded to earlier, the principal challenge is the evaluation of matrix element ๐‘ ๐‘š๐‘› . The order of the basis function (for both geometry and physics) is fourth order. Naively, the order of red integration rule for evaluating the inner products scale as ๐‘๐‘ž โˆ 16 per triangle, giving an effective scaling as ๐‘๐‘ž2 โˆ 256๐‘ฮ”2 where ๐‘ฮ” is the number of patches associated with the control vertices ๐‘š and ๐‘›; note that using ๐‘ฮ” is an abuse of notation, as ๐‘ฮ” can differ for ๐œ“ ๐‘š and ๐œ“ ๐‘› . As is apparent, this cost quickly dominates the overall cost. 44 Note that if the two nodes are in proximity to each other, one needs additional structure to evaluate the hypersingular integrals, this interaction is illustrated in Fig. 4.1. The approach we espouse is to construct an adaptive quadrature scheme that can be accelerated via wideband MLFMA. To set the stage for the computation, assume that each triangle in the domain ฮ“๐‘› of the ๐‘›th basis function is partitioned recursively into 4๐‘™ subpatches where ๐›พ ๐‘™ โˆˆ N. The domain of each of these partitions is now denoted by ฮ“๐‘› for ๐›พ = 1, ยท ยท ยท , 4๐‘™ ๐‘ฮ” . It follows that any matrix element can be computed in terms of its partial contributions such that ร•ร•โˆซ ๐›พ ๐›พ ๐‘ ๐‘š๐‘› = ๐œ“ ๐‘š (r)โ„’ ๐ต๐‘€ ๐œ“ ๐‘› (r), ฮ“๐‘› ๐‘‘๐‘†  ๐œ ๐œ ๐›พ ฮ“๐‘š (4.5) ๐œ๐›พ ร•ร• = ๐‘ ๐‘š๐‘› ๐œ ๐›พ ๐›พ ๐›พ where โ„’ ๐ต๐‘€ denotes the evaluation of the operator over the domain ฮ“๐‘š , and the summations is over the number of source and observation subpatches. In our implementation we have fixed the depth of recursion ๐‘™. We note the following: โ€ข The size of each sub-patch is typically sub-wavelength; often less than a tenth of a wavelength. ๐œ ๐›พ โ€ข If the interaction is a self interaction, i.e., ๐‘š = ๐‘›, and ฮ“๐‘š = ฮ“๐‘› , one needs to use Hadamard finite part integration[156]. Next, we prescribe a methodology to ameliorate the cost of evaluating (4.5). 4.3.1 Amelioration of Matrix Vector Product Costs To reduce the cost above, we take recourse to a tree based algorithm that is robust at small length scales. Fast multipole methods have long been in vogue in both acoustics and electromagnetics [161], their use has been to reduce the cost of evaluating a matrix vector product. As has been shown in Chew et. al. [41]. the manner which this approach has been used is to first enclose the entire object in a cubical domain (henceforth, called the root node), and recursively subdivide until one reaches a prescribed size of the smallest box (henceforth, termed the leaf node). Then at each level one partitions interactions between 45 Figure 4.1: Interaction between two nodes wherein the black subpatch denotes the source, red patches are its nearfield neighbors and all other patches lie in its farfield region. (a) (b) Figure 4.2: Comparison of interactions between FMM and wideband MLFMA. (a) Limit on traditional leaf box size means that interactions between large patches must be computed directly. (b) For wideband MLFMA, interactions between large patches may be split between nearfield and farfield such that favorable scaling is maintained. all boxes at that level as being in each others near or far field via the following rules: (a) any two boxes are in each others near field provided they share a geometrical region in space (node/edge/face), (b) boxes are denoted as far field pairs if they are in each others far field and their parents are in each others near field. This partition of interactions ensures that near field interactions are done at the leaf level, and all other interactions are accounted for using far field operators at different levels in the tree; Fig. 4.2 illustrates this computational infrastructure. In effect, using the above computational procedure one computes a matrix vector product by partitioning ๐’ต = ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ + ๐’ต ๐‘“ ๐‘Ž๐‘Ÿ . The matrix ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ is explicitly computed and stored, and ๐’ต ๐‘“ ๐‘Ž๐‘Ÿ is never explicitly computed. The product ๐’ต ๐‘“ ๐‘Ž๐‘Ÿ โ„ is effected via 46 (a) (b) Figure 4.3: Illustration of traditional and wideband MLFMA trees. (a) Traditional: MLFMA leaf box size (and therefore tree height) is limited by the size of the largest patch. (b) Wideband MLFMA: leaf box size and tree height are not limited by patch size, and may be chosen for optimal MLFMA scaling. five tree operations; charge to multipole (C2M), multipole to multipole (M2M), multipole to local (M2L), local to local (L2L) and local to observer (L2O). These denote operations mapping radiation of sources onto a tree, traversal up and across the tree, and finally, mapping these fields onto observers. The nuances of these operations are discussed in several papers, but notationally consistent definitions can be found in [112, 80]. Implicit in the above discussion is that one indirectly classifies interactions between basis functions. The starting point of building a tree is to choose a leaf box size, typically chosen to be 0.2๐œ† where ๐œ† is the wavelength corresponding to the wavenumber ๐‘˜ 0 . All basis functions whose support, ฮ“๐‘š , overlaps significantly with the domain of a leaf box are assigned to that box. In a typical discretization, the domain of each unknown is approximately 0.1๐œ†. As a result, each unknown can be mapped into one leaf box. This is not the case for higher order discretizations wherein the domain of each basis function is significantly larger [51]. As a result, their domain would span multiple leaf boxes. One alternative is to make the leaf box size significantly larger, to say largest ฮ“๐‘š , for ๐‘š = 1, ยท ยท ยท , ๐‘๐‘ . An alternative would be to exploit the tree infrastructure to reduce the cost of computing each ๐‘ ๐‘š๐‘› โˆˆ ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ . It follows that given (4.5), the near interactions at the leaf node may instead be classified 47 ๐œ ๐›พ ๐œ๐›พ in terms of ฮ“๐‘š and ฮ“๐‘› , and one can instead compute ๐‘ ๐‘š๐‘› and this collection constitutes ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ . But we can do better than previous work[51], by choosing the leaf box size to ๐œ (and correspondingly, the size of ฮ“๐‘š ) to be arbitrarily small so as to achieve the optimal tradeoff between total computational cost for both near and far field computation; Fig. 4.2b illustrates these advantages. To enable choosing leaf box sizes smaller than 0.2๐œ†, we use wideband FMM [112] that is stable at small lengths. With no loss in generality, we denote boxes whose side lengths are smaller than 0.2๐œ† as Accelerated Cartesian Expansion (ACE) [143] boxes, and those at that size or larger, as FMM boxes. To this end, consider two leaf boxes that interact in the far field of each other. The domain of these leaf boxes are denoted using ๐œ” ๐‘  and ๐œ” ๐‘œ , respectively. Assume that ๐›พ ๐œ ฮ“๐‘› โˆฉ ๐œ” ๐‘  โ‰  โˆ…. Analogously, ฮ“๐‘š โˆฉ ๐œ” ๐‘œ โ‰  โˆ…. Let r๐‘๐‘œ and r๐‘๐‘  denote the centers of ๐œ” ๐‘œ and ๐œ” ๐‘  . Then using the notation introduced in [143, 112] โˆซ ๐›พ ๐›พ ๐œ“ ๐‘š (r)โ„’ ๐ต๐‘€ ๐œ“ ๐‘› (r), ฮ“๐‘› = โ„’ ๐ด โ—ฆ ๐’ฏ๐ด โ—ฆ โ„ณ ๐ด .  (4.6a) ๐œ ฮ“๐‘š Where the operators are defined as โ„ณ ๐ด  M(0) , ยท ยท ยท , M(๐‘ƒ) ,  โˆซ (r0 โˆ’ r๐‘๐‘  )(๐‘) (4.6b) (๐‘) ๐‘ M = (โˆ’1) ๐‘‘r0 ๐‘›ห† 0 ยท โˆ‡0 ๐œ“ ๐‘› (r0), ฮ“๐‘› ๐›พ ๐‘! ๐’ฏ๐ด โ—ฆ โ„ณ ๐ด  L(0) , ยท ยท ยท , L(๐‘ƒ)  ร• ๐‘ƒ (4.6c) L (๐‘) = โˆ‡(๐‘›) ๐บ(r๐‘๐‘œ , r๐‘๐‘  ) ยท (๐‘› โˆ’ ๐‘) ยท M(๐‘›โˆ’๐‘) , ๐‘›=๐‘ and ร• ๐‘ƒ (1 โˆ’ ๐›ผ) + ๐›ผ๐›ฝ ๐‘›ห† ยท โˆ‡ (r โˆ’ r๐‘๐‘œ )(๐‘) ยท ๐‘ ยท L(๐‘) .   โ„’ ๐ด โ—ฆ {ยท} = โˆ’ (4.6d) ๐‘=0 Eqs. (4.6b)-(4.6c) define C2M, M2L and L2O operations. In the above expressions, A(๐‘›) and r(๐‘›) (whose component can be expressed in compressed form as r(๐‘›) (๐‘›1 , ๐‘›2 , ๐‘›3 ) = ๐‘ฅ ๐‘›1 ๐‘ฆ ๐‘›2 ๐‘ง ๐‘›3 ) denote an ๐‘› ๐‘ก โ„Ž rank tensor and C(๐‘›โˆ’๐‘š) = A(๐‘›) ยท ๐‘š ยท B(๐‘š) denotes an ๐‘š-fold contraction between tensors of rank ๐‘› and ๐‘š. It is apparent that (4.6c) depends only on centers of leafs boxes 48 and not on either source or observation sub-triangles. As a result, using (4.6c) for all source subtriangles in ๐œ” ๐‘  and observation subtriangles in ๐œ” ๐‘œ enables cost savings. While these equations describe a 1-level setting, it can be extended to a multi-level wideband setting by nesting these operations[112]; note, details of these operators, proofs of convergence, cost complexity, etc., can be obtained from [143, 112]. 4.3.2 Implementation details In what follows, we lay out the algorithmic structure necessary to implement the formulation that has been outlined. Specifically, we will assume that one has an available geometric description comprising of control points, and a connectivity map regarding the neighborhood of each control point. In addition, we assume that one has the available necessary tools to obtain limit quantities (surface/fields), their gradients, normals, and so on, given a barycentric coordinate. Given this initial set of tools, our algorithmic procedure is outlined in Alg. 2. 4.3.3 Complexity Analysis In order to accurately analyze the complexity, we forgo studying our system with respect to the traditional measure ๐‘๐‘ , total number of basis functions, and elect to cast the present scaling discussion in terms of the total number of unique quadrature points per patch. We choose to do so, because ๐‘๐‘ can obscure the true cost of setup and execution time, for a high order system, given the relatively large spatial extent of basis function support that can cause potential overlap. Furthermore, while asymptotic scaling is generally the metric of choice when discussing fast methods, one must be more careful when handling higher order methods as compared to lower order methods because the constants that appear in front of the asymptotic scalings are often significantly larger and can have a more pronounced impact on the relative cost of different formulations. Given the prior definitions, we make some assumptions that make asymptotic analysis tractable: these are (a) all basis functions are associated with the same number of patches 49 Algorithm 2: A description of the overall algorithm. 1 Define: 2 ๐’ฏ0 : control mesh 3 โ„Ž: Average edge length per patch 4 ๐‘„: Order of high order integration rule 5 โ„“ : level of refinement for integration scheme 6 ฮ”๐‘ฅ0 = โ„Ž/2โ„“ : initialize leaf box size 7 8 Initialize wideband MLFMA: 9 initialize FMM, construct oct-tree 10 Assign โ‰ˆ ๐‘„/4โ„“ integration points per subpatch 11 Assign centers of subpatches to leaf boxes 12 Construct near and far interaction list in terms of (sub) patches 13 14 isoBEM: 15 Compute and store partial sums in the nearfield using q.(4.5); ๐’ต ๐‘ ๐‘’ ๐‘Ž๐‘Ÿ = ๐’ต๐‘š๐‘› 16 Solve isoBEM system ๐’ตโ„ = ๐’ฑ using GMRES 17 Compute the matrix-vector product in GMRES using wideband MLFMA for remaining near-field contributions and far-field contributions 18 19 Scattering Cross Section: 20 Compute f(๐‘ฅ), f : ฮฉ โ†’ R3 , where ๐‘ฅ โˆˆ ฮฉ {Map surface pressure onto farfield} ๐œˆ, (b) the number of patches that share a nodal vertex is ๐œ‡, (c) the average size of the patch used is approximately 0.2๐œ† and (d) the number of quadrature points used in any patch is ๐‘๐‘ž . The overall costs comprises three distinct portions; (a) costs to evaluate ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ , (b) cost to evaluate a matrix vector product with ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ , and (c) with ๐’ต ๐‘“ ๐‘Ž๐‘Ÿ . Traditionally, the cost for filling ๐’ต costs ๐’ช(๐‘๐‘2 ๐‘๐‘ž2 ๐œˆ 2 ) The cost of a matrix vector product with ๐’ต is ๐’ช(๐‘๐‘2 ). To analyze costs associated the above scheme, lets start with traditional MLFMA. Assume that the leaf box size is 0.2๐œ†, and one assigns patches (and not basis) to each leaf box. This implies that there are ๐’ช(1) patches per leaf box or, equivalently, โ‰ˆ ๐‘๐‘ leaf boxes. This manner of constructing the tree implicitly assumes that some of the interaction between basis functions whose domain overlap, like that shown in Fig. 4.1, are constructed via MLFMA whereas other portions are constructed via direct integration. From these arguments, it follows that the cost of constructing ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ scales as ๐’ช(๐‘๐‘ ๐‘๐‘ž2 ). The cost for 50 each of the MLFMA steps scales as follows: ๐’ช(๐‘๐‘ ๐‘๐‘ž ) for C2M and L2O operations, and ๐’ช(๐‘๐‘ log ๐‘๐‘ ) for the rest of the tree operations. It is evident that both computing the near field as well as mapping on to the leaf level in a traditional MLFMA tree are the bottlenecks. The story is slightly different for wideband MLFMA. For simplicity assume that each patch is divided into ๐‘๐‘ž subpatches, and the size of each leaf box is such that one has ๐’ช(1) subpatches. This implies that there are approximately ๐‘๐‘ ๐‘๐‘ž leaf boxes. The cost of evaluating ๐’ต ๐‘›๐‘’ ๐‘Ž๐‘Ÿ scales as ๐’ช(๐‘๐‘ ๐‘๐‘ž ). The cost of C2M and L2O scales as ๐’ช(๐‘๐‘ ๐‘๐‘ž ) and traversal up and down the tree as approximately ๐’ช(๐‘๐‘ ๐‘๐‘ž log ๐‘๐‘ ๐‘๐‘ž ). As is apparent, the  key difference between the two methods is a factor of ๐‘๐‘ž in the cost for near field. This is significant as ๐‘๐‘ž needs to be high enough to integrate a higher order function. This section verifies the theoretical scaling for both complexity and storage, Figs. 4.6 and 4.7, in terms of total number of number of quadrature points per patch or ๐‘๐‘ ๐‘๐‘ž ; both experiments are conducted on a series of boxes of size 1.151๐œ† โˆ’ 9.212๐œ†. In Fig. 4.6, we demonstrate the scaling of the timing of a single matrix-vector product. The scaling is ๐’ช(๐‘๐‘ ๐‘๐‘ž ๐‘™๐‘œ ๐‘”(๐‘๐‘ ๐‘๐‘ž )), as expected. In our careful analysis, we find that it is evident that the scaling is primarily from the farfield contributions and in particular dominated by the M2M and L2L stages; this results is as expected due to our system being higher order, thus demanding a high level integration rule. 4.3.4 Accuracy of wideband MLFMA for nearfield computations As the rationale for using wideband MLFMA is to alleviate the computational complexity associated with nearfield computations, in what follows we examine the accuracy and efficacy of the proposed approach. The numerical experiments to evaluate the accuracy are conducted on a uniformly discretized sphere. As is done in traditional MLFMA, we first choose the largest size of the leaf box to be greater than the support of any basis function. Next, we choose the height of the tree such that all basis functions that share a supporting patch interact in the nearfield of each other. To fix the integration rule necessary to evaluate inner products, we choose the number of subdivisions per-patch 51 and use low-order quadrature on each sub-patch. Specifically, we assume that each patch is sub-divided into 16 sub-patches. The resulting matrix or matrix-vector product (with any random vector) forms the benchmark for our subsequent tests. Next, we increase the height of the tree. Now, we would have contributions to the matrix vector product from both near and โ€œfarโ€ interactions via ACE/MLFMA. Note, as ACE is low-frequency stable, it may be tempting to arbitrarily choose a very small leaf box. The optimal size of the leaf box is approximately the size of the largest sub-patch. Note, division of each patch is tantamount to defining a quadrature rule for each patch. This said, there are obvious tradeoffs. A very small leaf box implies that most of the computation is done in the farfield, i.e. numerically. This is inadvisable for hypersingular integrals as regularization within the FMM algorithm is not possible. Our scheme relies on the order of polynomials we need to evaluate over the surface (geometry and basis order), defining a quadrature rule, using a division of patches to achieve that rule, and then defining a leaf box size such that that number of near interactions are approximately ๐’ช(1). On the finest division of a patch, we use a 3 point rule. In our experiment, we choose a sphere of radius 2.63๐œ† and a tree with leaf box of size ฮ”๐‘ฅ0 = 0.4๐œ†; we note that each patch consists of 16 subpatches. We truncate the tree at the leaf level, level โ„“ and compute all nearfield interactions. Next, we choose the leaf box at ฮ”๐‘ฅ0 = 0.2๐œ†, which translates to roughly 4 subpatches per leaf box. We then compute all the nearfield interactions and one level of MLFMA farfield. Again, we adjust the leaf box size to ฮ”๐‘ฅ 0 = 0.1๐œ†, which allows 1 subpatch per leaf box and introduces two levels of MLFMA farfield. As all the tree computations using ACE, the number of terms in the expansion, ๐‘ƒ controls the accuracy. Fig. 4.4 demonstrates the computational speedup and Fig. 4.5 describes the controllable accuracy as a function of ๐‘ƒ. As is evident from these figures, one can control the accuracy and the corresponding speedup. 52 Figure 4.4: Comparison of isoBEM nearfield timings against the ๐’ช(๐‘๐‘2 ๐‘๐‘2 ๐‘๐‘ž2 ) direct fill algorithm on a 2.63๐œ† sphere, where P is the expansion coefficient. Figure 4.5: Convergence in relative ๐ฟ2 error of matrix using isoBEM versus the ๐’ช(๐‘๐‘2 ๐‘๐‘2 ๐‘๐‘ž2 ) direct fill algorithm on a 2.63๐œ† sphere, where P is the expansion coefficient. 53 Figure 4.6: Scaling of MLFMA matrix vector product timings for a series of boxes of size 1.151๐œ† โˆ’ 9.212๐œ†. 4.3.5 Computational Complexity and Storage of MFLMA BIE-subdivision Scaling in storage for the near and far fields are shown in Fig. 4.7. As expected, the nearfield memory scales linearly with respect to ๐‘๐‘ ๐‘๐‘ž . The farfield memory, for our given experiment, is mostly dominated by the precomputation of the aggregation/disaggregation and translation operator and as a result we find scaling to be ๐’ช(๐‘๐‘ ๐‘๐‘ž ๐‘™๐‘œ ๐‘”(๐‘๐‘ ๐‘๐‘ž )). 4.4 Numerical Examples In this section, we present a number of results validating our isoBEM with respect to accuracy and efficiency. Furthermore, we present solutions to scattering from conical and non-conical scatters and make the appropriate comparisons against analytic solutions and flatBEM; note, flatBEM is used to denoted a solution methodology that uses piecewise flat tesselations to represent the geometry, and constant basis functions to represent the potentials. In all cases analyzed here, the discretization of objects used for analysis via 54 Figure 4.7: Scaling in nearfield and farfield memory for a series of boxes of size 1.151๐œ† โˆ’ 9.212๐œ†. flatBEM is such that its surface area is approximately identical to that used for isoBEM (within 99%). Note, to achieve this agreement, one needs a higher number of flat tesselations to mimic the higher order nature of subdivision surfaces. Finally, we also remind the reader that the number of degrees of freedom for flatBEM corresponds to the triangles whereas those for isoBEM corresponds to vertices (more precisely, control points). Furthermore, for all examples in the proceeding sections, we note that the incident field propagates along kฬ‚ = โˆ’zฬ‚ direction. The scattering cross (SCS) section [46] is used to as a metric to compare simulations, unless it is specified otherwise. All the cases demonstrated below assume that the test objects are sound-hard and are immersed in a homogeneous medium. The speed of sound in the ambient medium is assumed to be 343 ๐‘š/๐‘ , and the dimensions of all objects analyzed in this chapter are normalized with respect to the wavelength of the incident field. Finally, when errors are reported, they are the ๐ฟ2 norm between computed and reference quantities. 55 Figure 4.8: Scattering cross section of the sphere (๐œ™ = 0 cut). 4.4.1 Accuracy of isogeometric analysis To validate and demonstrate the accuracy of the proposed approach, we consider scattering from a sphere discretized at multiple resolutions. To this end, we consider a sphere with radius of 1.0 meter that is modeled using an initial control mesh comprised of 642 vertices and 1280 faces. Starting from this control mesh and performing two levels of Loop subdivision results in two meshes; the coarser mesh with 2562 vertices and 5120 faces, and finer mesh having 10242 vertices and 20480 faces, respectively. All three meshes have the same limit surface. This implies that an increase in the number of control points is equivalent to โ„Ž-refinement as the support of the basis function changes and therefore, the approximation of the continuous operator โ„’ ๐ต๐‘€ {ยท}. The number of degrees of freedom (DoF) in isoBEM is equivalent to the number of vertices; consequently, the three tests have 642, 2562, and 10242 DOFs, respectively. By 56 comparing those obtained using isoBEM with analytic series approach, the relative ๐ฟ2 errors in surface pressure are 0.0591, 0.0180 and 0.00750, respectively, for these three different discretizations. As is evident, the accuracy in the potential in the ๐ฟ2 norm is excellent and it is convergent. Finally, Fig. 4.8 makes a direct comparison between farfield data obtained using isoBEM, on the sphere of one level of refinement corresponding to 5120 unknowns, and flatBEM on the sphere of two level of refinement corresponding to 20480 unknowns, and analytical solutions. We find there to be excellent agreement between all three data sets. Figure 4.9: Scattering cross section of a 8.1๐œ† by 8.1๐œ† by 2.3๐œ† torus (๐œ™ = 0 cut). 4.4.2 Analysis of non-conical geometries Next, to demonstrate the viability of using isoBEM we compare results obtained using this method against those obtained using flatBEM. In all cases analyzed here, the discretization of objects used for analysis via flatBEM is such that its surface area is approximately identical to that used for isoBEM (within 99%). Note, to achieve this agreement, one needs a higher number of flat tesselations to mimic the higher order nature 57 Figure 4.10: Scattering cross section of a 5.71๐œ† by 11.66๐œ† by 5.71๐œ† Cone (๐œ™ = 0 cut). Figure 4.11: Scattering cross section of a 4.57๐œ† by 13.80๐œ† by 4.57๐œ† warhead (๐œ™ = 0 cut). 58 of subdivision surfaces. Finally, we also remind the reader that the number of degrees of freedom for flatBEM corresponds to the triangles whereas those for isoBEM corresponds to vertices (more precisely, control points). We start with analyzing scatter from a torus that fits within a box of size 8.1๐œ† by 8.1๐œ† by 2.3๐œ†. The number of degrees of freedom for flatBEM is 61312 whereas those for isoBEM is 7664. As is evident from Fig. 4.9 the agreement between the two is excellent. Next, we consider a cone that fits in a 5.71๐œ† by 11.66๐œ† by 5.71๐œ† bounding box. Again, flatBEM uses 59424 degrees of freedom while the isoBEM has 7430 degrees of freedom; Fig. 4.10 demonstrates the agreement between the two. Finally, we consider a warhead that fits in a 4.57๐œ† by 13.80๐œ† by 4.57๐œ† bounding box. The flatBEM uses 85504 degrees of freedom whereas that using isoBEM uses 10690. As in all the previous cases, the two sets of data is shown in Fig. 4.11 where we see that they agree very well with each other. 4.5 Conclusion In this chapter, we detail an implementation of an isogeometric anlaysis tool for scattering from rigid bodies using the Burton-Miller formulation. The basis sets we have chosen arise from subdivision representation of the geometry; the benefits of this representation are that the resulting representation are ๐ถ 2 almost everywhere. Unfortunately, this high degree of regularity implies that high order functions are used for representing both the geometry and physics. This translates to high numerical quadrature that can overwhelm overall costs. Indeed this has been a persistent challenge for all higher order methods. Our secondary contribution in this chapter is extension of the wideband multi-level fast multipole method to efficiently evaluate both components of the matrices and inner products with them. Through numerous numerical results, we have shown the efficacy of this approach in terms of accuracy against analytical data, comparison against piecewise flat triangulation for non-conical objects, and cost complexity for acoustically large objects. Overall, isogeometric analysis tools with the necessary mathematical augmentation permit efficient analysis 59 directly from a CAD representation. This chapter, ยฉ2020 JASA, is reprinted, with permission, from A. M. A. Alsnayyan, J. Li, S. Hughey, A. Diaz, and B. Shanker, "Efficient isogeometric boundary element method for analysis of acoustic scattering from rigid bodies," The Journal of the Acoustical Society of America on, May 2020, with slight modifications to fit in this dissertation. 60 CHAPTER 5 AN EFFICENT ELECTROMAGNETIC ISO-GEOMETRIC INTEGRAL EQUATION SOLVER 5.1 Maxwellโ€™s Equations Over the past six decades, the state of the art boundary integral equation solvers have grown by leaps and bounds to become a powerful tool for electromagnetic analysis. A sequence of advancements have enabled this transition, starting from the development of integral equations (see [127] and references therein for a more complete historical background), to methods to appropriately discretize them [130], to higher order representa- tions [71], to overcoming computational bottlenecks [161, 147, 51, 112], to well-conditioned formulations [30, 40, 48], and more recently, to preconditioning techniques [13, 2]. However, despite the significant recent progress made, the technological drivers demand a more sophisticated and more feature rich solver, albeit at reduced cost. Computational analysis typically proceeds in three stages: (a) construct a geometric model using a computer aided design (CAD) tool, (b) define a discrete representation of the model, and (c) finally, choose a representation of the physics on the discrete representation of the geometry. Geometry is typically represented using bi-variate splines (Bezier splines, B-splines, or non-uniform rational B-splines (NURBS)) that can provide higher order continuity on the surface. From this surface representation, a mesh is generated that typically provides low order continuity on the manifold. As an example, piecewise flat Lagrangian elements are ๐ถ 0 , i.e., continuous at interfaces between triangular tesselations of control nodes, or patches, but with discontinuous normals. Furthermore, even higher order meshes are higher order within a patch/subdomain, but still ๐ถ 0 across patches. As a result, basis functions defined on these meshes must impose additional constraints. In this framework, a number of different approaches to electromagnetic analysis tools have been developed, including: RWG basis sets [130], its higher order variants [71], Buffa-Christansen basis [35]. In addition, there exists an in-depth analysis and study into 61 computational bottlenecks such as ill-conditioning, low-frequency breakdown, dense-mesh breakdown, topological breakdown, etc [48, 78, 13]. Two more relatively recent methods take a different approach; they still seek to obtain a higher order parameterization of the geometry and thereby, higher order basis for physics. The first overcomes item (a) above, and directly models the object using higher order polynomials [120]. Another approach, the generalized method of moments (GMM), starts with (b) and builds a framework that accommodates both large (> 4๐œ†) and small patches as well as different functions on each patch [51, 97] all stitched together within a partition of unity framework. This is done using a non-watertight representation of standard meshes. Other methods rely on different techniques to enrich function spaces to represent physics (for instance, macro-basis sets [158]). All seek to achieve an efficient representation of geometry, physics, or both. An alternative approach that is gaining currency is equipped with the infrastructure to do physics using the same basis function used to construct the geometry; this is known as isogeometric analysis (IGA). The advantages of such an approach are as follows: they (a) eliminate the error in translating between geometry and the mesh, (b) the number of degrees of freedom is limited to that used for geometry representation which is significantly smaller than a corresponding mesh, and (c) the rules used for adaptation and refinement are identical for both geometry and physics; a vivid illustration can be found in [35, 146, 97, 98]. One must highlight that in isogeometric methods, basis functions are co-located on control nodes used to describe the geometry. This is in contrast with parametric methods that require additional infrastructureโ€“for an example of using subdivision for geometry and GMM basis sets, see [97]. The genesis of IGA methods started with using NURBS for solid mechanics [79], and more recently, in electromagnetics [35, 56] and acoustics [146]. NURBS are geometric descriptions that are topologically either a disk, a tube or a torus. As a result, one of the major difficulties that arise with NURBs is that the patches have to be seamlessly sewed 62 together in order to handle complex surfaces which is time-consuming and complicated. Furthermore, stitching together these patches can result in surfaces that are not watertight and sometimes discontinuous. These complexities have proven to be quite a hindrance when handling complex geometries [23]. Other modalities that have gained currency in geometry representation are T-splines and Loop subdivision. While T-splines have been used in an IGA setting (see [23, 59] and references therein) our focus in this paper will be on Loop subdivision. Loop subdivision have been extremely popular in the computer graphics industry due to the ease with which one can represent complex topologies, its scalability, inherently multiresolution features, efficiency and ease of implementation. More importantly, the surface representation is ๐ถ 2 , or continuous twice differentiable surface, almost everywhere making it an attractive candidate for defining physical basis sets as it avoids the requirement of defining additional mathematical framework that is commonplace in other low order basis set [78, 130, 37]. There has been a concerted effort to develop IGA methods on subdivision surfaces in a number of fields, including electromagnetics [97, 98, 67], acoustics [6, 146] and shape reconstruction/optimization [9, 18, 20, 155, 39] This chapter builds on our earlier body of work on Loop subdivision based IGA for the electric field integral equations [98] and construction of Debye sources [67]. In both these cases, the objects analyzed were simply connected and electrically small. Further, they only discretized the electric field integral equation. The key bottleneck is the number of quadrature points required to evaluate all necessary inner products on higher order geometry (4๐‘ก โ„Ž ) and 3๐‘Ÿ๐‘‘ order basis. A principal goal of this paper is to alleviate this bottleneck for all methods that use higher order surface representation and higher order basis for physics. It is illustrated here for subdivision basis. To do so, we exploit wideband multilevel fast multipole algorithm to evaluate all interactions (self, near, and far) with leaf boxes as small as 0.025๐œ†. Next, we introduce the MHB, for electromagnetic field computation; here we briefly highlight some of the features found in chapter 3. The MHB 63 defined on the Loop subdivision surface is attractive as it inherits properties that are of interest when representing physical quantities on a surface: high-order continuity and multiresolution. For instance, given that a hierarchy of subdivision meshes point to same limit surface, it follows that if a function is associated with control points, then this will follow the refinement rules of subdivision. That is, one can compose a MHB on the limit surface via a composition of those computed on a series of refined control meshes. Of course, one needs the prolongation and restriction operators [144, 19] for communication of data defined on different levels of refinement. Furthermore, the MHB provides a smooth, orthogonal global basis set. While this does provide an excellent representation on objects that are sufficiently smooth, the number of basis sets can be high if on desires higher fidelity near sharp features (say tips). But it is possible to incorporate localised enrichment based on known local physics or heuristics, in tandem with the standard basis set. Such examples have been shown for polynomial basis sets [140, 17]. Here we explore development of both these types of basis sets and application to problems in electromagnetics, specifically to compress systems resulting from discretization of boundary integral equations in electromagnetics, and demonstrate its numerous benefits. Lastly, our previous work on IGA on subdivision surfaces has taken two different paths: (a) defining currents via a Helmholtz decomposition on simply connected surfaces [98, 7] and (b) defining a div-conforming basis sets [101]. Using the former, we have been able to demonstrate numerous benefits of smoothness and thereby a complete Helmholtz decomposition on simply connected surfaces. This permits seamless implementation of Calderรณn operators [98, 7], Debye sources [67, 62] and scalar integral equations [100]. They show excellent promise due to their accuracy. Exact Helmholtz decomposition on simply connected surfaces together with Calderon preconditioner[98, 7]. What is critically missing from this suite is analysis of multiply connected structures. The harmonic basis functions or global loops on multiply connected structures are not easily handled within the purview of defining basis sets via surface derivatives, as they are both divergence-free and curl-free. 64 One alternative is to avoid the need to define these decompositions by using divergence conforming basis [101]. But, we no longer will be able to take advantage of the properties described in the previous paragraph. Alternatively, we could use divergence conforming functions together with exact Helmholtz decomposition in a manner akin to that prescribed in [13] (i.e., use projectors), but at the expense of the number of degrees of freedom. In this paper, we present an alternate method for defining a harmonic basis by using vector fields wherein each component is a randomized polynomial [4]. In chapter 3 we demonstrated that (a) the representation is exact in that the basis sets are divergence-free, curl-free, and both divergence-free and curl-free, (b) they are orthogonal to each other, and (c) the representation captures analytical fields defined on a torus. Here, we will we will use this basis to discretize a regularized CFIE formulations for analysis of electromagnetic scattering. To this end, the principal contributions of this chapter are as follows: we present (a) scattering analysis on simply connected and multiply connected Loop subdivision surfaces using a Helmholtz decomposition, (b) formulate and implement the H-MHB for multi-resolution scattering analysis on complex Loop subdivision surfaces, and (c) formulate and implement the L-MHB for scattering analysis on Loop subdivision surfaces with challenging features; furthermore, we pair these approach with a well-conditioned combined field integral equation to analyze objects as large as 120๐œ† using the wideband multilevel fast multipole algorithm. 5.2 Electromagnetic Boundary Integral Equations We consider the analysis of scattered fields {E๐‘  , H๐‘  }, from a perfect electrically conduct- ing (PEC) object ฮฉ, due to fields {E๐‘– , H๐‘– } incident on its boundary ฮ“ โˆˆ ฮฉ. It is assumed that this surface is equipped with a unique outward pointing normal denoted by nฬ‚(r), r โˆˆ ฮ“. The region external to this volume {R3 \ ฮฉ} is occupied by free space. In the exterior region {R3 \ ฮฉ}, the total fields {E, H} = {E๐‘  + E๐‘– , H๐‘  + H๐‘– }, are subject to the Maxwell equations Instead of solving Maxwellโ€™s equations directly, we use a mixed potential formula 65 and solve for the surface current density J. The scattered wave (E๐‘  , H๐‘  ) is obtained using equivalence theorems leading to the following: nฬ‚(r) ร— E๐‘  (r) = ๐’ฏ๐œ… โ—ฆ J(r), (5.1) ๐‘  nฬ‚(r) ร— H (r) = ๐’ฆ๐œ… โ—ฆ J(r), where, โˆซ ๐’ฏ๐œ… โ—ฆ J(r) = โˆ’๐‘—๐œ‚๐œ…nฬ‚(r) ร— ๐บ๐œ… (r, r0) ยท J(r0)๐‘‘r0 โˆซฮ“ (5.2a) ๐œ‚ 0 0 0 0 + ๐‘— nฬ‚(r) ร— โˆ‡ ๐บ๐œ… (r, r )โˆ‡ ยท J(r )๐‘‘r , ๐œ… ฮ“ โˆซ ๐’ฆ๐œ… โ—ฆ J(r) = nฬ‚(r) ร— โˆ‡ ร— ๐บ๐œ… (r, r0) ยท J(r0)๐‘‘r0 , (5.2b) ฮ“ where ๐บ๐œ… (r, r0) = exp[โˆ’๐‘—๐œ…|r โˆ’ r0 |]/(4๐œ‹|r โˆ’ r0 |), ๐œ… is the free space wavenumber, ๐œ‚ is the free space impedance, and J(r0) is the equivalent current that is induced on the surface. In the above expressions, and what follows, we assume and suppress exp[๐‘—๐œ”๐‘ก] time dependence. Now utilizing the property that the tangential component of the electric field is continuous using and the tangential component of the magnetic field is discontinuous, we arrive at the requisite electric field and magnetic field integral equations (EFIE/MFIE) as   ๐‘– ๐‘  (5.3a) EFIE: = nฬ‚(r) ร— nฬ‚(r) ร— E (r) + E (r) = 0,   MFIE: = nฬ‚(r) ร— H๐‘– (r) + H๐‘  (r) = 0. (5.3b) For so-called irregular frequencies both the EFIE and MFIE suffer from non-unique solutions, but their linear combination yields a uniquely solvable formulation throughout the frequency spectrum denoted as the combined field integral equation (CFIE): (1 โˆ’ ๐›ผ)(โ„ โˆ’ ๐’ฆ๐œ… ) โ—ฆ J + ๐›ผ nฬ‚ ร— ๐’ฏ๐œ… โ—ฆ J = (5.4) ๐‘– ๐‘– (1 โˆ’ ๐›ผ)nฬ‚ ร— H โˆ’ ๐›ผnฬ‚ ร— nฬ‚ ร— E , 66 where ๐›ผ is a positive constant. It is well known that CFIE suffers from several breakdowns (low frequency, dense mesh, topology, etc.) [32, 30, 48]. There has been an extensive body of literature addressing these bottlenecks [98, 13]. In the next section we aim to employ our basis that completely satisfy Helmholtz decomposition to be used in a Calderรณn setting. We introduce a regularized CFIE formulation that overcomes several of these breakdowns including the non-uniqueness problem. In what follows, we detail this regularized CFIE. 5.2.1 Regularized Combined Field Integral Equations A regularized reformulation of (5.4) is typically based on the construction of the regularizing operators dervied from Calderรณn identities and complexification techniques. Such operators โ„› ๐œ… have been proposed and analyzed in the literature [48, 32, 30, 50, 114, 54]. This body of work the CFIER written as follows:   โ„ โˆ’ ๐’ฆ๐œ… โ—ฆ J + โ„› ๐œ… โ—ฆ ๐’ฏ๐œ… โ—ฆ J = nฬ‚ ร— H๐‘– โˆ’ โ„› ๐œ… โ—ฆ (nฬ‚ ร— E๐‘– ). (5.5) 2 Here, โ„› ๐œ… is chosen as a regularizing operator for ๐’ฏ๐œ… such that the integral operators on the left hand side of (5.5) are second kind Fredholm operators. In particular, the choice of regularization operators is provided in [30]. This formulation was found to showcase the superior performance of solvers based on the novel Calderรณn-Complex CFIER (CC-CFIER) formulations that involve the boundary integral operators   โ„ โˆ’ ๐’ฆ๐œ… โ—ฆ J โˆ’ 2๐’ฏ๐œ…0 โ—ฆ ๐’ฏ๐œ… โ—ฆ J = nฬ‚ ร— H๐‘– + 2๐’ฏ๐œ…0 โ—ฆ (nฬ‚ ร— E๐‘– ), (5.6) 2 where ๐œ…0 = ๐œ… โˆ’ ๐‘—0.4๐œ2/3 ๐œ…1/3 and ๐œ is the maximum of the absolute values of mean curvatures on surface ฮ“. 5.2.2 Field Solvers on Loop subdivision surfaces In this section we will detail the discretization of (5.6) on multiply connected Loop subdivision surfaces. We note that all basis function are defined for all values of r โˆˆ ฮ“(r) and the harmonic basis have global support. In particular, we use a Galerkin scheme to 67 discretize these equations. We note that efficient evaluation of inner products is discussed in the next section. Note that discretizing Calderรณn type operators requires intermediate spaces, effected through a Gram matrix. We define the required Gram-matrix [๐บ] using ๏ฃฎ 11 ๏ฃฏ๐’ข 0 0 0 ๏ฃบ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ 0 ๐’ข ๏ฃฏ 22 0 0 ๏ฃบ ๏ฃบ [๐บ] = ๏ฃฏ ๏ฃฏ ๏ฃบ, (5.7a) 33 ๐’ข 34 ๏ฃบ ๏ฃบ ๏ฃฏ 0 ๏ฃฏ 0 ๐’ข ๏ฃบ ๏ฃฏ 43 ๐’ข 44 ๏ฃบ ๏ฃบ ๏ฃฏ 0 ๏ฃฐ 0 ๐’ข ๏ฃป โˆซ ๐‘™๐‘˜ ๐‘˜ ๐บ ๐‘›๐‘š = hJ๐‘›๐‘™ (r), J๐‘š ๐‘˜ (r)i = J๐‘›๐‘™ (r) ยท J๐‘š (r)๐‘‘r, (5.7b) ฮ“ where ๐‘™, ๐‘˜ โˆˆ {3, 4} and ๐บ 33 ๐‘›๐‘š = ๐บ ๐‘›๐‘š , ๐บ ๐‘›๐‘š = โˆ’๐บ ๐‘š๐‘› . Using (2.7) in (5.6) and Galerkin testing 44 34 43 results in a matrix system that can be written as [๐‘] [๐ผ] = [๐‘‰] (5.8a) where, [๐‘] = [๐บ]โˆ’1 [[๐ฟ] + [๐พ]] (5.8b) with J๐‘˜   ๐‘™๐‘˜ [๐พ]๐‘›๐‘š = J๐‘›๐‘™ (r), ๐‘š (r) โˆ’ ๐’ฆ๐œ… โ—ฆ J๐‘š๐‘˜ (r) , (5.8c) 2 ฮ“๐‘› [๐‘‡]e๐œ…๐‘™ ๐‘˜,๐‘›๐‘š = J๐‘›๐‘™ (r), ๐’ฏe๐œ… โ—ฆ J๐‘š ๐‘˜ (r) ฮ“๐‘› , (5.8d) where e๐œ… โˆˆ {๐œ…0 , ๐œ…}, and, as defined earlier ๐œ…0 = ๐œ… + 0.4๐œ2/3 ๐œ…1/3 , and ๐œ is the maximum of the magnitude of the mean curvature of the object, and [๐ฟ] = โˆ’2[๐‘‡]๐œ…0 [๐บ]โˆ’1 [๐‘‡]๐œ… . (5.8e) Furthermore, we have ๐‘˜ ๐‘˜ [๐ผ]๐‘š = ๐‘Ž๐‘š , (5.9a) [๐‘‰]๐‘›๐‘˜ = [๐บ]โˆ’1 โˆ’2[๐‘‡]๐œ…0 [๐บ]โˆ’1 [๐‘‰๐‘‡ ]๐‘›๐‘˜ + [๐‘‰๐พ ]๐‘›๐‘˜ ,   (5.9b) 68 with [๐‘‰๐‘‡ ]๐‘›๐‘˜ = J๐‘›๐‘˜ (r), E๐‘– (r) ฮ“๐‘› , (5.10a) [๐‘‰๐พ ]๐‘›๐‘˜ = J๐‘›๐‘˜ (r), nฬ‚ ร— H๐‘– (r) ฮ“๐‘› , (5.10b) Here, we have defined ๐‘™, ๐‘˜ โˆˆ {1, 2, 3, 4}. Lastly, we note that the stabilizing properties of the Calderรณn preconditioner are local [1], which allows the use of a localized version of the preconditioner [๐‘‡]๐œ…0 . As such, we choose to omit all interactions of a distance greater than 1.25๐œ†. 5.2.2.1 Manifold Harmonic Transform Now we turn our attention to discretizing (5.6) using the MHB; this basis would be excellent candidate to create a reduced order representation of currents. Consider a reduced ๐‘€ orthogonal MHBs that span ๐‘Š ๐‘€๐ป (ฮ“) โŠ‚ ฮจ(ฮ“). This is tantamount to using ๐‘€ < ๐‘›๐‘ฃ for both the representation and measurement space in (3.14). As a result, one obtains a compressed impedance matrix. As such we can reconstruct (5.8a) as [๐ป][๐‘ ๐ป ][๐ผ๐ป ] = [๐ป][๐‘‰๐ป ] (5.11) where, [๐‘ ๐ป ] = [๐ฟ๐ป ] + [๐พ ๐ป ] , (5.12a) [๐ผ๐ป ] = [๐ป]๐‘‡ [๐ผ] , (5.12b) h i ๐‘‡ ๐‘‡ [๐‘‰๐ป ] = [๐ป] โˆ’2[๐‘‡]๐œ…0 [๐ป] [๐ป] [๐‘‰๐‘‡ ] + [๐‘‰๐พ ] , (5.12c) and, [๐ฟ๐ป ] = โˆ’2[๐ป]๐‘‡ [๐‘‡๐œ…0 ][๐ป][๐ป]๐‘‡ [๐‘‡๐œ… ][๐ป], (5.12d) [๐พ ๐ป ] = โˆ’2[๐ป]๐‘‡ [๐พ][๐ป]. (5.12e) Here [๐ป]๐‘‡ denotes the transpose of [๐ป]. 69 5.3 Wideband MLFMA for Evaluation of Inner Products The evaluation of the aformentioned inner products and matrix vector products are challenging due to the domain of support of each loop basis function being electrically large and are on average โ‰ˆ 0.9๐œ†. Furthermore, we note that the basis functions are third order and the geometry is fourth order. Both serve to exacerbate costs as one needs higher order quadrature rules over both test and source domains. To ameliorate these, we exploit the wide-band FMM introduced by the authors in [112, 51]. The framework we propose has been used to accelerate matrix evaluations as well as matrix vector products for the Generalized Method of Moments (GMM) wherein patch sizes can be several wavelengths long [51] using a mixed potential formulation. In order to motivate this framework we begin by briefly examining the fundamental rubrics involved in evaluating all inner products. As has been shown in [101], expected convergence rates for for the EFIE are obtained, albeit with RWG type basis functions, provided one uses an adaptive quadrature based on subdividing each parent triangle. This implies that it is possible to choose a refinement level such that one can use a piecewise flat approximation for evaluation of singular integrals. But the downside is clear; the cost of evaluating the inner-products is high. Indeed, if there are ๐‘๐‘ quadrature points per basis, the cost for evaluation of of matrix interactions scale as ๐’ช(4๐‘๐‘2 ๐‘๐‘ฃ2 ). Our goal is to reduce this to ๐’ช(2๐‘๐‘ ๐‘๐‘ฃ log ๐‘๐‘ ๐‘๐‘ฃ ).  In order to mitigate the high computational cost we employ MLFMA. We note that the following is true for all methods that use higher order modeling. To set the stage for the ๐‘™๐‘˜ discussion, consider a matrix element [๐‘‡]๐œ…,๐‘›๐‘š หœ . It comprises contributions from both the magnetic vector potential and electric scalar potential. Let us focus on the latter, specifically, just evaluating the scalar potential due to J1๐‘› , โˆ’๐‘— โˆซ ฮฆ(r) = โˆ‡0 ๐บ๐œ…หœ (r, r0) ยท J1๐‘› (r0)๐‘‘r0 . (5.13) ๐œ”๐œ€ ฮ“๐‘› Consider a 1-level MLFMA prescription for an alternative evaluation of (5.13). We denote the center of a leaf box by r๐‘  ; at observations points sufficiently far away, where 70 X = r โˆ’ r๐‘  , the potential ฮฆ(r) can also be evaluated using โˆ’๐‘— ๐œ…หœ 2 โˆซ     ฮฆ(r) = kฬ‚ ยท โ„ณ ๐œ…หœ kฬ‚, r๐‘  ๐’ฏ ๐œ…หœ kฬ‚, X ๐‘‘ 2 kฬ‚, (5.14a) 16๐œ‹2 ๐œ”๐œ€ ๐‘†2 where   โˆซ J1๐‘› (r0)๐‘’ โˆ’๐‘— ๐œ…หœ kฬ‚ยท(r๐‘  โˆ’r ) ๐‘‘r0 , 0 โ„ณ ๐œ…หœ kฬ‚, r๐‘  = (5.14b) ฮ“๐‘› is the source to multipole map, and   ร•โˆž   (โˆ’๐‘—)๐‘› (2๐‘› + 1)โ„Ž ๐‘› (๐œ…๐‘‹)๐‘ƒ (2) ๐’ฏ ๐œ…หœ kฬ‚, X  หœ ๐‘› kฬ‚ ยท Xฬ‚ , (5.14c) ๐‘›=0 is the translation operator. Here, ๐‘† 2 denotes the unit sphere, parametrized by (๐œƒ, ๐œ™) โˆˆ [0, ๐œ‹] ร— [0, 2๐œ‹]. We note that kฬ‚ = kฬ‚(๐œƒ, ๐œ™). Note, the gradient on the Greenโ€™s function is evaluated spectrally; furthermore, a traditional approach to using MLFMA would ensure that the entire support domain ฮ“๐‘› lies within a leaf box, i.e., for loop basis functions the size of the leaf box is ฮ”0 โ‰ˆ 0.9๐œ†. This means that that each leaf box has approximately ๐‘๐‘ quadrature points, and the cost of computing the near field interactions in the MLFMA scheme is ๐’ช(4๐‘๐‘2 ๐‘๐‘ฃ2 ). As ๐‘๐‘ is relatively high, this still untenable. We would like ฮ”0 to be as small as possible, such that it contains far fewer quadrature points. Consider instead Fig. 5.1 which shows three leaf level boxes within ฮ“๐‘› . Furthermore, assume that we have to compute the self interaction of basis ๐‘› using MLFMA, boxes (1,3) are in the far field of each other, and (2,3) and (1,2) are in the near field of each other. Lets re-examine the evaluation of ฮฆ(r) for r โˆˆ ฮ“3,๐‘› . Using (5.13) โˆ’๐‘— โˆซ ฮฆ(r) = โˆ‡0 ๐บ๐œ… (r, r0) ยท J1๐‘› (r0)๐‘‘r0 ๐œ”๐œ€ ฮ“1,๐‘› ๐‘— โˆซ โˆ’ โˆ‡0 ๐บ๐œ… (r, r0) ยท J1๐‘› (r0)๐‘‘r0 , (5.15) ๐œ”๐œ€ ฮ“2,๐‘› = ฮฆ13 (r) + ฮฆ23 (r). Here, ฮ“๐‘–,๐‘› denotes the intersection of box-๐‘– with ฮ“๐‘› . There are two possible ways of evaluating the farfield interaction, ฮฆ13 (r), using a variation of (5.14) with the understanding 71 that the the domain of integration in (5.14b) is confined to ฮ“1,๐‘› . Specifically, โˆ’๐‘— ๐œ…หœ 2 โˆซ     ฮฆ13 (r) = kฬ‚ ยท โ„ณ 1 ๐œ… หœ kฬ‚, r ๐‘  ๐’ฏ ๐œ… หœ kฬ‚, r โˆ’ r ๐‘  ๐‘‘ kฬ‚ 2 16๐œ‹2 ๐œ”๐œ€ ๐‘†2 (5.16a) ๐œ…หœ โˆซ     = โ„ณ 2 ๐œ… หœ kฬ‚, r ๐‘  ๐’ฏ ๐œ…หœ kฬ‚, r โˆ’ r ๐‘  ๐‘‘ kฬ‚ 2 16๐œ‹2 ๐œ”๐œ€ ๐‘†2 where   โˆซ ห† J1๐‘› (r0)๐‘’ โˆ’๐‘— ๐œ…หœ ๐‘˜ยท(r๐‘  โˆ’r ) ๐‘‘r0 , 0 โ„ณ1 ๐œ…หœ kฬ‚, r๐‘  = (5.16b) ฮ“1,๐‘›   โˆซ โˆ‡0 ยท J1๐‘› (r0)๐‘’ โˆ’๐‘— ๐œ…หœ kฬ‚ยท(r๐‘  โˆ’r ) ๐‘‘r0 0 โ„ณ2 ๐œ…หœ kฬ‚, r๐‘  = โˆ’ ฮ“1,๐‘› โˆซ (5.16c) J1๐‘› (r0)๐‘’ โˆ’๐‘— ๐œ…หœ kฬ‚ยท(r๐‘  โˆ’r ) ๐‘‘r0 , 0 0 + uฬ‚ฮ“1,๐‘› (r ) ยท ๐œ•ฮ“1,๐‘› and uฬ‚๐œ•ฮ“ is outward pointing normal to the boundary ๐œ•ฮ“. Consider next, the near field evaluation of ฮฆ23 (r). As we want the leaf box size to be small, the minimum distance between the box centers becomes very small, and as such the order of singularity, due to the gradient on the Greenโ€™s function, introduces near-singular integration challenges. The remedy that is typically taken is to transfer the derivative onto the basis function. Specifically, ๐‘— โˆซ ฮฆ23 (r) = ๐บ๐œ… (r, r0)uฬ‚๐œ•ฮ“2,๐‘› (r0) ยท J1๐‘› (r0)๐‘‘r0 ๐œ”๐œ€ ๐œ•ฮ“2,๐‘› (5.17) ๐‘— โˆซ โˆ’ ๐บ๐œ… (r, r0)โˆ‡0 ยท J1๐‘› (r0)๐‘‘r0 . ๐œ”๐œ€ ฮ“2,๐‘› The lessons we take from the above equations are as follows: (a) the aforementioned line integrals have to be accounted for in (5.17) as they are implicitly included in (5.16b) and should cancel on the shared interface; (b) unfortunately, finding intersections between higher order surfaces and boxes is non-trivial; (c) with all challenges considered, we have to use (5.17) to evaluate ฮฆ23 (r) and not (5.15); Alternatively, it can be proven that interior line integral should vanish. This implies that an ideal choice would be to use (5.16c) and (5.17) sans line integrals in both to evaluate ฮฆ13 (r) and ฮฆ23 (r), respectively. Note, this example is illustrative. Further complication arises in the evaluation of the electric field, as it calls for the gradient of the scalar potential. As a result, one needs additional line integrals to 72 reduce the singularity. As is evident from the above discussion, using a mixed potential formulation together with wideband MLFMA permits evaluation of all integrals, near and far, without the consideration of the troublesome line integrals, but at the cost of more tree traversals [112]. Indeed, the size of the leaf box can now be as small as computationally expedient. Leaf box sizes can be chosen such that it contains ๐’ช(1) quadrature points, reducing the cost of near field evaluation to ๐’ช(2๐‘๐‘ ๐‘๐‘ฃ ). We elucidate this process by applying it to equation (5.8d), as depicted in Fig. 5.1; (5.8d) contains four independent terms that must be computed in the inner integral: three corresponding to the vector potential, and one corresponding to the scalar potential. It follows that any matrix element can be computed in terms of its partial contributions such that ๐›พ ร•ร• D E [๐‘‡]e๐œ…๐‘™ ๐‘˜,๐‘›๐‘š = โˆ’๐‘—๐œ”๐œ‡ J๐‘›๐‘™ , nฬ‚ (r) ร— ๐’ฎe๐œ… โ—ฆ J๐‘š ๐‘˜ ฮ“๐‘›๐œ ๐œ ๐›พ (5.18) ๐‘—๐›ฟ ๐‘™2 ๐›ฟ 1๐‘˜ D ๐›พ E + โˆ‡ฮ“ ยท J๐‘› , ๐’ฎe๐œ… โ—ฆ โˆ‡ฮ“ ยท J๐‘š ๐œ , 1 1 ๐œ”๐œ– ฮ“๐‘› where, โˆซ ๐›พ ๐‘˜ ๐‘˜ ๐’ฎe๐œ… โ—ฆ J๐‘š = ๐บ๐œ…หœ (r, r0)J๐‘š (r0) ๐‘‘r0 , (5.19) ๐›พ ฮ“๐‘š and the indices ๐œ and ๐›พ are subpatches of ฮ“๐‘š and ฮ“๐‘› ; subpatches within each otherโ€™s farfield are constructed via MLFMA whereas nearfield patches are constructed via direct integration, see Fig. 5.1. 5.3.1 Accuracy of wideband MLFMA for adaptive interactions Herein, we study the the accuracy of using wideband MLFMA is to alleviate the computational complexity associated with nearfield computations. To test the controllable accuracy of the aforementioned scheme, we conduct a controlled test. As an aside, the support of a basis function is a one-ring associated with a control vertex, three basis functions are defined on a patch. The most efficient assembly of interactions is computing these in a patchwise manner. To that end, consider two patches that share an edge. The edge length of each patch is approximately 0.25๐œ†. We compute the patch to patch interaction by 73 Figure 5.1: The support of two basis (no. 7 and no. 8) is shown in green and purple, respectively. In computing the interaction between the two, patch (triangular tesselation of control nodes) is partitioned into sub-patches (shown in red) to create an adaptive quadrature. The entire object is embedded in an MLFMA tree, and the size of a leaf box (shown within the green region) is about the size of a sub-patch. using an integration rule developed by subdividing each patch into 16 sub-patches and using a 3-point rule in each. Note, we are not computing self-patch interactions. Next, we compute the same interaction, but through an MLFMA tree with leaf box sizes ฮ”0 = 0.125๐œ† and ฮ”0 = 0.0625๐œ† which results in 1-level and 2-level computation of the interactions. The standard tree partitioning of interactions is used; the leaf box of size 0.125๐œ† has about 4 subpatches, whereas 0.0625๐œ† has approximately 1 subpatch. Given the size of leaf box, interactions are computed using wideband MLFMA which invokes Accelerated Cartesian Expansions (ACE) for leaf box sizes smaller than 0.2๐œ†. Fig. 5.2 demonstrates the controllable accuracy of computing these interactions as a function of ๐‘, the expansion coefficient for our wideband MLFMA scheme. As is evident from this figure, one can control the accuracy to very fine precision. 5.4 Numerical Examples In this section, we present a collection of numerical results to demonstrate the efficacy of the proposed approach. The first two main contributions are (a) subdivision based isogeometric formulation for simply connected objects, and (b) employing manifold 74 Figure 5.2: Convergence in relative ๐‘™2 error of partial matrix element using wideband MLFMA versus the direct fill algorithm, where ๐‘ is the expansion coefficient. harmonics for EM analysis. To this end, the data presented in this section highlights the following: (i) the accuracy of the two proposed approaches when compared against analytical data; (ii) the improved spectral properties of the CC-CFIER by means of the reduced numbers of iterations required for convergence of the GMRES iterative solver for Loop and MHB, (iii) the high-accuracy and reduced DOF under the MHB, and (iv) application of both to analyzing complex targets. In following set of examples, we focus on numerical examples to demonstrate the convergence and effectiveness of the presented numerical scheme to electromagnetic analysis for (a) application of our wideband MLFMA scheme for analyzing complex multiply connected structures; (b) timings and iteration count for overall GMRES iterative solve. Lastly, we seek to employing the H-MHB and L-MHB for EM analysis on Loop subdivision surfaces through numerical examples. We begin with (a) scattering anaylsis from complex targets using H-MHB and L-MHB, and (b) the high-accuracy, reduction in number of degrees of freedom, and improved computational efficiency using H-MHB and L-MHB, relative to MHB. 75 Unless otherwise stated, we compute scattering due to a plane wave field propagating in ๐œ…ห† = โˆ’ห†๐‘ง and polarized along ๐‘ฅห† axis. Furthermore, we compare radar cross sections (RCS) in the ๐œ™ = 0 plane, using the proposed methods against either analytical data or a validated method of moments code that is based on RWG basis functions, otherwise referred to as RWG-CFIE. For every scattering experiment presented in the tables, the maximum relative far-field error, denoted by ๐œ– โˆž , is defined as ๐‘๐‘Ž๐‘™๐‘ ๐‘Ÿ๐‘’ ๐‘“ max |Eโˆž (xฬ‚) โˆ’ Eโˆž (xฬ‚)| xฬ‚ ๐œ–โˆž = ๐‘Ÿ๐‘’ ๐‘“ , (5.20) max |Eโˆž (xฬ‚)| xฬ‚ ๐‘Ÿ๐‘’ ๐‘“ where the reference solutions Eโˆž was computed by Mie series in the case of spherical scatterers, otherwise, by a Loop subdivision based CC-CFIER. All of the numerical results presented in the tables and graphs in this section were obtained by prescribing a GMRES residual tolerance equal to 10โˆ’4 for the overall system and 10โˆ’10 for inverting the gram matrix with a diagonal preconditioner. Finally, we note that for electrically larger scatterers, we provide the iteration count to reach the specified GMRES tolerance, the time taken to reach the prescribed tolerance, and the error relative to the benchmark data. 5.4.1 Accuracy of Regularized Combined Field Integral Equation In the first set of numerical results, we aim to compare the accuracy and high order nature of the proposed approaches for the analysis of EM scattering against an analytical solution, as well as the number of iterations required by the GMRES solver to reach the prescribed tolerance. To this end, we consider a sphere of diameter 8๐œ† that is modeled using an initial control mesh comprising of 642 vertices and 1280 faces. We consider two meshes generated by refining the initial control mesh once and thrice, respectively, using Loop subdivision. Note, unlike typical mesh refinement, under the rules of subdivision, the limit surface that all meshes point to is identical. More to the point, all the required numerics are carried out on the limit surface, NOT the Lagrangian geometric approximation. This refinement process leads to a coarser sphere of 2,562 vertices and 5,120 faces and a finer 76 one composed of 40,482 vertices and 80,960 faces. The main benefits in refining a mesh is better approximation of the physics on the limit surface. In the experiments discussed next, the finer discretization was used with RWG basis (together with a Lagrangian geometry description). We ensured that the surface areas of the Lagrangian mesh agree within 99% to the subdivision mesh. In Figure 5.3, we compare RCS data on an 8๐œ† sphere. The degrees of freedom are as follows: for CC-CFIER: MH we use 1200 MHs, leading to 2400 DoF; RWG-CFIE results in 122,880 DoF, and CC-CFIER: Loop contains 5124 DoF. As is evident from Fig. 5.3, the agreement between the three sets of numerical data to analytical solutions is excellent. In addition, we have analyzed a series of electrically larger spheres. These geometries are obtained via refinement of the initial mesh, such that at any frequency, the edge length is approximately 0.3๐œ†. The details of these experiments are presented in Table. 5.1. As is evident from this table, there is excellent agreement between the proposed methods and analytic data. The convergence of Loop and MH implementations of CC-CFIER is approximately the same as is the total solve time. The approximately four fold compression is not sufficient to affect the overall solve time due to the well-conditioned gram matrix for the sphere. 77 Figure 5.3: Radar cross section of the sphere (๐œ™ = 0 cut). CC-CFIER-Loop CC-CFIER-MH Size ๐‘๐ฟ /๐‘๐‘€๐ป It./Total Time ๐œ–โˆž It./Total Time ๐œ–โˆž 8๐œ† 5124/2000 7/0m 35s 5.99E-4 7/0m 33s 6.32E-4 16๐œ† 20484/6000 8/4m 31s 5.99E-4 8/4m 26s 9.29E-4 32๐œ† 81924/24000 9/25m 42s 2.26E-4 9/25m 47s 2.33E-3 Table 5.1: Convergence data for a spheres of different diameters: 8๐œ†-32๐œ†. 5.4.2 EM Scattering from Complex Objects In this section, we provide several examples to demonstrate the viability of using the formulations presented here for EM scattering on complex objects. We do so by comparing our results obtained from CC-CFIER: MH against those obtained using the CC-CFIER: Loop and RWG-CFIE. First, we consider the bumpy cube shown in Fig. 5.4, that fits in a 8๐œ† ร— 8๐œ† ร— 8๐œ† box. The number of DoFs for the RWG-CFIE is 122880, whereas for the CC-CFIER: Loop it is 5124 and 78 2400 for the CC-CFIER: MH formulation; Fig. 5.4 illustrates excellent agreement between all three. As before, we use mesh refinement to generate electrically larger structures. The results of these runs are presented in Table. 5.2, specifically, iteration count for CC-CFIER: Loop and CC-CFIER: MH formulation. Figure 5.4: Radar cross section of the bumpy cube (๐œ™ = 0 cut). CC-CFIER-Loop CC-CFIER-MH Size ๐‘๐ฟ /๐‘๐‘€๐ป It./Total Time It./Total Time ๐œ–โˆž 8๐œ† 5124/2400 11/1m 4s 11/1m 2s 2.13E-3 16๐œ† 20484/7200 12/7m 0s 12/6m 37s 5.13E-3 32๐œ† 81924/28000 12/38m 29s 13/36m 10s 4.44E-3 Table 5.2: Convergence data for a bumpy cube of sizes varying from 8๐œ† โˆ’ 32๐œ†. We report that the iteration count is low, approximately the same for both Loop and MH, and both took approximately the same time for the matrix solve. The agreement between Loop and the compressed MH system is also excellent. 79 Figure 5.5: Radar cross section of the shuttle (๐œ™ = 0 cut). CC-CFIER-Loop CC-CFIER-MH Size ๐‘๐ฟ /๐‘๐‘€๐ป It./Total Time It./Total Time ๐œ–โˆž 20๐œ† 7942/4000 78/11m 24s 39/5m 15s 2.43E-3 40๐œ† 31684/12000 38/30m 57s 29/19m 46s 2.00E-3 80๐œ† 126724/36000 28/187m 30s 29/105m 47s 2.70E-3 Table 5.3: Data for shuttle geometries from 20๐œ† โˆ’ 80๐œ†. Next, we consider a shuttle that that fits in a 20๐œ† ร— 12.22๐œ† ร— 7.22๐œ† box. The number of DoFs for the RWG-CFIE is 190080 whereas for the CC-CFIER: Loop is 31684 and for CC-CFIER: MH we have 6000. From Fig. 5.5 shows excellent agreement between all three. Again, we refine the geometry to consider electrically larger scatterers, in this case up to 80๐œ†. Table. 5.3 reports the iteration count, for CC-CFIER: Loop and CC-CFIER: MH basis, as we increase the frequency. We find that the iteration count is stable for both formulation, and they are in excellent agreement. Further, we note the significant compression achieved 80 via MHBs. Figure 5.6: Radar cross section of the jet airliner (๐œƒ = 90 cut). CC-CFIER-Loop CC-CFIER-MH Size ๐‘๐ฟ /๐‘๐‘€๐ป It./Total Time It./Total Time ๐œ–โˆž 30๐œ† 12132/7000 57/14m 59s 39/9m 45s 3.97E-3 60๐œ† 48516/21000 42/82m 12s 39/52m 10s 4.78E-3 120๐œ† 194052/63000 41/739m 46s 40/275m 28s 1.34E-2 Table 5.4: Data for jetliner geometries from 30๐œ† โˆ’ 120๐œ†. Finally, we consider a Jet airliner that fits in a 18๐œ† ร— 17๐œ† ร— 5๐œ† box. In this example, the plane wave propagating in the ๐‘ฆห† direction (incident on the nose) and polarized along ๐‘ฅห† direction. The number of DoFs for the RWG-CFIE is 72768 whereas for the CC-CFIER:Loop is 12132 and the CC-CFIER: MH is 5000. It is evident from Fig. 5.6 that all three data sets agree well with each other. In Table. 5.4, we report the iteration count, for CC-CFIER: Loop and CC-CFIER: MH basis, as we increase the electrical size of the object. We find that the 81 iteration count is stable for both formulation, as well as excellent agreement. Also, note the excellent compression produced by MHBs. 5.4.3 Scattering from Multiply connected Complex Objects In what follows, we analyze the performance of two different cases of the proposed approach for studying the radar cross section (RCS) and compare it with a conventional MoM CFIE solution technique that relies on RWG basis functions, otherwise referred to as RWG-CFIE. These two cases are (a) the prescribed formulation in this paper, which we refer to as CC-CFIER: Loop + GL and (b) the same prescribed formulation, but we omit the harmonic component of the surface current, which we refer to as CC-CFIER: Loop; in other words we set J๐‘–๐‘˜ = 0, for ๐‘˜ โˆˆ {3, 4} and ๐‘– โˆˆ {1, . . . , ๐‘”}. In particular, we compare all three generated RCSs against each other, we report the iteration count to reach the specified GMRES solver tolerance and the time taken to reach the prescribed tolerance. Unless otherwise stated, we are comparing RCS in the ๐œ™ = 0 plane, due to a plane wave field propagating in ๐œ…ห† = โˆ’ห†๐‘ง and polarized along ๐‘ฅห† axis. Furthermore, in the experiments discussed next, the finer discretization was used with RWG basis (together with a Lagrangian geometry description). We ensured that the surface areas of the Lagrangian mesh agree within 99% to the subdivision limit surface. All of the numerical results presented in the graphs in this section were obtained by prescribing a GMRES residual tolerance equal to 10โˆ’7 for the overall system and 10โˆ’11 for inverting the gram matrix. The first example we consider is a double torus of genus 2, see Fig. 5.7a, that fits in a 10๐œ† ร— 4.81๐œ† ร— 2.06๐œ† box. The number of DoF for the CFIE: RWG is 36864, CC-CFIE: Loop + GL results in 6146 DoF, and lastly the CC-CFIE: Loop has 6142 DoF, respectively. The CFIE: RWG converges in 77 iterations after 2 m 7 s, CC-CFIE: Loop + GL converges in 17 iterations after 1 m 8 s, and CC-CFIE: Loop reaches tolerance within 22 iterations, for a total of 1 m 7 s. From Fig. 5.8, we report excellent agreement between both the CC-CFIE: Loop + GL and the RWG-CFIE and we have highlighted the glaring disparity when omitting the 82 (a) (b) Figure 5.7: Comparison of interactions between FMM and wideband MLFMA. (a) Double torus with 2 global loops. (b) Chmutov surface with 28 global loops. Figure 5.8: Radar cross section of the double torus (๐œ™ = 0 cut). harmonic component in CC-CFIE: Loop between the other methods in the inset. The second example is a Chmutov surface, see Fig. 5.7b, that fits in a 5๐œ† box with a genus of 28. The number of DoF for the CFIE: RWG is 144000, CC-CFIE: Loop + GL has 23950 DoF, and lastly the CC-CFIE: Loop results in 23894 DoF, respectively. The CFIE: RWG converges after 350 iterations in 26 m 25s, CC-CFIE: Loop + GL converges after 36 83 iterations in 3 m 48 s, and lastly, the CC-CFIE: Loop converges after 104 iterations in 10 m 20 s. Fig. 5.9 shows excellent agreement between both the CC-CFIE: Loop + GL and the RWG-CFIE and lastly, we highlight the glaring disparity when omitting the harmonic component in CC-CFIE: Loop between the other methods in the inset. Figure 5.9: Radar cross section of the Chmutov (๐œ™ = 0 cut). 5.4.4 Scattering analysis with Localized Manifold Harmonics (a) (b) (c) Figure 5.10: Localized manifold harmonics for the red region shown on the left. We show a couple localized manifold harmonics. Next, we demonstrate the utility of the combination of MHB and L-MHB to an electromagnetic scattering problem. The currents can be defined in a manner identical 84 to (3.14) with both MHB and L-MHB used as the basis set in (5.8a) [6]. In this example, we consider scattering from a cone sphere, see Fig. 5.10, that fits in a 24๐œ† ร— 6.78๐œ† ร— 6.78๐œ† box. We compare three sets of far-field data, those obtained from (i) CC-CFIER: Loop, (ii) CC-CFIER: MHB and (iii) CC-CFIER: MHB + L-MHB. The aim of these comparisons is to demonstrate the efficiency and accuracy of the MHB + L-MHB basis set relative to MHB and Loop subdivision. The region near the tip is enriched with L-MHB as shown in Fig. 5.10. The rationale for doing so exploits the knowledge that near tips currents tend to be localized and are smoother away; hence, the need for enrichment via singular basis or dense discretization near tips [116, 52, 70]. Figure 5.11: Radar cross section of the cone sphere (๐œ™ = 0 cut). 85 Figure 5.12: Radar cross section of the tetrahedron (๐œ™ = 90 cut). (a) (b) Figure 5.13: The eigenvalues at different levels for the tetrahedron (Fig. 3.10), note the similar trends. The ratio ฮ›๐‘™+1 /ฮ›๐‘™ for the different levels, note that the similarity breaks later for higher levels, and the ratio is very close to 1. The reference solution, CC-CFIER: Loop, has 24198 DoF and the solution converges in 302 iterations for 7.7 seconds per iteration. For comparison, we begin with 3000 MHB 86 in CC-CFIER: MHB 1, which leads to 6000 DoF, converges in 17 iterations for 7.9 seconds per iteration; the maximum relative far-field error, ๐œ– โˆž , here is 6.47E-4. For the CC-CFIER: MHB + L-MHB there are 2474 MHB and 526 L-MHB for a total 6000 DoF. Overall, the system reaches tolerance after 29 iterations for 7.8 seconds per iteration. For this case, the maximum relative far-field error, ๐œ– โˆž , is 9.69E-5. We find that while the convergence rate suffers slightly for MHB + L-MHB, the system is significantly compressed and far more accurate in the given metric. In order to achieve the same level of accuracy between the two methods, we increase the size of the MHB basis to 10000, resulting in 20000 DoF; we report this data as CC-CFIER: MHB 2. The solution converges in 18 iterations for 8.7 second per iteration. Lastly, we plot the RCS generated from all methods discussed above in Fig. 5.11. We see excellent agreement between CC-CFIER: Loop, CC-CFIER: MHB 2, and CC-CFIER: MHB + L-MHB. The main discrepancy between all the methods is highlighted in the inset; it is apparent that the L-MHB resolves the information at the tip far better than the stand alone MHB. 5.4.5 Scattering from a Multi-resolution Manifold Next, we provide details on using this framework for electromagnetic analysis of scattering from a simply connected object. Consider a tetrahedron that fits in a 8๐œ† ร— 8๐œ† ร— 8๐œ† box, as shown in Fig. 3.10. The surface ฮ“ is parametrized by a series of three sequentially refined control meshes. We follow a hierarchical multiplicative projection procedure [138] in Alg. 2 to solve (5.8a). We compare three sets of results: CC-CFIER: Loop, CC-CFIER: MHB and CC-CFIER: H-MHB. For H-MHB, the three levels of refinement are used starting with 514 DoF Fig. 3.10a, 2050 DoF Fig. 3.10a, and lastly 8194 DoF Fig. 3.10a. We defined H-MHB over all 3 meshes starting with the first 100 MHB on Fig. 3.10a, the next 900 on Fig. 3.10b, and the last 2000 on the Fig. 3.10a for a total of 3000 MHB. The eigenspectrum of this basis is illustrated in Fig 5.13. The remaining proposed methods are solved on the finiest control mesh, with CC-CFIER: Loop having 16390 DoF and CC-CFIER: MHB having 6000 DoF. 87 Algorithm 3: Hierarchical Multiplicative Projection Procedure Algorithm. input : A series of Loop subdivision meshes ๐‘‡ ๐‘™ , ๐‘™ โˆˆ {0, . . . , ๐‘“ } H-MHB ๐ป ๐‘™ , ๐‘™ โˆˆ {0, . . . , ๐‘“ } Control error ๐œ– output : Surface current J defined on ๐‘‡ ๐‘“ 1 //Initialization 2 while ๐œ– โˆž โ‰ค ๐œ– do 3 for ๐‘– := 1 to ๐‘“ โˆ’ 1 do 4 w = [๐‘… ๐‘– ] ยท ยท ยท [๐‘… ๐‘“ โˆ’1 ][๐บ ๐‘“ ]โˆ’1 ([๐‘‰ ๐‘“ ] โˆ’ [๐‘ ๐‘“ ][๐ผ ๐‘“ ]) 5 Solve [๐‘ ๐‘– ]y = [๐ป ๐‘– ]๐‘‡ [๐บ ๐‘– ]w 6 [๐ผ ๐‘“ ] = [๐ผ ๐‘“ ] + ๐›ผ[๐‘ƒ ๐‘“ โˆ’1 ] ยท ยท ยท [๐‘ƒ ๐‘– ][๐ป ๐‘– ]y 7 end 8 Solve [๐‘ ๐‘“ ]y = [๐ป ๐‘“ ]๐‘‡ [๐บ ๐‘“ ]([๐‘‰ ๐‘“ ] โˆ’ [๐‘ ๐‘“ ][๐ผ ๐‘“ ]) 9 [๐ผ ๐‘“ ] = [๐ผ ๐‘“ ] + ๐›ผ[๐ป ๐‘“ ]y 10 compute ๐œ– โˆž 11 end We plot the RCS generated from all methods discussed above in Fig. 5.12. As is evident, there is excellent agreement between CC-CFIER: Loop, CC-CFIER: MH, and CC-CFIER: H-MHB. The reference solution, CC-CFIER: Loop, reaches the prescribed GMRES tolerance after 82 iterations. The maximum relative far-field error CC-CFIER: H-MHB is 6.93E-4 converging after 3 iterations; GMRES tolerance for the first band is reached after 3 iterations, 12 iterations for the second band, and 7 iterations for the last. The CC-CFIER: MHB has a maximum relative far-field error of 7.40E-5 and reaches GMRES tolerance after 13 iterations. As is evident from this experiment H-MHB permits reusability of data from coarser meshes. 5.5 Conclusion In this chapter, we have presented an integral equation based isogeometric analysis method for electromagnetic scattering from subdivision surfaces; in presenting this approach, we assumed a multiply connected structure, used a complete surface Helmholtz decomposition to effect a Calderรณn operator. In particular, this work extends the previous numerical infrastructure on to existing Loop subdivision IGA methods to enable analysis of 88 multiply connected structures with the help of random vector projections onto subdivision surfaces. As a result, we are able to prescribe a complete Helmholtz decomposition on the surface wherein orthogonality of different components is preserved. Furthermore, we address the main bottleneck of evaluating inner-products for higher order basis functions on higher order surfaces by employing a wideband MLFMA to evaluate โ€œallโ€ interactions. Finally, we introduce the notion of MHB as a means to represent the currents on the surface. These geometry basis can be used for compression of both the manifold and physics on the manifold. In addition, we have explored the benefits of MHB; specifically, we have demonstrated that (a) it is possible to use MHB for multiply connected objects, (b) it is possible to locally enrich MHB leading to better capture of localized behavior, and (c) using a composition of data from a hierarchy of refined geometries to compute results from the finest geometry. All aforementioned attributes have been demonstrated in numerous results using both the subdivision and MH basis, on a collection of electrically large geometries. This chapter, ยฉ2022 IEEE, in part is reprinted, with permission, from A. M. A. Alsnayyan and B. Shanker, "Iso-geometric Integral Equation Solvers and their Compression via Manifold Harmonics," IEEE Transactions on Antennas and Propagation on, April 2022, with slight modifications to fit in this dissertation. This chapter, ยฉ2022 IEEE, in part, is reprinted, with permission from A. M. A. Alsnayyan and B. Shanker, "A complete Helmholtz decomposition on multiply connected subdivision surfaces and its application to integral equations," IEEE Transactions on Antennas and Propagation on, August 2022, with slight modifications to fit in this dissertation. This chapter, ยฉ2022 IEEE, in part, is reprinted, with permission from A. M. A. Al- snayyan and B. Shanker, "Manifold Harmonics and their application to computational electromagnetics," IEEE Journal on Multiscale and Multiphysics Computational Techniques on, August 2022, with slight modifications to fit in this dissertation. 89 CHAPTER 6 SHAPE OPTIMIZATION 6.1 Introduction Research into inverse scattering dates back decades and has found applications in a number of wide ranging fields, including areas such as medical diagnostics, detection of buried objects, tomography, and non-destructive evaluation [42, 154, 166, 115, 64, 94]; in these problems, the goal is to retrieve the distribution of constitutive properties in a domain and/or geometry given a set of measured scattered field data. Along these lines of an inverse scattering problem, there are two subclasses of problems that can be considered: (a) can one modify shapes such that one obtains desired scattered fields, or (b) reconstruct the shape of a object given scattered field data (and boundary conditions on the surface). As is evident, both problems are closely related. The former is commonly referred to as shape optimization; it has a number of applications ranging from acoustics [21, 21, 102, 87, 81] to electromagnetics [141, 28] to medical imagining [166, 115] to design of horns [21], and a number of other applications [167, 61]. The computational techniques used therein have been integrated with local based optimization methods [19, 28, 81] as well as global optimization methods [129, 66]. This chapter is devoted to the latter problem - the recovery of the shape of scattering obstacles from phaseless far-field data. Interest in this class of problems is widespread as it finds application in a number of different disciplines [47]. The main challenges that arise in this problem are (a) the ill-posed nature of the problem, (b) the number of degrees of freedom in the optimization problem, and (c) the optimal minimization method used for reconstruction. In this chapter, the aim is to address challenge (b) as it arises in a shape reconstruction routine that is constrained by the forward scattering problem, with the goal of minimizing the geometrical parameters that define the boundary of the scatterer, or in other words minimizing the search space. In particular, we look at 90 constructing and then modifying a compressed parameterization of some starting surface until we reach a prescribed minimization of a cost functional that measures the discrepancy between the prescribed far-field data and the far-field pattern corresponding to the current approximation of the paramterized scatterer. As is apparent, the geometry parameterisation and its interplay with the analysis plays a crucial role in shape optimisation/reconstruction. When a geometric representation, distinct from the computer-aided design (CAD) model is used to represent the shape, additional errors are incurred as well as a large parameterisation space is needed to represent high fidelity shapes, drastically increasing the difficulty of the optimization scheme. In addition, it leads to non-physical oscillations in the optimised geometry [75] as well as severely distorted mesh requiring auxiliary mesh smoothing, further complicating the problem. To remedy these difficulties, the shape optimisation/reconstruction can be constructed in a high-order isogeometric framework. The effects of geometry representation and analysis on piecewise constant meshes are exacerbated when involved in shape reconstruction. Our approach. Specifically, with flat tesselation, both the number of degrees of freedom necessary to represent (and perturb) the geometry is significant as is the number of basis functions necessary to model the physics. In what follows, our approach is to use a higher order smooth representation of our initial guess as well higher order basis within a framework of isogeometric methods. Elsewhere [6], we have shown that these methods are effective and agree well with those computed using piecewise constant tesslations and basis (with the latter at higher cost). Isogeometric methods use the same underlying basis sets to represent both the geometry and physics on the geometry. This class of techniques was pioneered by Hughes [79], and has since been adapted for a number of different types of problems in structural mechanics [45, 102], electromagnetics [98, 100, 97, 67, 5, 6] and acoustics [6, 105, 39]. The research in using isogeometric analysis has treaded along two paths: use of non-unifrom B-splines and subdivision surfaces; our focus is on subdivision. Subdivision is a powerful geometric modelling technique for generating smooth surfaces on arbitrary connectivity meshes 91 which are the generalization of splines to arbitrary connectivity meshes. Specifically, we use the Loop scheme based on triangular meshes and quartic box-splines [106]. Subdivision surfaces are perhaps the ideal candidate for IGA shape reconstruction as one can exploit other facets of subdivisionโ€“hierarchical refinement, higher order continuity, and arbitrary topology, i.e., no other restricting assumptions on the geometry need be considered (star- shaped domains, for example), see [160, 45, 27] and references therein for examples on IGA shape optimization. The general approach in these papers is to exploit the hierarchy of the control mesh underlying a subdivision surface to do multiresolution editing [108, 145, 19]. The coarse control mesh vertex positions are modified to perform large-scale editing and the fine control mesh vertex positions are modified to add localized changes. This allows for a relatively compressed representation of geometry, allowing for a smaller search space and thus computationally feasible shape reconstruction problem. This approach works well. But, one of the questions that we ask is whether we can construct a better compression scheme. An inspiration for the answer to this question can be found in computer graphics [134, 157, 122, 96]. The nub of these ideas is to develop a mechanism that enables compression and morphing of manifolds in an efficient manner. These techniques rely on the manifold harmonic basis (MHB) set; a basis set constructed from the eigenfunctions of the Laplace-Beltrami operator (LBO). As is well known, LBO is a self-adjoint linear surface operator that captures all the intrinsic properties of the shape and is invariant to extrinsic shape transformations such as isometric deformations. Its eigenfunctions or the MHBs can be understood as a generalization of the Fourier spectrum for functions defined on a general surface manifold. They provide a unique, compact, elegant, and multi-resolution basis for spectral shape processing that is independent of the actual shape representation and parameterization. Several successful applications have been proposed that take advantage of these desirable properties, such as spectral geometry filtering, compression, and surface deformation [157, 29]. As subdivision basis are higher order and smooth it (a) enables further compression (of both the surface and 92 physics) and (b) has better convergence rates when used for analysis [121, 98]. Together, these features result in fewer degrees of freedom in an optimization setting. Indeed, while manifold harmonics can be defined on piecewise tesselation, additional geometric degrees of freedom required and the number of manifold harmonics required (due to lack of smoothness) make it a poor choice. Finally, we will leverage our earlier work on subdivision-based isogeometric methods for acoustics to develop a framework that carries over the benefits of subdivision meanwhile maintaining the favorable properties of MHs for use in both analysis and morphing the manifold for shape reconstruction. The specific contributions of this chapter are as follow: 1. develop a MHB based compression scheme for representation of and analysis on manifolds, 2. develop an inverse source based initial estimate for shape, 3. develop a multi-resolution framework for shape reconstruction, 4. and demonstrate the viability of this technique on a number of challenging targets. Note, the synthetic data is obtained using an isogeometric method based on subdivision surfaces. We have also shown elsewhere (for both acoustic and electromagnetic scatter- ing) that the difference between data obtained using a subdivision based isogeometric method and basis defined on piecewise triangulation is not significant but the latter comes at significantly higher cost. The methods developed herein are agnostic to the specific optimization technique used in the reconstruction process. In other words, the novelty/contributions of this work do not lie in the optimization scheme chosen, but the multi-resolution framework used to systematically update and refine the geometry and physics throughout the reconstruction process. And as such, we have used readily available libraries for the optimization procedure (Method of Moving Asymptotes (MMA) [82]) and use the straight-forward finite difference method to implement the necessary gradient. We note that there exists the more efficient adjoint based methods that can be implemented [93, 90, 128, 166, 95]. 93 The remainder of this chapter is organized as follows: In Section 6.2, the forward scattering constraint and the shape reconstruction problem are posed. In Section 6.3, we propose a backprojection scheme that initialize our shape reconstruction algorithm. In Section 6.4, we present a number of results on structurally challenging objects to demonstrate the salient features of our reconstruction algorithm. Finally, in Section 6.5, we summarize our contribution in this body of work. 6.2 The forward scattering and shape reconstruction problems Consider a soft scatterer embedded in ฮฉ โˆˆ R3 , whose boundary is denoted by ฮ“0 with a uniquely defined outward pointing normal nฬ‚, as depicted in Fig. 6.1. Figure 6.1: The shape reconstruction problem. The objective of this experiment is to determine the shape of the scatter ฮ“0 such that the measured scattered fields ฮฆ๐‘š ๐‘  (r) at ฮ“ , obtained for a collection of directions and ๐‘  wavenumber, d๐‘– , ๐œ… ๐‘— satisfy ฮฆ๐‘š ๐‘  (r) = ฮฆ๐‘  (r). This is equivalent to minimizing the following  ๐บ objective function โˆซ 1 ๐‘  ๐‘  find ๐‘Ž๐‘Ÿ ๐‘”๐‘š๐‘–๐‘› ๐ฝ(ฮ“) := |ฮฆ๐บ (r) โˆ’ ฮฆ๐‘š (r)| 2 ๐‘‘r. (6.1) ฮ“โˆˆR3 2 ฮ“๐‘  As stated, this is a constrained shape reconstruction problem with ฮ“ as the design variable and the forward far-field scattering problem as the constraint on the admissible exterior phaseless far-field data ฮฆ๐‘š ๐‘  . In this case the objective function is defined in the 94 limiting case in which ฮ“๐‘  goes to infinity; following standard practice, fields are normalized by distance. ๐‘  On ฮ“๐‘  , data due to fields scattered by the obstacle, ฮฆ๐บ (r)(d๐‘– ,๐œ… ๐‘— ) , is available for a combination of incident waves propagating in direction d๐‘– โˆˆ D, where D = {d1 , d2 , ยท ยท ยท , d๐‘› }, with wavenumbers ๐œ… ๐‘— โˆˆ K, where K = {๐œ…1 , ๐œ…2 , ยท ยท ยท , ๐œ… ๐‘˜ }, denoted as ฮฆ๐‘– (r)(d๐‘– ,๐œ… ๐‘— ) (in the what follows, the subscripts will be omitted for notational simplicity). For each incident field, the resulting total field ฮฆ๐‘ก (r) = ฮฆ๐‘  (r) + ฮฆ๐‘– (r) satisfies the following boundary value problem โˆ‡ฮฆ๐‘ก (r) + ๐‘˜ 2 ฮฆ๐‘ก (r) = 0 r โˆˆ ฮฉ, (6.2a) ฮฆ๐‘ก (r) = 0 r โˆˆ ฮ“0 , (6.2b) โˆš ๐œ•ฮฆ๐‘    lim ๐‘Ÿ โˆ’ ๐‘–๐œ…ฮฆ๐‘  = 0 r โˆˆ ฮฉ. (6.2c) ๐‘Ÿโ†’โˆž ๐œ•๐‘› Here, ฮฆ๐‘  denotes the scattered field, (6.2b) defines the sound-soft boundary condition and (6.2c) imposes the Sommerfeld radiation condition. 6.2.1 Boundary Element Formulation In this section, we will provide a general outline on applying boundary element method to solve (6.2) and obtain the measured far-field scattered fields ฮฆ๐‘š ๐‘  (r), where r โˆˆ ฮ“ , in ๐‘  order to compute the objective function and its gradient. In particular, the Burton-Miller formulation [36] provides a unique solution, and is expressed as: โ„’ ๐ต๐‘€ [ฮ›, ฮ“] (r) = ๐‘‰ ๐‘– (r) r โˆˆ ฮ“, โ„’ ๐ต๐‘€ [ฮ›, ฮ“] (r)  (1 โˆ’ ๐›ผ)๐’ฎ[ฮ›, ฮ“](r) + ๐›ผ๐›ฝ๐’Ÿ 0[ฮ›, ฮ“](r) r โˆˆ ฮ“, (6.3) ๐‘‰ ๐‘– (r)  (1 โˆ’ ๐›ผ)ฮฆ๐‘– (r) + ๐›ผ๐›ฝ nฬ‚ ยท โˆ‡ฮฆ๐‘– (r) r โˆˆ ฮ“. where ฮ›(r) denotes ๐œ•n ฮฆ(r), nฬ‚ is the outward unit normal vector to ฮ“, 0 โ‰ค ๐›ผ โ‰ค 1 is a coupling factor that makes the solution unique at all frequencies [36], and ๐›ฝ is the constant weighting factor. Note, setting ๐›ผ = 0 or ๐›ผ = 1 introduces a non-trivial null space at frequencies that correspond to the interior resonance of the structure; the reader is referred to [150] for a theoretical explanation. Here, ๐’ฎ : ๐ป โˆ’1/2 (ฮ“) โˆ’โ†’ ๐ป 1/2 (ฮ“) denotes the single layer boundary 95 integral operator โˆซ ๐’ฎ[ฮ›, ฮ“](r)  ฮ›(r0)๐บ(r, r0)๐‘‘r0, r โˆˆ ฮ“, (6.4a) ฮ“ and ๐’Ÿ 0 : ๐ป โˆ’1/2 (ฮ“) โˆ’โ†’ ๐ป โˆ’1/2 (ฮ“) denotes the adjoint double layer operator ๐œ•๐บ(r, r0) 0 โˆซ  ๐’Ÿ 0[ฮ›, ฮ“](r)  ฮ›(r0) ๐‘‘r , r โˆˆ ฮ“. (6.4b) ฮ“ ๐œ•n 0 where ๐บ(r, r0) = ๐‘’ โˆ’๐‘–๐œ…|rโˆ’r | /4๐œ‹|r โˆ’ r0 | is the free-space Helmholtz kernel in R3 , ๐œ… is the wavenumber, and ๐›ฝ = ๐‘–/๐œ…. An ๐‘’ ๐‘–๐œ”๐‘ก dependence is assumed and suppressed. We introduce the far-field operator โ„’ ๐‘“ ๐‘Ž๐‘Ÿ : ๐ป โˆ’1/2 (ฮ“) โˆ’โ†’ ๐ฟ2 (ฮ“๐‘  ) given by ๐‘  (6.5a) ฮฆ๐‘š (rฬ‚)  โ„’ ๐‘“ ๐‘Ž๐‘Ÿ [ฮ›, ฮ“](rฬ‚), โˆซ 1 ฮ›(r0)๐‘’ ๐‘–๐œ…rฬ‚ยทr ๐‘‘r0, rฬ‚ โˆˆ ฮ“๐‘  . 0 โ„’ ๐‘“ ๐‘Ž๐‘Ÿ [ฮ›, ฮ“](rฬ‚)  (6.5b) 4๐œ‹ ฮ“ The exact solution of (6.3) is generally not available. The numerical solution of the integral equations is effected in a discrete setting using an isogeometric method based on Loop subdivision; the exact details are provided in Chapter 4. 6.2.2 The objective function The strategy we use to minimize the cost functional (6.1) is a gradient-based approach which requires the derivative of the cost functional (6.1) with respect to the domain Figure 6.2: The perturbed shape according to design perturbation h(r). 96 perturbations. The domain perturbation, defined by ฮ“ โ†’ ฮ“ โ„Ž๐œ is described by a vector displacement field h(r) such that r โ†’ r + ๐œh(r), where r โˆˆ ฮ“ and ๐œ is a scalar parameter that denotes the amount of shape change, as illustrated in Fig. 6.2. The reader is referred to [16, 84], and references therein for a more detailed definition of shape perturbation as we have only outlined a brief sketch here. To this end, we define the gradient of the cost functional (6.1) as ๐ฝ(ฮ“ โ„Ž๐œ ) โˆ’ ๐ฝ(ฮ“) ๐ฝ 0(ฮ“) = lim . (6.6) ๐œโ†’0 โ„Ž Where ๐ฝ(ฮ“) and ๐ฝ(ฮ“ โ„Ž๐œ ) are the cost functional evaluated for the reference and perturbed domain. A critical step of this shape reconstruction algorithm is the evaluation of both the objective function and itโ€™s gradient which are required for minimizing the objective function. In this case, the objective function is dependent on the shape of the boundary and the perturbation vector field h is represented using shape coefficients ๐›ฝ = ๐›ฝ ๐‘– , wherein ฮ“(๐›ฝ)  described by a linear combination of manifold harmonics (see chapter 3). Therefore, we can now repose the reconstruction problem as โˆซ 1 ๐‘  ๐‘  find ๐›ฝ such that ๐ฝ(ฮ“(๐›ฝ)) := ๐‘š๐‘–๐‘› |ฮฆ๐บ (r) โˆ’ ฮฆ๐‘š (rฬ‚)| 2 ๐‘‘r. (6.7) 2 ฮ“๐‘  where ฮฆ๐‘š ๐‘  (rฬ‚) is the modulus of the far-field scattered data for some object ฮ“, described by a set of shape coefficients ๐›ฝ. The gradient of the objective function ๐ฝ with respect to the geometrical parameter ๐›ฝ ๐‘– can be done using a finite-difference method [159] ๐œ•๐ฝ(๐›ฝ) ๐ฝ(๐›ฝ 1 , ยท ยท ยท , ๐›ฝ ๐‘– + ๐œ๐‘– , ยท ยท ยท , ๐›ฝ ๐‘๐‘ฃ ) โˆ’ ๐ฝ(๐›ฝ1 , ยท ยท ยท , ๐›ฝ ๐‘– , ยท ยท ยท , ๐›ฝ ๐‘๐‘ฃ ) = . (6.8) ๐œ•๐›ฝ ๐‘– ๐œ๐‘– This finite-difference scheme is the most straightforward implementation, given that it still enables us to highlight the main objective of this implementation: the use of multiresolution MHB in shape optimization. The reader is turned to [166, 159] for more efficient approaches that can be adopted to this scheme. The evaluation of ฮฆ๐‘š ๐‘  in effect, the objective function ๐ฝ, for some shape is done by discretizing our forward problem (6.2) by the boundary element method. 97 6.2.3 Multi-resolution shape reconstruction In this application, we aim to utilize the natural multi-resolution framework of MHB detailed above for the shape reconstruction problem. Specifically, we denote the set of coefficients in (6.9) using ๐›ฝ = ๐›ฝ 1 , ยท ยท ยท , ๐›ฝ ๐‘๐‘ฃ . It is apparent that one can group    these in a collection of spatial frequency bands. In effect, ๐›ฝ = [โ„ฌ1 , ยท ยท ยท , โ„ฌ๐พ ] where  โ„ฌ๐‘– = [๐›ฝ ๐‘–,1 , ยท ยท ยท , ๐›ฝ ๐‘–,๐‘› ๐‘– ], where there are ๐‘› ๐‘– coefficients in each set. Note, the sets are a contiguous partition of ๐›ฝ . It follows, that one can implement a multi-resolution  analysis/optimization scheme by updating collection of โ„ฌ๐‘– one at a time, as opposed to the entire set. Thus, the magnitude of ๐›ฝ ๐‘– , dictate the relative importance of a given eigenfunction at a given spatial โ€œfrequencyโ€. This transformation makes it possible to define equivalent signal processing operations: this includes under/over sampling, filtering, multi-resolution, windowing, and so on. To wit, this representation makes it possible to enrich shapes (and functions defined on the shape) systematically. To illustrate, the effectiveness of this approach, consider the manifold ฮ“(r) again. Equation (6.9) can be rewritten as ร• ๐‘€ ๐‘๐‘ฃ ร• ฮ“(r) = ๐›ฝ ๐‘– ๐ป๐‘– (r) + ๐›ฝ ๐‘– ๐ป๐‘– (r) . (6.9) ๐‘–=1 ๐‘–=๐‘€+1 | {z } | {z } Low Frequency High Frequency In this expression, we have designated some cutoff coefficient ๐›ฝ ๐‘€ , or correspondingly, โˆš ๐œ† ๐‘€ to be a maximum cutoff โ€œfrequencyโ€. This effectively is a versatile tool for shape reconstruction that enables us to enrich data at different levels of fidelity. 6.3 Back-projection Based Initialization In cases of interest, a random guess to this shape reconstruction algorithm will unlikely return a satisfactory solution. As such, we need to obtain a good initial guess for ฮ“0 solely from the target scattered far-field data. In particular, we turn to a volume source reconstruction method to obtain a point cloud that can then be used to construct a subdivision mesh. Let ฮฆ๐‘  (๐œ… ๐‘› , ๐œ…ห† ๐‘š , rฬ‚) represent the target far-field data, where ๐œ… ๐‘› = ๐œ”๐‘› /๐‘ 98 is the wave number, ๐œ…ห† ๐‘š is an incident wave direction and rฬ‚ is the unit-vector along the direction of observation (alternatively defined in terms of ๐œƒ and ๐œ™). Following [10], the boundary of the object can be approximated by assuming there exists some source distribution ฮ›๐‘’ ๐‘ž (๐œ… ๐‘› , ๐œ…ห† ๐‘š , r) in a volume ๐‘‰ that satisfy the following minimization problem, ๐‘š๐‘–๐‘› k ฮฆ๐‘  (๐œ… ๐‘› , ๐œ…ห† ๐‘š , ๐‘Ÿห†) โˆ’ โ„’ ๐‘“ ๐‘Ž๐‘Ÿ ฮ›๐‘’ ๐‘ž (๐œ… ๐‘› , ๐œ…ห† ๐‘š , r) , ๐‘‰ k .   ๐‘›,๐‘š (6.10) Here, โ„’ ๐‘“ ๐‘Ž๐‘Ÿ ฮ›๐‘’ ๐‘ž , ๐‘‰ is a far-field projector similar to (6.5b); the domain of integration is   defined over ฮ“๐‘  and the function is defined over the volume ๐‘‰. We then coherently combine these volume equivalent sources in a manner following [111], ๐‘€ ๐‘–๐‘›๐‘ ๐‘ ๐‘“ ๐‘Ÿ๐‘’ ๐‘ž ฮฆ๐‘š,๐‘› ร• ร• ๐‘–๐œ… ๐‘› (rยทkฬ‚๐‘š ) ฮฆฬƒ๐‘’ ๐‘ž (r) = ๐‘’ ๐‘ž (r)๐‘’ , r โˆˆ ๐‘‰. (6.11) ๐‘š=1 ๐‘›=1 Now, we consider the cow in Fig. 6.3a as a candidate for constructing an initial Loop subdivision mesh using the volume source reconstruction method. We are given a set of target far-field scattering patterns obeying (6.2) due to 16 incident waves operating at a range of 200 Hz to 1200 Hz in 50 Hz increments. The reconstruction domain ๐‘‰ is a cubic region of 4 m x 4 m x 4 m discretized into 50 x 50 x 50 grid points. At each point in the reconstruction domain we evaluate the coherent summation of back propagated fields. The points that are of the highest intensity are closest to the boundary of the object ฮ“; this can be interrupted as an equivalent acoustics pressure in the volume. At this stage, we prescribe some threshold such that we isolate the points nearest to the boundary of the target object. This effect can be seen in Fig. 6.3b, wherein a cross-section of the reconstructed equivalent acoustic pressure along the ๐‘ง axis is plotted. The points that are left can be interrupted as a point cloud of the boundary, see Fig. 6.3d, which is used to construct the initial mesh, see Fig. 6.3c. 6.3.1 Shape Reconstruction Algorithm Thus far, we have developed/described a comprehensive set of tools that we employ for shape optimization. The reconstruction proceeds as follows: given initial farfield data, 99 (a) Original model. (b) Cross section of 3-D reconstructed object. (c) Point cloud. (d) Mesh generated from point cloud. Figure 6.3: Initial guess pipeline: VSRM. we (a) obtain an initial mesh, (b) construct MHs for this initial guess, i.e. the intial set of coefficients ๐›ฝ 0 , and (c) use a optimization routine to find the set of coefficients that minimizes the constrained optimization problem in (6.1) wherein (6.2) is the constraint on the admissible fields. The steps for our shape reconstruction algorithm are outlined in Alg. 4. In addition, we note that the optimization method used in this paper is MMA presented from the NLopt library [82]. However, the infrastructure presented is agnostic to the specific optimization methodology used. 6.4 Numerical Examples In this section, we present a series of examples that validate the approach presented in this body of work. Our goal is to reconstruct complex simply-connected Loop subdivision surfaces from synthetic data. As stated, we will assume that far field data is available for a set of incident illuminations and frequencies. An incident illumination, d๐‘– , is identified by 100 Algorithm 4: Shape reconstruction algorthim. 1 Define: 2 D: Directions of propagation 3 K: Wavenumbers 4 ฮฆ๐‘– : Incident fields for a set {D, K} 5 ฮ“: The boundary to be determined ๐‘  6 ฮฆ๐บ : Goal fields at ฮ“๐‘  due to incident fields ฮฆ๐‘  impinging on ฮ“ 7 Termination tolerance ๐œ– for optimization routine 8 9 Method: 10 VSRM: The Initial Guess ๐‘  11 Given ฮฆ๐บ , construct the volumetric acoustics source ฮฆ๐‘’ ๐‘ž,๐‘‰ using Eq. 6.11 12 e Generate an analysis-ready subdivision mesh โ„ณ0 , from the ฮฆ๐‘’ ๐‘ž,๐‘‰ point cloud using [?] 13 Solve Eq. 3.5 on โ„ณ0 to generate the MHs 14 Construct ๐›ฝ0 by performing a MHT on โ„ณ0 15 Optimization routine: 16 while ๐ฝ(๐›ฝ ๐‘– ) โ‰ฅ ๐œ– do 17 Solve the forward problem (Eq. 6.2) on โ„ณ ๐‘– 18 Compute objective function ๐ฝ(๐›ฝ) and corresponding sensitivity analysis Update design variables using an optimization algorithm: ๐›ฝ ๐‘– โ†’ ๐›ฝ ๐‘–+1   19 20 Update the subdivision mesh: โ„ณ ๐‘– โ†’ โ„ณ ๐‘–+1 21 end the pair (๐œƒ๐‘– , ๐œ™ ๐‘– ). The frequency of the incident field with wavenumber ๐œ… ๐‘— is denoted using ๐‘“ ๐‘— . In all cases presented, the solution outline is as follows: (a) we obtain synthetic data, or goal fields, from a subdivision description of the scatterer, (b) we add noise to this data defined by ๐‘  ๐‘  ๐‘  ฮฆฬ‚๐บ (r) = ฮฆ๐บ (r) + ๐›ฟฮฆ๐บ (r) (6.12) ๐‘  where ฮฆฬ‚๐บ (r) is the set of perturbed goal fields. The term ๐›ฟ > 0 is related to the signal-to-noise ratio (SNR) by the equation 10โˆ’๐‘†๐‘ ๐‘…/20 ๐›ฟ= โˆš , (6.13) ๐‘€ and is a unit vector with ๐‘€ random Gaussian entries. The ๐‘†๐‘ ๐‘… (in dB) typically fluctuates in real applications between 10 dBโ€“40 dB. For all our experiments, we choose a noise level of 10 dB ๐‘†๐‘ ๐‘… to demonstrate the robustness of our algorithm to noise; a ๐‘†๐‘ ๐‘… of 10 dB is 101 considered a high level of noise. Once we have our new perturbed goal fields, we employ our shape reconstruction routine following the steps presented in Fig. ??. Recall that at every iteration of this shape reconstruction algorithm we are physically perturbing the scattererโ€™s surface. Therefore, we must account for large deformations of the surface throughout the evolution of the algorithm. In order to do so, we require a high fidelity mesh. Given that we are using subdivision surface, we have โ‰ˆ ๐œ†/7 patches per wavelength. Furthermore, we note that in the case of a good initial guess, large mesh deformation will be far less prevalent and therefore easing the restriction on the fidelity of the mesh. Recall, that the goal is to achieve a reduction in the cost functional and not to find the global minimum of the non-convex optimization problem. Lastly, in all examples we provide both the surface area error, ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘› ), the Hausdorff measure ๐ป (ฮ“๐บ , ฮ“๐‘› ) and the value of the cost functional at the given iteration count. Furthermore, to illustrate the difference between the reconstructed and the goal surfaces we added to the figures their projections onto the ๐‘ฅ-๐‘ง, ๐‘ฆ-๐‘ง, and ๐‘ฅ-๐‘ฆ planes, which are indicated as black meshed and gray solid regions, respectively. 6.4.1 Shape Reconstruction of a Bean In the first example, we consider a bean that fits in a 2.71 m ร— 1.27 m ร— 1.18 m bounding box as our target shape, see Fig. 6.4e. The main point of this example is to highlight the comparison between an arbitrary starting point against a VSRM initial guess for our shape reconstruction algorithm. In particular, we conduct two different experiments: the first is to reconstruct the target shape starting from an arbitrary starting point, in this case a sphere of radius 1.0 m and the second, is to use VSRM to generate an initial guess, see Fig. 6.4a. When starting with a sphere, our goal fields are generated using the procedure described earlier using 9 incident fields at a set of directions D1 = (๐œƒ1 , ๐œ™1 ) . The frequency of these  fields are ๐‘“1 = 540๐ป๐‘ง. The initial starting mesh for the sphere is constructed using 1802 vertices and 3600 faces. The reconstruction is performed using 10 MHs. The second experiment is to start from a VSRM guess; we use 16 incident plane 102 waves defined by D2 = (๐œƒ2 , ๐œ™2 ) , and frequencies {100๐ป๐‘ง, 125๐ป๐‘ง, 150๐ป๐‘ง, ยท ยท ยท , 700๐ป๐‘ง}  to generate this initial guess. The starting mesh using the initial guess is discretized at 1202 vertices and 2400 faces. Similarly, the goal field used for the reconstruction process is generated from a single field operating at ๐‘“2 = 540๐ป๐‘ง traveling in the โˆ’yฬ‚ direction. Lastly, we use 20 MHs for the reconstruction routine. As seen in Figs. 6.4c and 6.4d, the algorithm was able to reconstruct the object and its features accurately from noisy data; notably the concavity of the object were successfully recovered in both cases. It should be noted that the algorithm was able to recover the correct size and width of the scattering object. Their relative placement and generic dimensions agree very well with those of the target. No. of iter. Initial 10 20 25 33 ๐ป (ฮ“๐บ , ฮ“๐‘› ) 5.37E-1 1.96E-2 8.21E-2 4.23E-2 4.87E-2 ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘› ) 5.75E-1 8.06E-2 3.15E-2 5.15E-3 8.97E-3 ๐ฝ(๐›ฝ ๐‘› ) 6.08 6.04E-1 1.42E-2 9.15E-2 3.96E-3 Table 6.1: Evolution of the error metrics with respect to iteration for the bean (Fig. 6.4e) starting from a sphere (Fig. 6.4b). No. of iter. Initial 10 20 25 33 ๐ป (ฮ“๐บ , ฮ“๐‘› ) 1.91E-1 1.73E-1 8.26E-1 8.49E-2 8.03E-2 ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘› ) 1.24E-1 8.28E-2 3.45E-2 3.07E-2 2.28E-2 ๐ฝ(๐›ฝ ๐‘› ) 2.57E-1 1.09E-1 1.06E-2 8.14E-2 6.60E-3 Table 6.2: Evolution of the error metrics with respect to iteration for the bean (Fig. 6.4e) starting from the initial guess (Fig. 6.4a). 103 (a) Initial Guess using VSRM. (b) Initial Guess Sphere. (c) Final shape (VSRM). (d) Final shape (Sphere). (e) Target shape. Figure 6.4: Shape reconstruction of a bean (6.4e): starting from the initial guess (6.4a) it converges to the reconstructed shape (6.4c) after 33 iterations; starting from a sphere (6.4b) it converges to the reconstructed shape (6.4d) after 33 iterations. 6.4.2 Shape Reconstruction of a Bumpy Cube Next, we examine the performance of this algorithm on a more complex/challenging target illustrated in Fig. 6.5c. This geometry contains both convex and non-convex features which demonstrates the efficiency of MHs representation as well present an overall challenging reconstruction problem. In this case, we only use VSRM to obtain an initial 104 guess. This is done using 16 incident planes waves in directions D2 = (๐œƒ2 , ๐œ™2 ) and  frequencies {100๐ป๐‘ง, 125๐ป๐‘ง, ยท ยท ยท , 800๐ป๐‘ง}. The goal fields are generated for a set of 9 incident fields D1 = (๐œƒ1 , ๐œ™1 ) operating at ๐‘“1 = 400 Hz. As such, the initial mesh is  discretized with 1202 vertices and 2400 patches. The reconstruction routine was done using 40 MHs. The final result is shown in Fig. 6.5b; the algorithm was able to reconstruct the object and its features accurately from noisy data. Notably, the concave features of the object, i.e., the protruding lobes, were successfully recovered. Furthermore, the algorithm was able to recover the correct size and width of the scattering object. Again, we find their relative placement and generic dimensions agree very well with the target. No. of iter. Initial 5 10 15 19 ๐ป (ฮ“๐บ , ฮ“๐‘› ) 2.48E-1 1.68E-1 9.96E-2 5.72E-2 4.18E-2 ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘› ) 1.15E-1 2.76E-3 1.90E-2 2.17E-2 1.41E-2 ๐ฝ(๐›ฝ ๐‘› ) 3.21 1.63 1.45E-1 7.36E-2 4.44E-2 Table 6.3: Evolution of the error metrics with respect to iteration for the bumpy cube (Fig. 6.5c) starting from the initial guess (Fig. 6.5a). 6.4.3 Multiresolution Shape Reconstruction In this final example, we consider the reconstruction of a cow that fits into a 1.7 m x 1.7 m x 1.0 m box, see Fig. 6.6c. This geometry is multi-scale, implying there are fine-scale surface features as well as coarse features, which can be utilized for demonstrating the efficacy of the multi-resolution feature of our reconstruction scheme. Furthermore, given the multi-scale surface features, both the initial guess and the reconstruction routine require a richer set of data relative to the previous examples. In this case, the initial guess, see Fig. 6.6a, is constructed from 16 incident plane waves in directions D2 = (๐œƒ2 , ๐œ™2 ) , for the  frequencies {100๐ป๐‘ง, 150๐ป๐‘ง, ยท ยท ยท , 1200๐ป๐‘ง๐ป๐‘ง}. The shape is optimized at two independent frequencies, 550 Hz and 900 Hz. Again, to maintain the prescribed number of samples over the mesh as ๐œ†/7, we discretize the mesh at 1102 vertices and 2200 faces for the first frequency and and then 2402 vertices and 4800 faces for the second frequency. Furthermore, 105 (a) Initial Guess. (b) Final Shape. (c) Target Shape. Figure 6.5: Reconstruction of the bumpy cube (Fig. 6.5c) starting from an initial guess (Fig. 6.5a) converging to the reconstructed surface (Fig. 6.5b) after 19 iterations. 106 at 550 Hz, we use 100 MHs and at 900 Hz we use 150 MHs. At each of these frequencies, we optimize in bands of 50 MHs. This implies that we choose ๐‘‘๐‘–๐‘š {โ„ฌ๐‘– } = 50 for any ๐‘–. Reconstruction is done sequentially, in that we choose a threshold error for โ„ฌ1 and once that is reached, we optimize โ„ฌ2 and so on. Once we reach a local minima for the last band, we take our current optimized shape and refined it, such that we can optimize for the next frequency at the proper sampling rate and so on. As is evident from Fig. 6.6b, the algorithm recovers the major significant features of the object. While it does recover all features, locations and shape, it does not completely capture the legs or ear to high fidelity. This is largely due to reconstruction at two frequencies only that are not sufficient to resolve rather small features. No. of iter. Initial 29 39 52 63 ๐ป (ฮ“๐บ , ฮ“๐‘› ) 1.20E-1 1.09E-1 1.01E-1 1.29E-1 9.00E-2 ๐‘† ๐‘’๐‘Ÿ๐‘Ÿ (ฮ“๐บ , ฮ“๐‘› ) 3.41E-2 1.05E-2 5.84E-3 2.95E-2 7.87E-3 ๐ฝ(๐›ฝ ๐‘› ) 6.69E-1 3.15E-2 1.93E-2 5.10E-1 2.44E-1 Table 6.4: Evolution of the error metrics with respect to iteration for the multiresolution target surface (Fig. 6.6c) starting from the initial guess (Fig. 6.6a). 6.5 Conclusion In this chapter, we proposed a novel optimization scheme for shape reconstruction; the crux of our contributions rely on using manifold harmonics as a foundation for shape reconstruction. We have exploited the compression and multi-resolution framework that these harmonics provide to address every facet of the reconstruction process; from representing geometry to physics on the manifold. These harmonics rely on an underlying isogeometric analysis framework built on subdivision basis sets. A number of results demonstrate the viability and efficiency of the proposed approach on a set of challenging targets. MHB proves to be a robust tool with a number of applications for inverse design; these include developing maps directly between geometry and far-fields, multi-resolution editing, integration with machine learning, and so on. 107 (a) Initial Guess. (b) Final Shape. (c) Target Shape. Figure 6.6: Reconstruction of target surface (Fig. 6.6c) starting from an initial guess (6.6a) converging to the reconstructed surface (Fig. 6.6c) after 63 iterations. 108 This chapter, ยฉ2022 JASA, is reprinted, with permission, from A. M. A. Alsnayyan and B. Shanker, "Laplace-Beltrami based multi-resolution shape reconstruction on subdivision surfaces," The Journal of the Acoustical Society of America on, March 2022, with slight modifications to fit in this dissertation. 109 BIBLIOGRAPHY [1] Robert J Adams. Physical and analytical properties of a stabilized electric field integral equation. IEEE Transactions on Antennas and Propagation, 52(2):362โ€“372, 2004. [2] Simon Adrian. Rapidly converging boundary integral equation solvers in computational electromagnetics. PhD thesis, Ecole nationale supรฉrieure Mines-Tรฉlรฉcom Atlantique Bretagne Pays de la Loire, 2018. [3] Yonathan Aflalo, Haim Brezis, and Ron Kimmel. On the optimality of shape and data representation in the spectral domain. SIAM Journal on Imaging Sciences, 8(2):1141โ€“1160, 2015. [4] Dhwanit Agarwal, Michael Oลƒeil, and Manas Rachh. Fast multipole acceler- ated solvers for the laplace-beltrami problem in three dimensions. arXiv preprint arXiv:2111.10743, 2021. [5] AMA Alsnayyan, L Kempel, and B Shanker. A complete helmholtz decomposition on multiply connected subdivision surfaces and its application to integral equations. arXiv preprint arXiv:2201.13223, 2022. [6] AMA Alsnayyan, J Li, S Hughey, A Diaz, and B Shanker. Efficient isogeometric boundary element method for analysis of acoustic scattering from rigid bodies. The Journal of the Acoustical Society of America, 147(5):3275โ€“3284, 2020. [7] AMA Alsnayyan and B Shanker. Iso-geometric integral equation solvers and their compression via manifold harmonics. arXiv preprint arXiv:2106.11907, 2021. [8] AMA Alsnayyan and B Shanker. Manifold harmonics for computational electromag- netics. In IEEE International Conference on Microwaves, Antennas, Biomedical Engineering & Electronic Systems (COMCAS), 2021. [9] AMA Alsnayyan and B Shanker. Laplace-beltrami based multi-resolution shape reconstruction on subdivision surfaces. The Journal of the Acoustical Society of America, 151(3):2207โ€“2222, 2022. [10] Yuri Montesinos Alvarez, Fernando Las-Heras, and Cebriรกn Garcรญa. The sources reconstruction method for antenna diagnostics and imaging applications. 2012. [11] S Amini and PJ Harris. A comparison between various boundary integral formulations of the exterior acoustic problem. Computer methods in applied mechanics and engineering, 84(1):59โ€“75, 1990. [12] Siamak Amini. On the choice of the coupling parameter in boundary integral formulations of the exterior acoustic problem. Applicable analysis, 35(1-4):75โ€“92, 1990. [13] Francesco P Andriulli, Kristof Cools, Hakan Bagci, Femke Olyslager, Annalisa Buffa, Snorre Christiansen, and Eric Michielssen. A multiplicative calderon preconditioner for the electric field integral equation. IEEE Transactions on Antennas and Propagation, 56(8):2398โ€“2412, 2008. 110 [14] Xavier Antoine and Marion Darbas. An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems. Multiscale Science and Engineering, 3(1):1โ€“35, 2021. [15] Gregory Arden. Approximation properties of subdivision surfaces. University of Wash- ington, 2001. [16] Jasbir S Arora. An exposition of the material derivative approach for structural shape sensitivity analysis. Computer methods in applied mechanics and engineering, 105(1):41โ€“62, 1993. [17] Ivo Babuลกka and Jens M Melenk. The partition of unity method. International journal for numerical methods in engineering, 40(4):727โ€“758, 1997. [18] Kosala Bandara and Fehmi Cirak. Isogeometric shape optimisation of shell structures using multiresolution subdivision surfaces. Computer-Aided Design, 95:62โ€“71, 2018. [19] Kosala Bandara, Fehmi Cirak, Gรผnther Of, Olaf Steinbach, and Jan Zapletal. Bound- ary element based multiresolution shape optimisation in electrostatics. Journal of computational physics, 297:584โ€“598, 2015. [20] Kosala Bandara, Thomas Rรผberg, and Fehmi Cirak. Shape optimisation with multiresolution subdivision surfaces and immersed finite elements. Computer Methods in Applied Mechanics and Engineering, 300:510โ€“539, 2016. [21] Erik Bรคngtsson, Daniel Noreland, and Martin Berggren. Shape optimization of an acoustic horn. Computer methods in applied mechanics and engineering, 192(11-12):1533โ€“ 1571, 2003. [22] Pieter J Barendrecht. Isogeometric analysis for subdivision surfaces. Eindhoven University of Technology: Eindhoven, The Netherlands, 2013. [23] Y. Bazilevs, V. M. Calo, J. A. Cottrell, J. A. Evans, T. J. R. Hughes, S. Lipton, M. A. Scott, and T. W. Sederberg. Isogeometric analysis using t-splines. Computer methods in applied mechanics and engineering, 199(5):229โ€“263, 2010. [24] Yuri Bazilevs, Victor M Calo, John A Cottrell, John A Evans, Thomas Jr R Hughes, S Lipton, Michael A Scott, and Thomas W Sederberg. Isogeometric analysis using t-splines. Computer Methods in Applied Mechanics and Engineering, 199(5-8):229โ€“263, 2010. [25] Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373โ€“1396, 2003. [26] Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani. Manifold regularization: A geometric framework for learning from labeled and unlabeled examples. Journal of machine learning research, 7(11), 2006. [27] K-U Bletzinger, Stefan Kimmich, and Ekkehard Ramm. Efficient modeling in shape optimal design. Computing Systems in Engineering, 2(5-6):483โ€“495, 1991. 111 [28] Anders Bondeson, Yongtao Yang, and Per Weinerfelt. Shape optimization for radar cross sections by a gradient method. International journal for numerical methods in engineering, 61(5):687โ€“715, 2004. [29] Mario Botsch, Leif Kobbelt, Mark Pauly, Pierre Alliez, and Bruno Lรฉvy. Polygon mesh processing. CRC press, 2010. [30] Yassine Boubendir and Catalin Turc. Well-conditioned boundary integral equation formulations for the solution of high-frequency electromagnetic scattering problems. Computers & Mathematics with Applications, 67(10):1772โ€“1805, 2014. [31] Cesare Bracco, Annalisa Buffa, Carlotta Giannelli, and Rafael Vรกzquez. Adaptive isogeometric methods with hierarchical splines: an overview. Discrete And Continuous Dynamical Systems, 39(ARTICLE):241โ€“261, 2019. [32] Oscar Bruno, Tim Elling, Randy Paffenroth, and Catalin Turc. Electromagnetic integral equations requiring small numbers of krylov-subspace iterations. Journal of Computational Physics, 228(17):6169โ€“6183, 2009. [33] Oscar Bruno, Tim Elling, and Catalin Turc. Well-conditioned high-order algorithms for the solution of three-dimensional surface acoustic scattering problems with neumann boundary conditions. J. Numer. Meth. Eng, 91(10):1โ€“23, 2009. [34] Oscar Bruno, Tim Elling, and Catalin Turc. Fast high-order algorithms and well- conditioned integral equations for high-frequency sound-hard scattering problems. 2012. [35] Annalisa Buffa and Snorre Christiansen. A dual finite element complex on the barycentric refinement. Mathematics of Computation, 76(260):1743โ€“1769, 2007. [36] AJ Burton and GF495032 Miller. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 323(1553):201โ€“210, 1971. [37] Wei Cai, Tiejun Yu, Han Wang, and Yฤณin Yu. High-order mixed rwg basis functions for electromagnetic applications. IEEE Transactions on Microwave Theory and Techniques, 49(7):1295โ€“1303, 2001. [38] Edwin Catmull and James Clark. Recursively generated b-spline surfaces on arbitrary topological meshes. Computer-aided design, 10(6):350โ€“355, 1978. [39] Leilei Chen, Chuang Lu, Wenchang Zhao, Haibo Chen, and Changjun Zheng. Subdi- vision surfacesโ€”boundary element accelerated by fast multipole for the structural acoustic problem. Journal of Theoretical and Computational Acoustics, 28(02):2050011, 2020. [40] Weng Cho Chew. Vector potential electromagnetics with generalized gauge for inhomogeneous media: Formulation. Progress In Electromagnetics Research, 149:69โ€“84, 2014. 112 [41] Weng Cho Chew, Eric Michielssen, JM Song, and Jian-Ming Jin. Fast and efficient algorithms in computational electromagnetics. Artech House, Inc., 2001. [42] L Chommeloux, Christian Pichot, and J-C Bolomey. Electromagnetic modeling for microwave imaging of cylindrical buried inhomogeneities. IEEE Transactions on Microwave Theory and Techniques, 34(10):1064โ€“1076, 1986. [43] Paolo Cignoni, Marco Callieri, Massimiliano Corsini, Matteo Dellepiane, Fabio Ganovelli, Guido Ranzuglia, et al. Meshlab: an open-source mesh processing tool. In Eurographics Italian chapter conference, volume 2008, pages 129โ€“136. Salerno, Italy, 2008. [44] Fehmi Cirak, Michael Ortiz, and Peter Schrรถder. Subdivision surfaces: a new paradigm for thin-shell finite-element analysis. International Journal for Numerical Methods in Engineering, 47(12):2039โ€“2072, 2000. [45] Fehmi Cirak, Michael J Scott, Erik K Antonsson, Michael Ortiz, and Peter Schrรถder. Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision. Computer-Aided Design, 34(2):137โ€“148, 2002. [46] David Colton and Rainer Kress. Integral equation methods in scattering theory. SIAM, 2013. [47] David L Colton, Rainer Kress, and Rainer Kress. Inverse acoustic and electromagnetic scattering theory, volume 93. Springer, 1998. [48] Harry Contopanagos, Benjamin Dembart, Michael Epton, John J Ottusch, Vladimir Rokhlin, John L Visher, and Stephen M Wandzura. Well-conditioned boundary integral equations for three-dimensional electromagnetic scattering. IEEE Transactions on Antennas and Propagation, 50(12):1824โ€“1830, 2002. [49] J Austin Cottrell, Alessandro Reali, Yuri Bazilevs, and Thomas JR Hughes. Isogeo- metric analysis of structural vibrations. Computer methods in applied mechanics and engineering, 195(41-43):5257โ€“5296, 2006. [50] Marion Darbas. Generalized combined field integral equations for the iterative solution of the three-dimensional maxwell equations. Applied Mathematics Letters, 19(8):834โ€“839, 2006. [51] D. Dault and B. Shanker. A mixed potential MLFMA for higher order moment methods with application to the generalized method of moments. IEEE Trans. Antennas Propag., 64(2):650โ€“662, 2016. [52] Daniel L Dault, Naveen V Nair, Jie Li, and Balasubramaniam Shanker. The general- ized method of moments for electromagnetic boundary integral equations. IEEE transactions on antennas and propagation, 62(6):3174โ€“3188, 2014. [53] Fernando De Goes, Mathieu Desbrun, Mark Meyer, and Tony DeRose. Subdivision exterior calculus for geometry processing. ACM Transactions on Graphics (TOG), 35(4):1โ€“11, 2016. 113 [54] Alexandre Dรฉly, Adrien Merlini, Simon B Adrian, and Francesco P Andriulli. On preconditioning electromagnetic integral equations in the high frequency regime via helmholtz operators and quasi-helmholtz projectors. In 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA), pages 1338โ€“1341. IEEE, 2019. [55] Tony DeRose, Michael Kass, and Tien Truong. Subdivision surfaces in character animation. In Proceedings of the 25th annual conference on Computer graphics and interactive techniques, pages 85โ€“94, 1998. [56] Jรผrgen Dรถlz, Stefan Kurz, Sebastian Schรถps, and Felix Wolf. Isogeometric boundary elements in electromagnetism: rigorous analysis, fast methods, and examples. SIAM Journal on Scientific Computing, 41(5):B983โ€“B1010, 2019. [57] Daniel Doo and Malcolm Sabin. Behaviour of recursive division surfaces near extraordinary points. Computer-Aided Design, 10(6):356โ€“360, 1978. [58] Michael R Dรถrfel, Bert Jรผttler, and Bernd Simeon. Adaptive isogeometric analysis by local h-refinement with t-splines. Computer methods in applied mechanics and engineering, 199(5-8):264โ€“275, 2010. [59] Michael R. Dรถrfel, Bert Jรผttler, and Bernd Simeon. Adaptive isogeometric analysis by local h-refinement with t-splines. Computer methods in applied mechanics and engineering, 199(5):264โ€“275, 2010. [60] Michael Driscoll, Katherine A Yelick, and Armando Fox. Subdivision surface evaluation as sparse matrix-vector multiplication. PhD thesis, Masterโ€™s thesis, EECS Department, University of California, Berkeley, 2014. [61] Mehdi Ebrahimi and Alireza Jahangirian. Aerodynamic optimization of airfoils using adaptive parameterization and genetic algorithm. Journal of Optimization Theory and Applications, 162(1):257โ€“271, 2014. [62] Charles L Epstein and Leslie Greengard. Debye sources and the numerical solution of the time harmonic maxwell equations. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 63(4):413โ€“ 463, 2010. [63] EJ Evans, MA Scott, X Li, and DC Thomas. Hierarchical t-splines: Analysis-suitability, bรฉzier extraction, and application as an adaptive basis for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 284:1โ€“20, 2015. [64] Reza Firoozabadi, Eric L Miller, Carey M Rappaport, and Ann W Morgenthaler. Subsurface sensing of buried objects under a randomly rough surface using scattered electromagnetic field data. IEEE Transactions on Geoscience and Remote Sensing, 45(1):104โ€“117, 2006. [65] Matthias Fischer, Ute Gauger, and Lothar Gaul. A multipole galerkin boundary element method for acoustics. Engineering analysis with boundary elements, 28(2):155โ€“ 162, 2004. 114 [66] PC Fourie and Albert A Groenwold. The particle swarm optimization algorithm in size and shape optimization. Structural and Multidisciplinary Optimization, 23(4):259โ€“ 267, 2002. [67] Xin Fu, Jie Li, Li Jiang, and Balasubramaniam Shanker. Generalized debye sources- based efie solver on subdivision surfaces. IEEE Transactions on Antennas and Propaga- tion, 65:5376โ€“5386, 2017. [68] LOTHAR GAUL, MATTHIAS FISCHER, and UTE GAUGER. Accuracy and efficiency of the multipole galerkin bem for acoustics. In Theoretical And Computational Acoustics 2003, pages 97โ€“106. World Scientific, 2004. [69] Tristan Goodwill and Michael Oโ€™Neil. On the numerical solution of the laplace- beltrami problem on piecewise-smooth surfaces. arXiv preprint arXiv:2108.08959, 2021. [70] Roberto D Graglia and Guido Lombardi. Singular higher order complete vector bases for finite methods. IEEE transactions on antennas and propagation, 52(7):1672โ€“1685, 2004. [71] Roberto D Graglia, Donald R Wilton, and Andrew F Peterson. Higher order interpolatory vector bases for computational electromagnetics. IEEE transactions on antennas and propagation, 45(3):329โ€“342, 1997. [72] Massimo Guiggiani and A. Gigante. A general algorithm for multidimensional cauchy principal value integrals in the boundary element method. Journal of Applied Mechanics-transactions of The Asme - J APPL MECH, 57, 1990. [73] Massimo Guiggiani, Guna Krishnasamy, Thomas Rudolphi, and F. Rizzo. A general algorithm for the numerical solution of hypersingular boundary integral equations. Journal of Applied Mechanics, 59, 1992. [74] Nail A Gumerov and Ramani Duraiswami. A broadband fast multipole accelerated boundary element method for the three dimensional helmholtz equation. The Journal of the Acoustical Society of America, 125(1):191โ€“205, 2009. [75] Raphael T Haftka and Ramana V Grandhi. Structural shape optimizationโ€”a survey. Computer methods in applied mechanics and engineering, 57(1):91โ€“106, 1986. [76] Roger F Harrington. Time-harmonic electromagnetic fields. McGraw-Hill College, 1961. [77] Behrooz Hassani, S Mehdi Tavakkoli, and NZ Moghadam. Application of isoge- ometric analysis in structural shape optimization. Scientia Iranica, 18(4):846โ€“852, 2011. [78] George C Hsiao and Ralph E Kleinman. Mathematical foundations for error es- timation in numerical solutions of integral equations in electromagnetics. IEEE transactions on Antennas and Propagation, 45(3):316โ€“328, 1997. 115 [79] Thomas JR Hughes, John A Cottrell, and Yuri Bazilevs. Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement. Computer methods in applied mechanics and engineering, 194(39-41):4135โ€“4195, 2005. [80] S. Hughey, H. M. Aktulga, V. Melapudi, L. Mingyu, B. Shanker, and E. Michielssen. Parallel wideband MLFMA for analysis of electrically large, nonuniform, multiscale structures. IEEE Trans. Antennas Propag., 67(2):1094โ€“1107, 2019. [81] Olha Ivanyshyn and Rainer Kress. Identification of sound-soft 3d obstacles from phaseless data. Inverse Problems & Imaging, 4(1):131, 2010. [82] Steven G. Johnson. The Nlopt nonlinear-optimization package. [83] Erik Jorgensen, John L Volakis, Peter Meincke, and Olav Breinbjerg. Higher order hierarchical legendre basis functions for electromagnetic modeling. IEEE Transactions on Antennas and Propagation, 52(11):2985โ€“2995, 2004. [84] Sokolowski J Zolรฉsio JP. Introduction to shape optimization, shape sensitivity analysis, 1992. [85] Bert Jรผttler, Angelos Mantzaflaris, Ricardo Perl, and Martin Rumpf. On isogeometric subdivision methods for pdes on surfaces. arXiv preprint arXiv:1503.03730, 2015. [86] Zachi Karni and Craig Gotsman. Spectral compression of mesh geometry. In Proceedings of the 27th annual conference on Computer graphics and interactive techniques, pages 279โ€“286, 2000. [87] Nam H Kim and Jun Dong. Shape sensitivity analysis of sequential structuralโ€“ acoustic problems using fem and bem. Journal of sound and vibration, 290(1-2):192โ€“208, 2006. [88] Seung-Goo Kim, Moo K Chung, Seongho Seo, Stacey M Schaefer, Carien M van Reekum, and Richard J Davidson. Heat kernel smoothing via laplace-beltrami eigenfunctions and its application to subcortical structure modeling. In Pacific-Rim Symposium on Image and Video Technology, pages 36โ€“47. Springer, 2011. [89] Stephen Martin Kirkup. The boundary element method in acoustics. Integrated sound software, 2007. [90] ANDREAS Kirsch. The domain derivative and two applications in inverse scattering theory. Inverse problems, 9(1):81, 1993. [91] RE Kleinman and GF Roach. Boundary integral equations for the three-dimensional helmholtz equation. SIAM review, 16(2):214โ€“236, 1974. [92] Alena Kopanicฬ†ร kovร , Rolf Krause, and Rasmus Tamstorf. Subdivision-based nonlin- ear multiscale cloth simulation. SIAM Journal on Scientific Computing, 41(5):S433โ€“S461, 2019. 116 [93] Rainer Kress and Axel Zinn. On the numerical solution of the three-dimensional inverse obstacle scattering problem. Journal of computational and applied mathematics, 42(1):49โ€“61, 1992. [94] K. J. Langenberg, K. Mayer, P. Fellinger, and R. Marklein. Imaging and Inverse Scattering in Nondestructive Evaluation with Acoustic and Elastic Waves, chapter 22, pages 165โ€“172. Springer US, 1993. [95] Frรฉdรฉrique Le Louรซr. A spectrally accurate method for the direct and inverse scattering problems by multiple 3d dielectric obstacles. ANZIAM Journal, 59:1โ€“49, 2017. [96] Bruno Lรฉvy. Laplace-beltrami eigenfunctions towards an algorithm that" under- stands" geometry. In IEEE International Conference on Shape Modeling and Applications 2006 (SMIโ€™06), pages 13โ€“13. IEEE, 2006. [97] J. Li, D. Dault, and B. Shanker. New trends in computational electromagnetics, chapter New trends in geometric modeling and discretization of integral equations, pages 315โ€“372. Institute of Engineering and Technology, 2019. [98] Jie Li, Daniel Dault, Beibei Liu, Yiying Tong, and Balasubramaniam Shanker. Subdi- vision based isogeometric analysis technique for electric field integral equations for simply connected structures. Journal of Computational Physics, 319:145โ€“162, 2016. [99] Jie Li, Daniel Dault, and Balasubramaniam Shanker. A quasianalytical time domain mie solution for scattering from a homogeneous sphere. arXiv preprint arXiv:1311.1237, 2013. [100] Jie Li, Xin Fu, and Balasubramaniam Shanker. Formulation and iso-geometric analysis of scalar integral equations for electromagnetic scattering. IEEE Transactions on Antennas and Propagation, 66(4):1957โ€“1966, 2018. [101] Jie Li and Balasubramaniam Shanker. Isogeometric analysis of em scattering on multiply-connected subdivision surfaces. In 2017 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, pages 1557โ€“1558. IEEE, 2017. [102] Haojie Lian, P Kerfriden, and SPA361274807185767 Bordas. Shape optimization directly from cad: An isogeometric boundary element approach using t-splines. Computer Methods in Applied Mechanics and Engineering, 317:1โ€“41, 2017. [103] Shujin Lin, Fang You, Xiaonan Luo, and Zheng Li. Deducing interpolating subdivision schemes from approximating subdivision schemes. ACM Transactions on Graphics (TOG), 27(5):1โ€“7, 2008. [104] Tzu-Chu Lin. A proof for the burton and miller integral equation approach for the helmholtz equation. Journal of mathematical analysis and applications, 103(2):565โ€“574, 1984. 117 [105] Zhaowei Liu, Musabbir Majeed, Fehmi Cirak, and Robert N. Simpson. Isogeometric fem-bem coupled structural-acoustic analysis of shells using subdivision surfaces. International Journal for Numerical Methods in Engineering, 113(9):1507โ€“1530, 2018. [106] C. Loop. Smooth subdivision surfaces based on triangles. Masterโ€™s thesis, University of Utah, Department of Mathematics, 1987. [107] Charles Loop. Smooth spline surfaces over irregular meshes. In Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 303โ€“310, 1994. [108] Michael Lounsbery, Tony D DeRose, and Joe Warren. Multiresolution analysis for surfaces of arbitrary topological type. ACM Transactions on Graphics (TOG), 16(1):34โ€“73, 1997. [109] Filippo Maggioli, Simone Melzi, Maks Ovsjanikov, Michael M Bronstein, and Emanuele Rodolร . Orthogonalized fourier polynomials for signal approximation and transfer. In Computer Graphics Forum, volume 40, pages 435โ€“447. Wiley Online Library, 2021. [110] Steffen Marburg. The burton and miller method: Unlocking another mystery of its coupling parameter. Journal of Computational Acoustics, 24(01):1550016, 2016. [111] Jose Angel Martinez-Lorenzo, Carey M Rappaport, and Fernando Quivira. Physical limitations on detecting tunnels using underground-focusing spotlight synthetic aperture radar. IEEE Transactions on Geoscience and Remote Sensing, 49(1):65โ€“70, 2010. [112] V. Melapudi, H. Huang, B. Shanker, and T. Van. A novel wideband FMM for fast integral equation solution of multiscale problems in electromagnetics. IEEE Trans. Antennas Propag., 57(7):2094โ€“2104, 2009. [113] Simone Melzi, Emanuele Rodolร , Umberto Castellani, and Michael M Bronstein. Localized manifold harmonics for spectral shape analysis. In Computer Graphics Forum, volume 37, pages 20โ€“34. Wiley Online Library, 2018. [114] Rajendra Mitharwal and Francesco P Andriulli. On the multiplicative regularization of graph laplacians on closed and open structures with applications to spectral partitioning. IEEE Access, 2:788โ€“796, 2014. [115] Farhan Mohamed and C Vei Siang. A survey on 3d ultrasound reconstruction techniques. Artificial Intelligenceโ€”Applications in Medicine and Biology, 2019. [116] NV Nair, B Shanker, and L Kempel. Generalized method of moments: A boundary integral framework for adaptive analysis of acoustic scattering. The Journal of the Acoustical Society of America, 132(3):1261โ€“1270, 2012. [117] Jean-Claude Nรฉdรฉlec. Mixed finite elements in r3. Numerische Mathematik, 35(3):315โ€“ 341, 1980. [118] Jean-Claude Nรฉdรฉlec. Acoustic and electromagnetic equations: integral representations for harmonic problems. Springer Science & Business Media, 2001. 118 [119] Vinh Phu Nguyen, Cosmin Anitescu, Stรฉphane PA Bordas, and Timon Rabczuk. Iso- geometric analysis: an overview and computer implementation aspects. Mathematics and Computers in Simulation, 117:89โ€“116, 2015. [120] Branislav M Notaros. Higher order frequency-domain computational electromagnet- ics. IEEE Transactions on Antennas and Propagation, 56(8):2251โ€“2276, 2008. [121] Qing Pan, Timon Rabczuk, Gang Xu, and Chong Chen. Isogeometric analysis for surface pdes with extended loop subdivision. Journal of Computational Physics, 398:108892, 2019. [122] Giuseppe Patanรจ. Laplacian spectral basis functions. Computer Aided Geometric Design, 65:31โ€“47, 2018. [123] Joris Peeters, Kristof Cools, Ignace Bogaert, Femke Olyslager, and Daniรซl De Zutter. Embedding calderรณn multiplicative preconditioners in multilevel fast multipole algorithms. IEEE transactions on antennas and propagation, 58(4):1236โ€“1250, 2010. [124] Niklas Peinecke, Franz-Erich Wolter, and Martin Reuter. Laplace spectra as finger- prints for image recognition. Computer-Aided Design, 39(6):460โ€“476, 2007. [125] Jรถrg Peters and Xiaobin Wu. On the local linear independence of generalized subdivision functions. SIAM journal on numerical analysis, 44(6):2389โ€“2407, 2006. [126] Andrew F Peterson. Mapped vector basis functions for electromagnetic integral equations. Synthesis Lectures on Computational Electromagnetics, 1(1):1โ€“124, 2005. [127] Andrew F. Peterson, Scott L. Ray, and Raj Mittra. Computational Methods for Electro- magnetics. Wiley-IEEE Press, 1997. [128] Roland Potthast. Frรฉchet differentiability of boundary integral operators in inverse acoustic scattering. Inverse Problems, 10(2):431, 1994. [129] Yahya Rahmat-Samii and Eric Michielssen. Electromagnetic optimization by genetic algorithms. Microwave Journal, 42(11):232โ€“232, 1999. [130] Sadasiva Rao, Donald Wilton, and Allen Glisson. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Transactions on antennas and propagation, 30(3):409โ€“ 418, 1982. [131] Ulrich Reif. A unified approach to subdivision algorithms near extraordinary vertices. Computer Aided Geometric Design, 12(2):153โ€“174, 1995. [132] Ulrich Reif and Peter Schrรถder. Curvature integrability of subdivision surfaces. Advances in Computational Mathematics, 14(2):157โ€“174, 2001. [133] Martin Reuter, Silvia Biasotti, Daniela Giorgi, Giuseppe Patanรจ, and Michela Spag- nuolo. Discrete laplaceโ€“beltrami operators for shape analysis and segmentation. Computers & Graphics, 33(3):381โ€“390, 2009. 119 [134] Martin Reuter, Franz-Erich Wolter, and Niklas Peinecke. Laplaceโ€“beltrami spectra as โ€˜shape-dnaโ€™of surfaces and solids. Computer-Aided Design, 38(4):342โ€“366, 2006. [135] RN Rieben, GH Rodrigue, and DA White. A high order mixed vector finite element method for solving the time dependent maxwell equations on unstructured grids. Journal of Computational Physics, 204(2):490โ€“519, 2005. [136] Emanuele Rodolร , Luca Cosmo, Michael M Bronstein, Andrea Torsello, and Daniel Cremers. Partial functional correspondence. In Computer graphics forum, volume 36, pages 222โ€“236. Wiley Online Library, 2017. [137] Steven Rosenberg and Rosenberg Steven. The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds. Number 31. Cambridge University Press, 1997. [138] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003. [139] Stefan A Sauter and Christoph Schwab. Boundary element methods. In Boundary Element Methods, pages 183โ€“287. Springer, 2010. [140] C Schwab. Numerical mathematics and scientific computation, 1998. [141] Johannes Semmler, Lukas Pflug, Michael Stingl, and Gรผnter Leugering. Shape optimization in electromagnetic applications. In New Trends in Shape Optimization, pages 251โ€“269. Springer, 2015. [142] Kubilay Sertel and John L Volakis. Incomplete lu preconditioner for fmm implemen- tation. Microwave and Optical Technology Letters, 26(4):265โ€“267, 2000. [143] B. Shanker and H. Huang. Accelerated cartesian expansions - A fast method for computing of potentials of the form ๐‘…โˆ’๐œˆ for all real ๐œˆ. JCP, 226(1):732โ€“753, 2007. [144] Meged Shoham, Amir Vaxman, and Mirela Ben-Chen. Hierarchical functional maps between subdivision surfaces. In Computer Graphics Forum, volume 38, pages 55โ€“73. Wiley Online Library, 2019. [145] Frutuoso GM Silva and Abel Joรฃo Padrรฃo Gomes. Interactive editing of multiresolu- tion meshes. In Proceedings. 17th Brazilian Symposium on Computer Graphics and Image Processing, pages 202โ€“209. IEEE, 2004. [146] Robert N Simpson, Michael A Scott, Matthias Taus, Derek C Thomas, and Haojie Lian. Acoustic isogeometric boundary element analysis. Computer Methods in Applied Mechanics and Engineering, 269:265โ€“290, 2014. [147] J. M. Song, C. C. Lu, and W. C. Chew. Mlfma for electromagnetic scattering by large complex objects. IEEE Transactions on Antennas and Propagation, 45:1488โ€“1493, 1997. [148] Andrew P Spann, Hong Zhao, and Eric SG Shaqfeh. Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow. Physics of Fluids, 26(3):031902, 2014. 120 [149] Jos Stam. Evaluation of loop subdivision surfaces. SIGGRAPH 99 Course Notes, 06 2001. [150] Olaf Steinbach. Numerical approximation methods for elliptic boundary value problems: finite and boundary elements. Springer Science & Business Media, 2007. [151] Olaf Steinbach. Boundary integral equations for helmholtz boundary value and transmission problems. Direct and inverse problems in wave propagation and applications, 14:253โ€“292, 2013. [152] Olaf Steinbach and Wolfgang L Wendland. The construction of some efficient precon- ditioners in the boundary element method. Advances in Computational Mathematics, 9(1):191โ€“216, 1998. [153] Gilbert W Stewart. A krylovโ€“schur algorithm for large eigenproblems. SIAM Journal on Matrix Analysis and Applications, 23(3):601โ€“614, 2002. [154] Walid Tabbara, Bernard Duchรชne, Ch Pichot, Dominique Lesselier, Luc Chommeloux, and Nadine Joachimowicz. Diffraction tomography: contribution to the analysis of some applications in microwaves and ultrasonics. Inverse Problems, 4(2):305, 1988. [155] Toru Takahashi, Tatsuro Yamamoto, Yuta Shimba, Hiroshi Isakari, and Toshiro Matsumoto. A framework of shape optimisation based on the isogeometric boundary element method toward designing thin-silicon photovoltaic devices. Engineering with Computers, 35(2):423โ€“449, 2019. [156] T Terai. On calculation of sound fields around three dimensional objects by integral equation methods. Journal of Sound and Vibration, 69(1):71โ€“100, 1980. [157] Bruno Vallet and Bruno Lรฉvy. Spectral geometry processing with manifold harmonics. In Computer Graphics Forum, volume 27, pages 251โ€“260. Wiley Online Library, 2008. [158] Alexandre Vion, Ruth V Sabariego, and Christophe Geuzaine. A model reduction algorithm for solving multiple scattering problems using iterative methods. IEEE transactions on magnetics, 47(5):1470โ€“1473, 2011. [159] Curtis R Vogel. Computational methods for inverse problems. SIAM, 2002. [160] Wolfgang A Wall, Moritz A Frenzel, and Christian Cyron. Isogeometric structural shape optimization. Computer methods in applied mechanics and engineering, 197(33- 40):2976โ€“2988, 2008. [161] Stephen Wandzuraz. The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription. IEEE Antennas and Propagation Magazine, 35(3):7โ€“12, 1993. [162] Daniel S Weile and Xiaobo Wang. Strong singularity reduction for curved patches for the integral equations of electromagnetics. IEEE Antennas and Wireless Propagation Letters, 8:1370โ€“1373, 2009. 121 [163] Yair Weiss, Antonio Torralba, Robert Fergus, D Koller, D Schuurmans, Y Bengio, and L Bottou. Spectral hashing. in advances in neural information processing systems 21. In Proc. of the Twenty-Second Annual Conf. on Neural Information Processing Systems, Vancouver, British Columbia, Canada, pages 8โ€“11, 2008. [164] Bai Xiao, Edwin R Hancock, and Richard C Wilson. Geometric characterization and clustering of graphs using heat kernel embeddings. Image and Vision Computing, 28(6):1003โ€“1021, 2010. [165] Chenkai Xu, Hongwei Lin, Hui Hu, and Yaqi He. Fast calculation of laplace-beltrami eigenproblems via subdivision linear subspace. Computers & Graphics, 97:236โ€“247, 2021. [166] Athanasios D Zacharopoulos, Simon R Arridge, Oliver Dorn, Ville Kolehmainen, and Jan Sikora. Three-dimensional reconstruction of shape and piecewise constant region values for optical tomography using spherical harmonic parametrization and a boundary element method. Inverse Problems, 22(5):1509, 2006. [167] Jan Zapletal and Jiล™รญ Bouchala. Shape optimization and subdivision surface based approach to solving 3d bernoulli problems. Computers & Mathematics with Applications, 78(9):2911โ€“2932, 2019. [168] Hao Zhang, Oliver Van Kaick, and Ramsay Dyer. Spectral mesh processing. In Computer graphics forum, volume 29, pages 1865โ€“1894. Wiley Online Library, 2010. [169] Chang-Jun Zheng, Hai-Bo Chen, Hai-Feng Gao, and Lei Du. Is the burtonโ€“miller formulation really free of fictitious eigenfrequencies? Engineering Analysis with Boundary Elements, 59:43โ€“51, 2015. [170] Denis Zorin. Subdivision for modeling and animation. SIGGRAPH2000 course notes, 2000. [171] Denis Zorin, Peter Schrรถder, and Wim Sweldens. Interpolating subdivision for meshes with arbitrary topology. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pages 189โ€“192, 1996. [172] Denis Zorin, Peter Schrรถder, and Wim Sweldens. Interactive multiresolution mesh editing. In Proceedings of the 24th annual conference on Computer graphics and interactive techniques, pages 259โ€“268, 1997. 122