A Refinement-by-Superposition -Method for (curl)- and (div)-Conforming Discretizations

We present refinement-by-superposition (RBS) hp-refinement infrastructure for computational electromagnetics (CEMs), which permits exponential rates of convergence. In contrast to dominant approaches to hp-refinement for continuous Galerkin methods, which rely on explicit constraint equations, the multilevel strategy presented drastically reduces the implementation complexity. Through the RBS methodology, enforcement of continuity occurs by construction, enabling arbitrary levels of refinement with ease, and without the practical (but not theoretical) limitations of constrained-node refinement. We outline the construction of the RBS hp-method for refinement with <inline-formula> <tex-math notation="LaTeX">${H}$ </tex-math></inline-formula>(curl)- and <inline-formula> <tex-math notation="LaTeX">${H}$ </tex-math></inline-formula>(div)-conforming finite cells. Numerical simulations for the 2-D finite element method (FEM) solution of the Maxwell eigenvalue problem demonstrate the effectiveness of RBS hp-refinement. As an additional goal of this work, we aim to promote the use of mixed-order (low- and high-order) elements in practical CEM applications.

Geometric discretization by quadrilaterals and hexahedra, while significantly more accurate with respect to the number of degrees of freedom (DoFs) than modeling with triangles or tetrahedra [6], presents significant challenge to fully dynamic mesh adaptivity.In refinement with triangles and tetrahedra, local adaptivity directives propagate to a small set of neighboring elements, enabling the insertion of new DoFs without modification to the entire element structure.Similar refinement approaches with quadrilateral and hexahedral cells, however, would dictate global refinement, thereby destroying the utility of h-adaptivity.Inserting transition elements also poses significant challenge to the approximation quality, particularly for vectorial shape functions, which are highly sensitive to the loss of linear independence in the physical domain.Finally, application of a discontinuous Galerkin method, while evading this problem, among its other difficulties, generally requires more DoFs for the same level of accuracy.
As shown in [7], [8], [9], and [10], when the solution satisfies certain regularity conditions, local enrichment of function spaces, also known as p-refinement, enables exponential convergence.When such conditions are not satisfied, however, the benefit of p-refinement is heavily degraded, reduced instead to algebraic convergence as in the case of pure spatial subdivision, or h-refinement.Pure p-refinement, while effective in certain situations (e.g., [11]), is therefore insufficient in general, and as such, a combined approach with both hand p-adaptivity is necessary to achieve exponential convergence for solutions with singularities or non-smooth behavior, motivating the need for more advanced and versatile approaches to h-refinement across application domains.
Previous works in computational electromagnetics (CEMs), for example, have demonstrated the potential of hp-adaptivity through hybrid meshes, for example, in [12].Most typically, however, FEM codes overcome the limitations of hpadaptivity by performing the insertion of constrained nodes, which-in contrast to true DoFs-are constrained to enforce continuity conditions with neighboring elements.Such approaches in CEM have shown significant performance increase and exponential convergence in the presence of singular solutions [13], [14], [15], [16], [17], [18], but at the cost of high implementation complexity, impeding wide-scale adoption.Furthermore, such methods are usually limited in implementation to 1-irregular mesh (i.e., only one hanging node per edge), which, not a severe limitation in practice, prevents arbitrary local refinement steps.Open-source libraries-such as deal.II [19]-have significantly simplified 0018-926X © 2023 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the implementation of hp-refinement codes; yet in some cases it might be inconvenient or undesirable to utilize third-party finite element method (FEM) libraries.
In response, we opt to extend the refinement-bysuperposition (RBS) approach introduced in [20], [21], [22], and [23] for hierarchical basis functions, which demonstrated exponential convergence for scalar problems with C 0 finite elements, to H (curl)-and H (div)-conforming finite elements.Additional studies with C 0 finite elements and the RBS hp-method with adaptivity in [24] further motivate extensions of the method to CEM.Naturally, any informative numerical approximation must satisfy characteristics that, in the limit of their resolution, dictate the convergence of the discrete problem and its solution to the continuous one.In this article, we study, in particular, the Maxwell eigenvalue problem, both for its computational challenges [12], and its utility in forming and consistently evaluating numerical benchmarks.For this model problem, significant numerical analysis undergirds the choice of discretization, namely the proven suitability of H (curl)-conforming finite elements via the discrete compactness property for h- [25], p- [26], and hp-refinement modalities [27].Additional studies, such as in [28], [29], consider H (div)-conforming discretizations and their convergence properties in a similar manner.
While the proposed approach significantly simplifies the implementation of hp-refinement infrastructure for applications in CEM, in contrast to more traditional refineby-replacement (RBR) strategies, the proposed approach decreases the sparsity of the system matrices under both hand p-refinements, whereas RBR only reduces sparsity under p-refinements.The proposed approach is, therefore, less suitable for application to very large problems with vast differences in scales.For applications of the approach to boundary element method (BEM) problems, for example, surface integral equation (SIE) problems in CEM discretized by the method of moments (MoMs), such considerations do not apply given the global nature of the Green's function but would instead concern the increase in integration time due to the overlap of refinement layers.
The remainder of this article is organized as follows.Section II details the construction of the RBS hp-method, covering the enforcement of the required continuity conditions (tangential or normal continuity) and ensuring linear independence after the insertion of descendant refinement layers.Section III examines application to an H (curl)-conforming discretization of the Maxwell eigenvalue problem.We examine a challenging eigenpair with a singular eigenfunction as studied in [30].The presented approach yields exponential convergence of the eigenvalue with respect to the number of DoFs (NDoFs), which, along with the ease of implementation, illustrates the practical value for applications in CEM.

