Adaptive Semi-Structured Mesh Reﬁnement Techniques for the Finite Element Method

Featured Application: This work can improve the computational resources employed to achieve a given accuracy in electromagnetic simulation. Abstract: The adaptive mesh techniques applied to the Finite Element Method have continuously been an active research line. However, these techniques are usually applied to tetrahedra. Here, we use the triangular prismatic element as the discretization shape for a Finite Element Method code with adaptivity. The adaptive process consists of three steps: error estimation, marking, and reﬁnement. We adapt techniques already applied for other shapes to the triangular prisms, showing the differences here in detail. We use ﬁve different marking strategies, comparing the results obtained with different parameters. We adapt these strategies to a conformation process necessary to avoid hanging nodes in the resulting mesh. We have also applied two special rules to ensure the quality of the reﬁned mesh. We show the effect of these rules with the Method of Manufactured Solutions and numerical results to validate the implementation introduced. Author Contributions: Conceptualization, A.A.-M. and L.E.G.-C.; software, A.A.-M.; validation, A.A.-M. and L.E.G.-C.; investigation, A.A.-M. and L.E.G.-C.; data curation, A.A.-M.; writing—original preparation, A.A.-M.; writing—review and editing, L.E.G.-C.; supervision, L.E.G.-C.; resources,


Introduction
The Finite Element Method (FEM) is a mature tool used to obtain the numerical solution of partial differential equations (PDEs) used in multiple engineering fields and physics [1][2][3][4]. Apart from the solid mathematical foundation of FEM, one of the main advantages of this technique is the use of adaptive refinement to improve the error convergence in the approximation of the field. This leads to better use of the computational resources since we obtain automatically more accurate solutions with fewer unknowns to solve. Additionally, the use of these techniques provides an estimation of the error in the problem. In FEM, we can perform h refinements (decreasing the size of the elements used for discretizing the original domain) and p refinements (increasing the order of basis functions that approximate the field under study). The combination of these two refinements leads to the so-called hp refinement [5][6][7], which might lead to exponential convergence when appropriate estimators are used.
The p refinements are especially effective when approximating smooth solutions, whereas h (also called mesh adaptive) refinements are more general and adapt to nonsmooth solutions (e.g., bends, corners, or other special geometries).
The use of different shapes depending on the geometry of the problem is advantageous from the computational point of view. We should use tetrahedra for irregular and unstructured geometries and hexahedra for rectangular-like geometries. The triangular prism can be considered as a hybrid of these two, and it can be used to connect hexahedral and tetrahedral meshes [25], whereas it is also well suited for planar geometries with an irregular pattern on the surface, e.g., laminates and sandwich structures [26,27] integrated circuit package designs [28], planar microwave circuits and antennas [29], waveguide sections [30], or coating for irregular volumes meshed with tetrahedra. Additionally, the generation of volumetric meshes for prisms is typically easy since we only need to extrude a 2D triangular mesh.
We need to introduce adaptivity techniques to all the different shapes to choose the most suitable shape for each problem. Unfortunately, the application of adaptivity techniques to triangular prisms is not so frequent in the literature apart from [31,32]. Whereas the formulation and the definition of the estimator do not depend on the discretization shape, the development of marking strategies and the application of different h refinements need to be particularized for the different shapes.
Here, we suggest to take the fundamental ideas from 2D refinement and apply them rigorously to the adaptivity with triangular prisms, solving the particularities (especially related to maintaining admissible meshes in the sense of having conformal solutions from the FEM point of view straightforwardly) that arise from the application of these ideas to volumetric meshes. First, we use the Method of Manufactured Solutions, MMS [33], to show the significance of the mesh quality indicators we have introduced and to experiment with five different marking strategies of the elements to be refined. Then, we validate our implementation with a propagation problem (specifically, a WR-90 rectangular waveguide). Finally, we show with an L-shaped waveguide section the effectiveness of the adaptivity algorithm.
The rest of the paper is structured as follows: In Section 2, we detail each different step of the adaptivity algorithm. In Section 3, we provide meaningful experiments with MMS and propagation problems to show the performance of the adaptivity refinement developed in the previous section. Finally, we draw the conclusions that can be extracted from this work in Section 4.

