1. 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
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 non-smooth solutions (e.g., bends, corners, or other special geometries).
A good compromise is the use of second-order basis functions [
8], as in [
9,
10]. The field of application of
h and
adaptivity, [
11,
12,
13,
14,
15,
16], is quite broad: multigroup neutron diffusion [
17], brittle fracture modeling [
18], tomography [
19], inverse scattering problems in electromagnetics [
20], electromagnetic cloaking [
21], acoustics [
22], metamaterials [
23], or mechanics [
24]. However, these techniques are specific to hexahedral and tetrahedral meshes.
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.
2. 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.
2.1. 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
as the field to approximate. We decompose the boundary of the domain
into the sets
,
, and
, corresponding to Dirichlet, Neumann, and Cauchy boundary conditions, respectively. We set
,
, and
as the vacuum wavenumber, relative electric permittivity, and relative magnetic permeability, respectively. The vector
is the outward unit normal vector regarding
. We assume isotropic and homogeneous materials without loss of generality. The terms
,
,
, and
are used for the volumetric excitation and non-homogeneous boundary conditions. If we are not using MMS,
,
, and
are set to
.
We use the space of functions [
1],.
to approximate the electric field due to its nature. The space
is the space of square-integrable functions over
. We also use the subspace
, 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:
Find
such that
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 , assuming that is a differential operator and is a source term, the first step is to manufacture employing the introduction of an analytic solution . Afterwards, we solve the differential equation obtaining the approximate solution that can be compared to 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.
2.2. 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 residual
, defined through
and a
face residual
that changes if a boundary condition is present in that face. Thus, we define the residuals
,
, and
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
for faces that belong to both elements
m and
n, i.e.,
so
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
is
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.
2.3. 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
) [
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
the set of elements of the original domain
, whereas
is the set of elements marked to be refined, with
. 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
after a first uniform refinement.
Algorithm 1 Marking strategy using the next-step estimator inspired by [38] |
- 1:
proceduremark_nextstep(,,) - 2:
← 1 - 3:
i← 0 - 4:
if then - 5:
← - 6:
else - 7:
for all do - 8:
← - 9:
end for - 10:
end if - 11:
while do - 12:
← - 13:
← - 14:
for all do - 15:
if then - 16:
- 17:
end if - 18:
end for - 19:
i← - 20:
end while - 21:
← - 22:
end procedure
| ▹ Estimator of the last step ▹ Process of conformation to avoid hanging nodes ▹ Set of original elements (of size N) ▹ Set of elements marked to be refined (of size )
▹ Number of iterations of mark_nextstep ▹ To force a first round of uniform refinement ▹ Go through each element ▹ Adjustment of ▹ Threshold to be marked ▹ Add to the set of marked elements ▹ Increase the iteration number ▹ Conform marked elements
|
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 introduced in the fixed-energy fraction and SER algorithms, included in Algorithms 4 and 5, respectively.
Algorithm 2 Quantile marking strategy based on [39] |
- 1:
proceduremark_quantile(,,,) - 2:
← - 3:
← - 4:
while do - 5:
for all do - 6:
- 7:
end for - 8:
← - 9:
← - 10:
end while - 11:
end procedure
| ▹ Process of conformation to avoid hanging nodes ▹ Order elements from highest to lowest local residual ▹ Percentage of elements to be refined ▹ Set of original elements (of size N) ▹ Set of elements marked to be refined (of size NR)
▹ Get an ordered set of elements ▹ Add to the set of marked elements ▹ Conform marked elements ▹ Try with a smaller set
|
Algorithm 3 Maximum marking strategy |
- 1:
proceduremark_maximum(,,,) - 2:
← - 3:
for all do - 4:
if then - 5:
- 6:
end if - 7:
end for - 8:
← - 9:
end procedure
| ▹ Process of conformation to avoid hanging nodes ▹ Percentage to set threshold ▹ Set of original elements (of size N) ▹ Set of elements marked to be refined (of size )
|
Algorithm 4 Fixed-energy fraction marking strategy |
- 1:
proceduremark_fixedenergy(,,,) - 2:
← - 3:
← - 4:
← 1 - 5:
while do - 6:
for all do - 7:
- 8:
end for - 9:
← - 10:
← - 11:
end while - 12:
end procedure
| ▹ Process of conformation to avoid hanging nodes ▹ Order elements from highest to lowest ▹ Percentage to set amount of energy of ▹ Set of original elements (of size N) ▹ Set of elements marked to be refined (of size ) ▹ Threshold to be marked in energy ▹ Get an ordered set of elements ▹ Initialize number of elements to mark. ▹ is the size of ▹ Add to the set of marked elements ▹ Conform marked elements ▹ Increase the number of elements to refine.
|
Algorithm 5 Slight variation of the marking strategy in [42] |
- 1:
proceduremark_ser(,,,) - 2:
← - 3:
← - 4:
← 1 - 5:
← 1 - 6:
while do - 7:
← - 8:
for all do - 9:
- 10:
end for - 11:
← - 12:
end while - 13:
end procedure
| ▹ Process of conformation to avoid hanging nodes ▹ Order elements from highest to lowest ▹ Percentage to set amount of energy of ▹ Set of original elements (of size N) ▹ Set of elements marked to be refined (of size )
▹ Threshold to be marked in energy ▹ Get an ordered set of elements ▹ Initialize number of elements to mark. ▹ Initialize threshold to mark elements ▹ is the size of ▹ Decrease threshold in each iteration ▹ Add to the set of marked elements ▹ Conform marked elements
|
2.4. 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 1
x elements as 2
x 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 1
x, we check the divided angle and we mark it as 2
x 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.
4. 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.