II. RBS: DESCRIPTION AND CONSTRUCTION
Let ⊂ R d be a polyhedral domain, d ∈ {2, 3} (the case of d = 1 being trivial), with its boundary ∂ enclosing a volume of homogeneous or inhomogeneous material.We assume that the problems considered are governed by the source-free (timeharmonic) Maxwell's equations along with accompanying boundary conditions, where E, H, ω, µ, and ε denote, respectively, the electric field, the magnetic field, frequency, permeability, and permittivity.The objective is to generate a sequence of overlapping meshes where on the lth level denotes a nonoverlapping, potentially incomplete covering of via elements K l i associated with an index set I l such that the resulting discretization imparts exponential convergence on the assemblage with B denoting the appropriate solution space, namely H (curl; ) or H (div; ), where These two function spaces, both vital to computational electromagnetic problems in different contexts, are tied via the de Rham complex or sequence [12], shown below for d = 3, which identifies the coincidence of ranges and kernels of the constitutive operators This relationship between H (curl; ) and H (div; ) also ties the approximation structure closely, as discussed throughout the remainder of this section.In fact, in 2-D, the shape functions for H (div)-conforming discretizations may be computed simply as rotations of the H (curl)-conforming shape functions of equivalent order.The objective of this multilevel discretization is in some sense distinct from that of multigrid methods [31], as we consider only the generation of H (curl)-and H (div)-conforming approximation spaces of hp character, and not coupling of coarse-to-fine solution ensembles.However, the following method and the adaptivity structure it affords can augment adaptive multigrid approaches.
On a given level l, we enforce a series of constraints on the cells K l i and K l j , i, j 2) Between elements K l i and K l j , hanging nodes are prohibited.
Regarding the hierarchical meshes, we require only that T 0 forms a complete covering of .Subsequent levels may form a complete covering, which we define as a global refinement.
We also state the following property, which, for the numerical solution of the aforementioned problems, impacts integration and sparsity: For Finally, for levels l and k, l ̸ = k, hanging nodes are permitted.
With an underlying hierarchical H (curl)-or H (div)conforming basis, such as introduced in [32], exponential convergence of the solution given in (4) may be achieved with suitable refinements by a collection of overlay meshes.This RBS approach yields the desired discretization by imposing homogeneous Dirichlet boundary conditions on the boundaries of the inserted descendant cells (i.e., the collection of overlay meshes) [20], [21], [22], [23].Continuity requirements, therefore, may be easily enforced for arbitrary levels of refinements (i.e., n-irregular meshes) and heterogeneity in the chosen orders of the hierarchical basis throughout the mesh.This enables an algorithmically straightforward and low-cost method to add hp-refinement capabilities for H (curl)-and H (div)-conforming discretizations by enforcing, respectively, tangential and normal continuity.
We approach the description explicitly from a 2-D perspective; however, the process generalizes trivially to 3-D.First, we classify each shape function in the following manner.In the case of an H (curl)-conforming discretization, we assign each shape function according to the properties of the nonzero tangential components at the boundary of the cell.These shape functions are classified into three categories: the nodefunctions, that is, those functions with non-zero tangential components at only one node; the edge-functions, that is, those functions with nonzero tangential components along one and only one edge; and the cell-functions, that is, those functions that have no nonzero tangential components on the boundary.The existence of node-functions is only necessary when in 2.5-D, the axial component (i.e., perpendicular to the 2-D plane of the geometry) of the solution is nonzero, as in the examples in [11].The same classification strategy, albeit according to the nonzero normal components at the boundary, is applied for seeking solutions in H (div). Naturally, the difficulty in inserting unknowns rests in the treatment of the edge-functions (and potentially the node-functions), while the cell-functions, which introduce no DoFs influencing the boundary, may be inserted or excised without global considerations.
While we focus on isotropic refinement infrastructure, the superposition-based approach supports anisotropy in p-refinement directives trivially and, with modification, to anisotropic h-refinements [33].For example, for a single h-refinement step applied to one cell in 2-D, four new child cells are inserted one refinement layer above the parent cell according to the constraint of isotropic refinement only.Note that geometrically, the child cells lie in the same physical space as the parent; rather, the designation of "above" is purely conceptual.
We have two simultaneous considerations in the refinement process: continuity and linear independence.Continuity must be enforced due to potential nonuniformity in the polynomial degree of the basis on neighboring cells and the existence (or lack thereof) of child cells on the various refinement levels.Linear independence, on the other hand, is guaranteed by the proper delegation of DoFs between the refinement layers descended from the origin cell.DoFs must be deactivated on the parent cells and activated on the child cells, depending on the refinement levels of the cell and its neighbors.
The refinement and coarsening directives are applied for both cases through the assignment of the DoFs to the geometrical structures as mentioned above (the cells, the edges, the nodes, and, in 3-D, the faces).