Methods
Every adaptive algorithm needs four steps: • Solve the electromagnetic domain under study. • Estimate how accurate the approximation of the fields is. • Mark the elements we want to refine to improve the accuracy. • Refine the defined elements in the previous step.
We repeat these four steps until we achieve a threshold of the estimated error or after a number of iterations is run. This section explores each of these steps separately.

Variational Formulation
We assume that the domain Ω of the electromagnetic problem to be solved is a smooth domain for which we have the following boundary value problem: taking the electric field E as the field to approximate. We decompose the boundary of the domain ∂Ω into the sets Γ D , Γ N , and Γ C , corresponding to Dirichlet, Neumann, and Cauchy boundary conditions, respectively. We set k 0 , ε r , and µ r as the vacuum wavenumber, relative electric permittivity, and relative magnetic permeability, respectively. The vectorn is the outward unit normal vector regarding Ω. We assume isotropic and homogeneous materials without loss of generality. The terms F, Ψ D , Ψ N , and Ψ C are used for the volumetric excitation and non-homogeneous boundary conditions. If we are not using MMS, F, Ψ D , and Ψ N are set to 0.
We use the space of functions [1],.
to approximate the electric field due to its nature. The space L 2 (Ω) is the space of squareintegrable functions over Ω. We also use the subspace H 0 (curl, Ω), whose functions enforce Dirichlet boundary conditions. We define the sesquilinear forms: using * to denote complex conjugation, Ω as a given volume, and Γ as a specific surface. Applying the Galerkin method yields: The solution to this problem is obtained by solving the system of equations obtained after the discretization of the domain into second-order triangular prisms from [10].

Method of Manufactured Solutions
The MMS consists of the manufacture of an analytical solution to a differential equation by solving the problem backward. If we take an equation D(E) = f , assuming that D is a differential operator and f is a source term, the first step is to manufacture f MMS employing the introduction of an analytic solution E MMS . Afterwards, we solve the differential equation obtaining the approximate solution E FEM that can be compared to E MMS to assess the quality of the approximation.
The excitations that we use in (7) are We use MMS in Section 3.1 to show the effect of the different marking strategies and the impact of the rules introduced in the algorithm to improve the quality of the adapted mesh.

Estimator
For computational purposes, the estimator of the areas where the accuracy of the solution would benefit more to have smaller elements should be cheap to compute. We use a straightforward local estimator based on [34], although other good alternatives are [35][36][37]. We compute the estimator individually for each element m with a volumetric and a face residual R (m) face that changes if a boundary condition is present in that face. Thus, we define the residuals R D , R N , and R C for the Dirichlet, Neumann, and Cauchy boundary conditions, respectively, as We also use the tangential continuity of the magnetic field (not imposed explicitly) as the face residual R neigh for faces that belong to both elements m and n, i.e., neigh for each face in the m-th element. Adding up all the residuals, we introduce the residual for each element m, which is introduced through where K k is 1 2 and 1 for triangular and rectangular faces, respectively. Moreover, h is the diameter of the entity (volume or face). Note that the sesquilinear forms are local to the element m or a given face k.

Marking Strategies
We use five different marking strategies, i.e., algorithms that decide which elements are refined in terms of the local residual defined in Section 2.2: • Next-step estimator: inspired by [38], it is based on the estimation of the error that an element would have when refined. • Quantile: here, we order the elements from highest to lowest residual and we refine a given percentage, as shown in [39]. • Maximum: we define a threshold relative to the maximum residual in the mesh (using a parameter θ ∈ [0, 1]) [40]. • Fixed-energy fraction: we order the elements by their residual from highest to lowest, and we refine the first elements which constitutes a given percentage of the energy of the residual in the whole problem [41]. • SER (Solve-Estimate-Refine): we combine here the fixed-energy fraction and maximum strategies [42]. We mark sets of elements (ordered from highest to lowest residuals) which are smaller and smaller (determined by a varying Heaviside step function) to achieve a given percentage of the energy of the whole problem. If the step function is forced to select only one element, we have the fixed-energy fraction strategy.
All of these strategies have been tested with unstructured meshes, i.e., with triangles or tetrahedra. We need to introduce modifications to these techniques for them to be applied to semi-structured meshes. In particular, we use a process of conformation to avoid the well-known issue of hanging nodes [6]. We have a hanging node when a direct correspondence of basis functions cannot be established between neighbor elements: e.g., when the edge of one element starts in the middle of the edge of the neighbor element. We could use special strategies to remedy this problem [6], but they are out of the scope of this work. This conformation process makes undetermined the number of prisms to be refined in the marking step, since the refinement has to be propagated horizontally and vertically.
We call P the set of elements of the original domain P, whereas P R is the set of elements marked to be refined, with P R ⊂ P. In Algorithm 1, we suggest a variation (introducing an adaptive threshold to guarantee that, in each round, some elements are marked) of the algorithm proposed in [38], assuming that the residual follows asymptotically R = ch λ after a first uniform refinement.

Algorithm 1
Marking strategy using the next-step estimator inspired by [38] Require: R for all e ∈ P do Go through each element 8: end for 10: end if 11: while N R = 0 do 12: θ ← θ/10 i Adjustment of ς cut 13: k+1 } Threshold to be marked 14: for all e ∈ P do 15: if R (e) k ≥ ς cut then 16: e ∈ P R Add to the set of marked elements 17: end if 18: end for 19: Increase the iteration number 20: end while 21: We include the algorithm for the quantile marking strategy in Algorithm 2. We adapt the implementation to the presence of the conformation process: we perform an iteration considering fewer and fewer elements until the number of elements to be refined is lower than the given percentage (specified through parameter θ) of elements to be refined.
The maximum marking strategy does not require special treatment due to the conformation of the mesh, as shown in Algorithm 3, in contrast to the slight modifications intro-duced in the fixed-energy fraction and SER algorithms, included in Algorithms 4 and 5, respectively.

Algorithm 2
Quantile marking strategy based on [39] Require: conform_refinement(P R ) Process of conformation to avoid hanging nodes Require: order_residual(R k ) Order elements from highest to lowest local residual R Get an ordered set of elements 4: while N R < θN do 5: for all e ≤ θ R · N R do 6: e ∈ P R Add to the set of marked elements 7: end for 8: Threshold to be marked in energy 3: Get an ordered set of elements 4: Initialize number of elements to mark. 5: while for all e ≤ N ref do 7: e ∈ P R Add to the set of marked elements 8: end for 9: P R ← conform_refinement(P R ) Conform marked elements 10: Increase the number of elements to refine.
Threshold to be marked in energy 3: Get an ordered set of elements 4: Initialize number of elements to mark.

5:
θ R ← 1 Initialize threshold to mark elements 6: while ∑ Decrease threshold in each iteration 8: for all e ≤ N ref do 9: e ∈ P R Add to the set of marked elements 10: end for 11: P R ← conform_refinement(P R ) Conform marked elements 12: end while 13: end procedure