A. Enforcing Continuity Requirements
The activation and deactivation of the DoFs on the boundaries of the cells follow a unified procedure based on [20], [21], [22], and [23].For each layer in the discretization, from the origin layer T 0 to the highest refinement layer T n , DoFs are assigned to each of the geometrical elements according to the continuity requirements desired, in this case, tangential continuity for seeking a solution in H (curl) and normal continuity for H (div).Each cell, edge, node (when necessary in 2.5-D), and face (exclusively in 3-D) collects a list of DoFs, both active and inactive.
As the DoFs are accumulated, each one is matched as necessary with the associated shape functions on neighboring cells according to the vectorial direction and multiindex of the associated shape function.As opposed to the boundaryfunctions (i.e., node, face, and edge), the cell-functions automatically satisfy continuity requirements and therefore no special considerations are necessary except for those related to ensuring linear independence across overlapping levels.The boundary-function DoFs are marked as active according to the existence of neighbors in the refinement level and the expansion order of those cells.Both restrictions are handled seamlessly and without distinction as the overall process amounts to traversing the geometrical entities in the discretization, which are assigned their maximal sets of associated DoFs, and activating only the DoFs according to the above compatibility conditions.
We summarize the procedure for activating DoFs based on the continuity requirements as follows.
1) For each cell, edge, node, and face in each refinement layer, collect the associated DoFs.a) For H (curl), associate the DoFs with the relevant geometric entity based on the nonzero tangential components.b) For H (div), associate the DoFs with the relevant geometric entity based on the nonzero normal components.2) Iterate through each refinement layer and each edge, node, and face.a) If a suitable refinement neighbor exists, match the shape functions associated with the adjacent Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.cells, activating only the fully matched DoFs and deactivating the rest.In summary, with efficient data structures for relating the DoFs on the boundaries of cells for the cases of H (curl) and H (div), enforcing continuity requirements for seeking solutions in either function space amounts to straightforward queries between neighboring cells.

B. Eliminating Linear Dependence of the Hierarchical Refinements
To ensure linear independence between overlapping shape functions, we prioritize the highest feasible refinement level possible.For example, the cell-functions, which by definition satisfy the continuity requirements automatically, require deactivation on the parent cell and activation of the DoFs on the child cells.Unlike the handling of the edge-and node-functions, this transfer occurs without any queries to the discretization other than checking if the descendant cells exist.Now, for the edge-and node-functions, additional care is necessary.In this case, the preference to delegate DoFs to the child cells is constrained by the refinement state of one or more neighbors of the cell.In other words, in 2-D, if a parent cell shares an edge with another refined parent cell, the DoFs on the parent edge may be transferred to the corresponding edges on the child refinement layer.Likewise, the deactivation of a node-function on the parent refinement layers requires that the corresponding node is surrounded by refined cells.In other words, as in [20], [21], [22], and [23], active geometrical components may not "overlap" with respect to the refinement layers.
We summarize the activation and deactivation of DoFs as follows.
1) On an h-refinement step, deactivate cell-functions on the parent cell and activate the cell-function DoFs on the child cell.2) If a geometrical component (node, edge, or face) on the descendant layer is active (i.e., it has associated active DoFs), deactivate the corresponding component on the parent layer.According to this procedure, a parent cell sufficiently surrounded by refined cells may be entirely deactivated to ensure linear independence and maximize the resolution of the approximation.In such cases, the sparsity of the resulting system is enhanced.