Refinement
Our refinement strategy uses two-dimensional strategies for triangles (horizontal direction) applying afterwards an extrusion in the segment (vertical) direction, being split if required. To distinguish between the different kinds of refinements that appear in the algorithm, we use a two-digit number, being the first digit related to the horizontal direction, whereas the second digit includes the information for the vertical direction.
For the horizontal part, we use a red-green algorithm, [43]. When we mark a triangle to be refined, we mark it as. red, so four new triangles are created joining the middle points of each edge. We set the neighbor triangles (when they are not marked themselves) to green, i.e., we create two new triangles joining the middle point of the red edge to the opposite vertex. The triangles with more than one red neighbor are automatically marked also as red (although other variations, as the red-green-blue algorithm, are possible [40]). We denote red and green refinements by 2 and 1, respectively. For the vertical direction, we note if there is a refinement with a 1, leading to five different refinements as shown in Figure 1.
To avoid the problem of hanging nodes, we apply a conformation process that regularizes the resulting mesh. We need to propagate the refined elements in the horizontal and vertical direction, i.e., when we refine the triangular faces of a given prism, we have to refine the triangular faces of the top and bottom neighbors. Let us assume that only one element is marked as 21 (refinement in both directions), as in Figure 2. Then, the prisms in the same layer are marked as 11, since we mark the neighbors triangles as green and we need to refine the vertical direction to avoid hanging nodes. For the same reason, all the prisms in the same layer without triangle refinement are marked as 1, whereas the vertical neighbors of the 21 elements are marked as 20, and the horizontal neighbors of these new 20 elements are marked as 10.  However, the use of green refinement might lead to meshes of poor quality, which might introduce additional error [44]. One problem is the appearance of eyes in the mesh if green elements are used in one vertex shared by many elements (and the vertex whose angle is divided is this one vertex). This is a non-reversible problem since the error would accumulate in this point, and we cannot improve the shape of the triangles introducing red elements. We use an accumulator for all the vertices in the mesh that counts the number of elements that have this vertex, setting a threshold and marking 1x elements as 2x for these eyes not to appear. The other problem is the appearance of very deformed triangles (i.e., with very small angles). To remedy that, when an element is marked as 1x, we check the divided angle and we mark it as 2x when it is lower than 23 degrees. (which experimentally has been proven as a good threshold). We show the effectiveness of these two quality criteria in Section 3.1.

Results
We use three sets of results: first, we use MMS to assess the effect of the different marking strategies and the two criteria for ensuring the quality of the mesh; second, we employ a WR-90 waveguide to validate the refinement strategy; and, third, we analyze an L-shaped waveguide to show the effectiveness of the whole algorithm.

MMS
We solve a cube with dimensions [0, 1] m × [0, 1] m × [0, 1] m with vacuum as material. We use a working frequency of 50 MHz. We use to manufacture the excitation and solve the problem. This choice is motivated by the fact of having a polynomial out of the second-order space of basis functions which produces a higher error with a higher value of y. We show a mesh with an eye in Figure 3, whereas, in Figure 4, we have avoided this eye with the procedure detailed in Section 2.4. We can see an abnormal increase in the error distribution around these eyes, therefore showing the importance of avoiding them in real applications. In these figures and all that follow, the distribution of the error is computed at each point as E FEM − E MMS 2 .    The effect of the second criterion to improve the quality of the refined mesh is present in Figure 5. If we forced that the 1x elements have a maximum angle of 0.4 rad in the problematic vertex, we increase the number of elements of the refinement but we don't obtain abnormal error areas, as shown in Figure 6. However, in this last figure, we can see another abnormal area related to the appearance of non-treated eyes. Thus, if we also apply the criterion present in Figure 4, we can see that the error is higher where it is expected as in Figure 7.      Finally, we assess the effect of the five marking strategies detailed in Section 2.3 with the parameters shown in Table 1. We apply different rounds of refinement and compare the maximum value of the error in Figure 8. The difference between the marking strategies is small, although all of the strategies outperform the uniform refinement. In addition, the best strategy with this metric is Algorithm 3 as expected since we are using the maximum value of the point error.  To provide a better insight into the effect of each strategy, the meshes and the error for the five marking strategies after three rounds of refinement are shown in Figures 9-13. The differences between strategies are small and strongly depend on the chosen parameter. The strategies based on the energy of the residual (fixed-energy and SER) show smaller elements especially around the areas of stronger error (even when the variation between residuals is minimal), whereas the maximum strategy is more uniform. On the other hand, the quantile strategy might lead to some irregular meshes (since we may apply the same refinement to elements with very different residuals). Finally, the next-step marking strategy needs to refine all the mesh in the first step, leading to more structured meshes.   Figure 13. Approximation of (13) using the SER marking strategy.

WR-90 Waveguide
To show the performance of the algorithm with real problems, we simulate a WR-90 waveguide with length l = 1λ, using a working frequency of f = 7.5 GHz. We use a structured mesh to show clearly the effect of the refinement, and the maximum marking strategy with θ = 0.65. We show the initial mesh in Figure 14. The estimator shows that the higher error is located close to the electric walls, so the marked elements are in that area. The first round of refinement is included in Figure 15, where we can see that the error is lower in the refined elements. Thus, the main contribution to the error is located in the middle of the waveguide, where we mark the elements to be refined. Finally, the same effect is seen in Figure 16, noting that. we have achieved the behavior of the algorithm that we would expect.

L-Shaped Waveguide
To show the utility of an adaptivity algorithm, we solve an electromagnetic problem with a singularity. We bend a WR-90 waveguide along its transverse face (E-plane), where the electric field is maximum. We show the geometry of the problem in Figure 17, using a working frequency of f = 7.5 GHz. For the adaptivity refinement, we use the maximum marking strategy with θ = 0.65. We show the initial mesh and the results of applying the estimator and the marking strategy in Figure 18. In addition, in Figures 19 and 20, we show the same. outputs for two additional steps of refinement. We see that the estimator yields the highest error around the singularity, with a closer and closer location of the red area with smaller and smaller elements.    We show the evolution of the electric field in Figure 21, whereas in Figure 22 we use only uniform refinement. We can see that the singularity is much better represented when using the adaptivity algorithm in contrast with uniform refinement. On the one hand, the value of the field at the singularity is higher when using the adaptivity algorithm. In this sense, for quantitative comparison, we show in Figure 23 the evolution of the value of the singularity with successive levels of refinement for both strategies (adaptive as in Figure 21 and uniform as in Figure 22) with respect to the number of unknowns. We observe that the singularity is better represented (since it has a higher value) with fewer unknowns than the uniform strategy, as we should expect from the previous results. On the other hand, the singularity of the field is constrained to the corner when the adaptivity algorithm is used, i.e., the adaptivity improves the representation of the strong decaying of the field out of the singularity.

Conclusions
There exist many different possibilities in the FEM community when choosing the shape of the finite element. For versatility and convenience, specific techniques are sometimes developed only for tetrahedra or hexahedra. To effectively choose the discretization shape in terms of the geometry of the problem, we need to develop the same techniques for all the shapes.
Specifically, the use of adaptive techniques is very popular to save computational resources and introduce smart meshes where the areas that contribute more to the error are further refined. However, the use of these techniques applied to triangular prisms is not so common in the literature, whereas the use of these elements is quite advantageous in some scenarios (e.g., in planar circuits or thin material layers). For this reason, we have developed an adaptivity technique specifically designed for triangular prisms.
We can summarize the main contributions of this paper with the following points: • We have introduced an estimator inspired by [34], with minor modifications to adjust the use of triangular prisms. • We have experimented with five different marking strategies, modified to be adapted to the semi-structured nature of the triangular prism. • We have suggested a refinement technique based on the red-green algorithm from [43]. • We have introduced specific criteria to increase the quality of the generated meshes: specifically, we avoid the eyes in the mesh and the appearance of highly deformed triangles. • We have developed a conformation process to avoid the problem of the hanging nodes.
We have validated with numerical results the performance of the algorithm. We have shown the different behavior of each marking strategy: strategies based on the energy of the residual for each element show smaller elements around areas of stronger error, whereas strategies based on a given threshold of the element with maximum error yield more uniform meshes. In contrast, strategies based on refining a given percentage of the elements lead to more irregular meshes. All of the strategies provide smaller elements in higher error areas for both MMS and propagation problems. In addition, we have shown the performance of the algorithm with a singularity due to the geometry of the problem to be solved. This algorithm provides smaller elements in the vicinity of the singularity, as expected giving a better approximation of the singularity with fewer unknowns compared to the uniform refinement.
As future research lines, the treatment of hanging nodes would eliminate the process of conformation, easing the step of refinement. Moreover, we could detect in which direction the estimator is more needed and apply the refinement only in the horizontal or vertical direction as detected. Funding: This work has been financially supported by TEC2016-80386-P.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

FEM
Finite Element Method MMS Method of Manufactured Solutions SER Solve-Estimate-Refine 2D Two dimensions