C. Summary of the Overall Approach
Fig. 1 shows examples of how basis functions are distributed across refinement layers as generated from the linearity and continuity enforcement procedures summarized in Sections II-A and II-B.Similar to the descriptions of the RBS process in 1-D in [20], [21], [22], and [23], Fig. 1(a) summarizes the procedure for a 1-D domain, including the transfer of DoFs associated with lower refinement levels to the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.descendant layers and the ability to choose the expansion order p arbitrarily.Note that in 1-D, we have only the linear class of boundary-functions, that is, at the boundary between two cells (across all the refinement levels), only one active boundary-DoF exists.The hierarchical basis functions illustrated in Fig. 1(a) (and those used in Section III) are based on the maximally orthogonalized basis functions [32].
In 2-D (and 3-D), however, many active DoFs exist on the cell boundaries as a result of employing higher order boundary-functions.Depicted in Fig. 1(b), we demonstrate a similar refinement model as in the 1-D case.Unlike in the 1-D case, the depicted refinement in 2-D results in the enforcement of the domain boundary conditions propagating to the higher refinement levels when available.Furthermore, in this case, many of the parent cells retain a large number of DoFs assigned to the boundary due to the higher order boundary-functions.In Fig. 1(b), such occurrences are denoted by the cells with solid boundaries and transparent interiors, in addition to the matching designations related to the active geometrical components as seen in Fig. 1(a).
Finally, for each cell located on a new refinement layer, an additional mapping between integration space and reference space is introduced, resulting in a succession of mapping operations.Note, however, that regardless of the curvature of the origin cell, all subsequent mappings from the reference cell to the child cell have constant Jacobian determinants and may be handled with ease during integration.What remains is whether the resulting solution space belongs to a subset of H (curl) or H (div).To elucidate the preservation of this property under the RBS paradigm, consider the collection of subdomains { i } associated with each of the leaf cells in T .For both H (curl) and H (div), over each subdomain (and therefore the entire domain), we have a linear combination of shape functions in L 2 that satisfy, depending on the desired approximation space, certain conditions on their partial derivatives.As a result, this finite sum remains in L 2 .Furthermore, between each subdomain, as a result of the construction in Section II-A, the constraints on the curl and divergence, respectively, are maintained.In particular, regardless of the difference in refinement level between subdomains, tangential continuity or normal continuity is enforced for seeking solutions in H (curl) and H (div), respectively.As a result, the theoretical provisions of, for example, [27], [29], for the efficacy of hp-refinements in many computational electromagnetic applications of interest is unchanged as a result of this streamlined RBS structure for H (curl)-and H (div)-conforming approximations.Potential convergence rates with respect to the NDoFs are, therefore, still limited by the discretization properties (both local spatial resolution and local function space fidelity) and the behavior of the underlying, unknown exact solution.Moreover, as in the case of RBR, effective exploitation of the RBS infrastructure depends on error estimation and adaptivity to assemble the approximation space for solution objectives or requirements.

D. Handling Integration in the Multilevel Structure
With the approximation structure complete, we now outline the procedure for handling the filling of the Galerkin system matrices.Recall from Section II and (8), that the RBS approach, as opposed to the more typical RBR perspective, reduces sparsity under both hand p-refinements, as opposed to just p-refinements.This reduction of sparsity directly results from the integration requirements demanded by RBS-type infrastructure.
As in the organization of the DoFs relative to cell neighbors for enforcing the constraints in Sections II-A and II-B, integration procedures benefit from an inheritance tree.Based on this inheritance tree and the difference in levels between cells of overlapping subdomains, integrations are performed only over cell intersections.For testing and trial shape functions associated with the same cell, numerical integration is entirely unchanged from the conventional case.When the difference in levels is nonzero, then one cell is an active descendant that is a strict subset of an active ancestor cell; in such case, (mechanically) the integration is performed over the entire subdomain of the descendant cell and the equivalent region on the ancestor cell.Conceptually, of course, the integration may be treated as an integration over the full ancestor cell, noting that the shape function associated with the descendant cell is nonzero only over its subdomain; however, this interpretation has little practical value considering the impact on quadrature.Specifically, this cell intersection driven integration perspective preserves the efficacy (and in polynomial integrand cases, a priori convergence rates) of conventional quadrature rules (typically Gauss-Legendre, but not exclusively) experienced with uniform or RBR meshes.
In comparison to conventional matrix filling, therefore, each interaction between trial and testing shape functions relies on query to the inheritance tree for the active cells, with the integrals appropriately constrained by the overlapping subdomains.

III. NUMERICAL RESULTS
We now demonstrate the suitability of the RBS hprefinement methodology by solving the following Maxwell eigenvalue problem (in variational form): for a finite dimensional subspace B hp ⊂ H (curl; ), where m(u hp , φ hp ) = ⟨u hp , φ hp ⟩, and a(u hp , φ hp ) = ⟨∇ t × u hp , ∇ t × φ hp ⟩.We further assert that the domain ⊂ R 2 is terminated by the Dirichlet boundary condition n × u hp = 0 on ∂ .Finally, u hp is purely transversal (meaning that node-type DoFs do not appear).
While not exclusively applicable to eigenvalue problems with singularities, we study the approach for a 2-D cross section of an L-shaped waveguide, shown in Fig. 2(a), which features many singular eigenfunctions, to demonstrate the capability to achieve exponential convergence in the presence of solution irregularity.We focus our analysis on the convergence of the smallest eigenvalue to an accurate numerical computation [30] of the benchmark problem originally proposed by [34].The eigenfunction associated with this eigenvalue exhibits a singularity in the field at the reentrant corner, as seen in Fig. 2(b).
Following the procedure outlined in Section II, the initial discretization is successively refined about the reentrant corner.New refinement layers are inserted in groups with p = 1 and the expansion orders of each preexisting cell are increased by one each iteration, resulting in an emphasis on h-refinements closer to the reentrant corner and an emphasis on p-refinements away from the reentrant corner.We note that this illustrative a priori refinement strategy is neither optimal nor adaptive.Adaptive strategies, such as in [30], may be applied in place of the illustrative refinement approach Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
presented in this manuscript.A collection of discretizations with h-refinements targeting the reentrant corner (from L = 0 to L = 8 refinement levels) and global (i.e., uniform) increments in p serve as the comparison approach.
Example discretizations from each approach with five refinement layers are illustrated in Fig. 3. Fig. 3(a) depicts the progression from third-order field expansion to first-order while undergoing simultaneous h-refinements and Fig. 3(b) features the same level of h-refinement with homogeneous third-order field expansion.
Fig. 4 shows the convergence results for the first eigenvalue with the two approaches to refinement.RBS hp-refinement achieves exponential convergence while the successively p-refined discretizations at various levels of h-refinement (L = 0 to L = 8) provide only algebraic convergence.The linear trend with respect to (NDoFs) 1/3 as in Fig. 4(b) indicates the strong consistency of the exponential convergence.

IV. CONCLUSION
We have demonstrated the capability to achieve exponential convergence through an RBS hp-method in CEM.At the cost of reducing sparsity in FEM applications, the significant reduction in implementation complexity facilitates straightforward adoption of hp-refinement techniques with arbitrary levels of refinement.
When applied to the computation of the eigenvalue associated with a singular eigenfunction for H (curl)-conforming elements, the method delivers perfect exponential convergence while enforcing the tangential continuity requirements by construction rather than through constraint equations.Finally, the entire procedure directly applies to enforcement of normal continuity when H (div)-conforming elements are required and also extends to 3-D applications easily.
Finally, an open-source library based on this approach was further developed in [35] and is available at [36].

Fig. 1 .
Fig. 1.RBS hp-refinement activation and deactivation procedure.(a) Depiction of the process in 1-D.(b) Depiction of the process in 2-D.(c) Overhead perspective of the distribution of h-refinements applied in the 2-D example.

Fig. 2 .
Fig. 2. Model and problem under study.(a) Initial discretization for the L-shaped domain.(b) Field magnitude of the first eigenfunction, illustrating the singularity at the reentrant corner.

Fig. 3 .
Fig. 3. Example discretizations for the RBS hp-method and the selectively h-refined comparison method with uniform p.(a) RBS hp-method discretization with maximum and minimum expansion orders of three and one, respectively.(b) RBS h-method with a uniform expansion order of three.The two discretizations have L = 5 refinement levels.

Fig. 4 .
Fig. 4. Convergence of the first eigenvalue for the RBS hp-method and the comparison approach with h-refinement levels from L = 0 to L = 8 and increasing uniform expansion orders.(a) Double logarithmic representation.(b) Log-cube-root representation.