A High Resolution PDE Approach to Quadrilateral Mesh Generation

TL;DR

This paper introduces a high-order spectral element method for generating high-resolution quadrilateral meshes in complex 2D domains, leveraging spectral convergence and advanced field-guided techniques.

Contribution

It presents a novel high-order PDE-based approach using spectral elements and open source tools for precise quadrilateral mesh generation without post-processing.

Findings

Spectral convergence achieves high accuracy in guiding fields.

Corner handling with DG allows flexible discretization.

Meshes have naturally curved quadrilaterals without additional curvature adjustments.

Abstract

We describe a high order technique to generate quadrilateral decompositions and meshes for complex two dimensional domains using spectral elements in a field guided procedure. Inspired by cross field methods, we never actually compute crosses. Instead, we compute a high order accurate guiding field using a continuous Galerkin (CG) or discontinuous Galerkin (DG) spectral element method to solve a Laplace equation for each of the field variables using the open source code Nektar++. The spectral method provides spectral convergence and sub-element resolution of the fields. The DG approximation allows meshing of corners that are not multiples of $\pi/2$ in a discretization consistent manner, when needed. The high order field can then be exploited to accurately find irregular nodes, and can be accurately integrated using a high order separatrix integration method to avoid features like limit cycles. The result is a mesh with naturally curved quadrilateral elements that do not need to be curved a posteriori to eliminate invalid elements. The mesh generation procedure is implemented in the open source mesh generation program NekMesh.

Tables

Table 1. Table 1 : Value of I c subscript 𝐼 𝑐 I_{c} , ( 12 ), and the associated valence.

$I_c$	Valence ($𝒱$)
$-1$	5
$0$	4
$+1$	3

Equations

$$C\left(\psi \right)=\left\{{\vec{c}}_k={\left(c_x,c_y\right)}_k=\left(\cos \left(\psi +k\frac{\pi }{2}\right),\sin \left(\psi +k\frac{\pi }{2}\right)\right)\right\},$$

$$\psi =\frac{1}{4}atan2\left(\mathrm{v},\mathrm{u}\right)\in \left[-\frac{\pi }{4},\frac{\pi }{4}\right],$$

$${\vec{v}}_b=\left(u_b,v_b\right)=\left(\cos \left(4{\theta }_b\right),\sin (4{\theta }_b)\right),$$

$$\{\begin{array}{c}{\nabla }^2\vec{v}=0,\vec{x}\in \Omega \text{hfill}\\ \vec{v}={\vec{v}}_b,\vec{x}\in \partial \Omega \text{hfill}\end{array}$$

$$a(\vec{v},\vec{w})=l(\vec{w})\forall \vec{w}\in H^1(\Omega )$$

$$H^1(\Omega )=\{\vec{w}\in L^2(\Omega )\mid D^{\alpha }\vec{w}\in L^2(\Omega )\forall \left|\alpha \right|\le 1\}.$$

$$a({\vec{v}}^{\delta },{\vec{w}}^{\delta })=l({\vec{w}}^{\delta })\forall {\vec{w}}^{\delta }\in V_N,$$

$${\vec{v}}^{\delta }(\vec{x})=\sum _n{\hat{v}}_n{\varphi }_n({\mathcal{X}}_e^{-1}(\vec{x})),$$

$$V:=\{\vec{v}\in H^1(\Omega )\mid {\vec{v}|}_{{\Omega }_e}\in {\mathcal{P}}_P({\Omega }_e)\},$$

$$i_{\gamma }=\frac{1}{2\pi }{\oint }_{\gamma }d\varphi ,$$

$$\varphi =atan2\left(\mathrm{v},\mathrm{u}\right)=4\psi .$$

$$I_c=\frac{1}{\pi /2}{\oint }_c\frac{d\psi }{d\theta }d\theta ,$$

$$\mathcal{V}=4-I_c.$$

$$I\left({\theta }_0,{\theta }_f\right)=\frac{1}{\pi /2}{\int }_{{\theta }_0}^{{\theta }_f}\frac{d\psi }{d\theta }d\theta ,$$

$$\mathcal{V}=\frac{\Delta \theta }{\pi /2}-I\left({\theta }_0,{\theta }_f\right),$$

$$\frac{d\vec{x}}{dt}=\tilde{v}\left({\psi }^{\prime }(\vec{x})\right)$$

$$\vec{x}(t_0)={\vec{x}}_0$$

Full text

††journal: Journal of Computational Physics

1 Introduction

A quadrilateral subdivision of a two dimensional domain into a minimal number of subdomains can serve as the starting point from which to generate meshes for finite element [1], finite volume, (block structured) finite difference [2] and spectral element methods [3, 4]. To generate finite element and finite volume meshes, one would subdivide each subdomain into quadrilateral cells sufficiently small to resolve the geometry and expected solution features. The resulting mesh would have the desired property of only a small number of irregular nodes where the valence, i.e., the number of elements sharing a point, is not equal to four. The subdivision of the subdomains could also be used as nodes for high order block-structured finite difference methods if the subdivision within each block is structured. Finally, the subdomains could be used as is with high order geometry information intact, or subdivided in a structured or unstructured manner, to be used as elements for a quadrilateral based spectral element method.

Fully automatic quadrilateral mesh generation is still difficult [5], especially for high order (curved) elements. Most previous work has focussed on generating large numbers of small straight sided elements. The extension to high order is usually done by curving the elements a posteriori, e.g. [6]. Unstructured quadrilateral mesh generators tend to generate meshes with large numbers of irregular nodes [7]. Extraneous irregular nodes are generated even for simple geometries, where mesh topologies can be generated by hand with a minimum number of them. As an example, Fig. 1 shows two meshes generated with a subdivision algorithm [8, 9], which generates more subdivisions than necessary.

One promising approach for generating quadrilateral subdivisions in arbitrary domains has come out of the computer graphics community. It uses cross fields [1, 5] in a field-guided procedure. Cross fields resemble cross hatching used in drawings. Crosses, which are composed of two direction vectors and their negatives at a point, are invariant to rotations of $\pi/2$. Cross field methods generate meshes with fewer irregular nodes. Fig. 2 shows the same two geometries as Fig. 1 decomposed with a cross field method. The method produces the same topology as an experienced user might generate by hand, and has a minimum number of irregular nodes.

Details of cross field procedures vary, but the basic idea is that the crosses are oriented, usually at the nodes of an existing triangular mesh, according to a minimization or smoothing procedure. Triangles in which the crosses cannot smoothly vary contain singularities, which serve as the irregular nodes of the quadrilateral decomposition to be computed. The internal edges of the quadrilaterals are found by integrating separatrices (represented as streamlines of the cross field) [10], starting from the irregular nodes until they reach another irregular nodes or a physical boundary. The quadrilateral subdivision is usually further subdivided into smaller quadrilateral elements.

Low order techniques using crosses do not precisely locate irregular nodes within the domain. Irregular nodes can be located only to within the size of the element in which they occur. Spurious singularities can be generated, which may have been coalesced in an ad hoc fashion. For instance, if the underlying mesh is too coarse it is possible that two or more singularities can fall within one element. The valence of the irregular points computed from a coalesced singularity could lead to the incorrect number of separatrices and hence the incorrect valence of the associated interior quadrilateral mesh node.

Low order tracing of the separatrices also often leads to multiple curves that need to be coalesced. Associated with this is the well-known limit cycle problem where separatrices, instead of meeting as they should, spiral indefinitely [11].

Finally, boundary conditions are difficult to apply at corners where a jump discontinuity in the cross field will occur when a standard continuous Galerkin (finite element) approximation is used for smoothing. Rather than putting a singularity there, the corner may be effectively smoothed by averaging the cross field at such points [1], at the expense of adding more interior singular points.

We address issues found in the use of low order techniques by using a high resolution approach to compute a guiding field, accurately locating the interior irregular nodes, and accurately integrating the streamlines used to subdivide the domain. We use a high order continuous or discontinuous Galerkin spectral element method (SEM or DGSEM) on a triangular mesh, the approximation depending on the regularity of the boundary conditions needed to solve Laplace problems. Unlike previous work, we perform all operations on the original highly-resolved guiding field rather than a cross field so that accuracy is not lost. We use an accurate locator for the irregular points that exploits the high resolution solution, and an accurate method for computing their valence. The streamline integration then uses the high order information available to obtain precise streamlines, which reduces the need to coalesce extraneous lines.

2 Mathematical Formulation

The problem is to subdivide any two dimensional domain, $\Omega$, that has piecewise smooth boundaries into quadrilateral subdomains. The domain might be simply or multiply connected. The decomposition will be regular if all the corners of the quadrilaterals have valence four, that is, each corner node connects four edges. In general, it will not be possible to generate a purely regular mesh, especially in multiply connected domains. Instead, some of the nodes will be irregular, where the valence will be larger or smaller than four, like those observed in Fig. 3.

2.1 Crosses and Cross Field

To motivate the high resolution approach to the quadrilateral decomposition of a domain, $\Omega$, and to put it into context, we recall the notion of a cross. A cross can be represented by

[TABLE]

for, $k=0,1,2,3$. Here, the tangent angle, or phase, $\psi$, is computed from the four quadrant inverse tangent,

[TABLE]

so that a cross is represented by four unit vectors $\vec{c}_{k}$ at any point $\vec{x}\in\Omega$ given a vector field $\vec{v}=\left(u\left(\vec{x}\right),v\left(\vec{x}\right)\right)$ at which $\vec{v}\neq 0$ from the principal direction. See Fig. 4. The important property of a cross is its 4-way rotational symmetry. Rotation by an angle $\pi/2$ does not change the cross.

The function $\psi$ is tangent or orthogonal (due to the jump in the atan2 function) to the streamlines of the field $\vec{v}$. Hence, it is parallel to one of the branches of a cross. It is undefined at $\vec{v}=(0,0)$. Since $\psi$ is computed from the arctangent, there will be a jump of of value $\pi/2$ in $\psi$ depending on the signs of $u$ and $v$ even if $\vec{v}$ is smooth. We call lines in the field across which $\psi$ jumps jump lines.

A full discussion of cross fields and methods based on them is beyond the scope of this paper. See, for example, references [1, 11, 12, 13] for more comprehensive discussions.

2.2 The Guiding Field, $\vec{v}$

It is necessary, then, to find the guiding field $\vec{v}$ from which to find the irregular nodes and trace separatrices. The only values that can be set a priori are on the boundary $\partial\Omega$ of $\Omega$. We align the field with the boundary to ensure that the mesh is aligned there. Therefore, we set

[TABLE]

where $\theta_{b}$ is the tangent angle of the boundary where the field is being computed. The vector $\vec{v}_{b}$ at any point along the boundary defines a cross (1) at that point. The factor of four in eq. (3) ensures the same $u,v$ values for each $90^{\circ}$ rotation of the angle, and hence the rotational symmetry. Fig. 5 shows example crosses on the boundary of a quarter circle domain.

We make three observations from eq. (3) about the boundary field $\vec{v}_{b}$:

Along any smooth closed boundary curve (interior or exterior), there are eight zeros each for $u$ and $v$, and those zeros are not coincident.

• 2.

Along smooth portions of the boundary curves and at corners whose angle is a multiple of $\pi/2$, $u$ and $v$ are continuous.

• 3.

At boundary corners with angles not divisible by $\pi/2$ there is a jump discontinuity in $\vec{v}_{b}$ and the cross field.

Once the boundary values are specified using expression (3), they are smoothly propagated to the interior. Smoothing procedures that have been used in the past include marching the associated crosses to the interior [14], minimizing an energy functional [1], or solving a Laplace problem [15]. Solving a Laplace problem for $u$ and $v$

[TABLE]

guarantees that the field $\vec{v}$ is smooth in the interior of the domain and satisfies the maximum principle. The solution of the Laplace problem is equivalent to the minimization of the Dirichlet energy. We choose to solve (4) because it allows us to use a high resolution approach such as the spectral element methods. In particular, the use of a discontinuous Galerkin approximation allows us to handle geometries where the corner angles generate jumps in the guiding field, see Sec. 2.2.1.

2.2.1 Continuous and Discontinuous Galerkin Spectral Element Methods

A key difference between our work here and that of previous work is that we compute high resolution solutions for the BVPs (4) with either a continuous Galerkin (CG) or a discontinuous Galerkin (DG) spectral element method on a triangular mesh [3]. Spectral element methods are spectrally accurate, meaning that the convergence rate depends only on the smoothness of the solution. They are high resolution in that they use a large number of degrees of freedom within an element. Unlike traditional low order finite element methods, the high order polynomial expansion of the solution inside each element allows us to locate irregular nodes, identify their valence, and finally trace separatrices with high order accuracy.

In this work, we use the spectral/hp element methods formulation described in detail in reference [3] and implemented in Nektar++ [16, 17]. We briefly describe the fundamentals of the methods in what follows. We consider the numerical solution of PDEs of the form $\mathcal{L}u=0$ over a domain $\Omega$. The domain $\Omega$ is taken as a set of finite elements, $\Omega_{e}$ — the mesh — such that $\Omega=\cup\Omega_{e}$ and $\Omega_{e_{1}}\cap\Omega_{e_{2}}=\partial\Omega_{e_{1}e_{2}}$ is either an empty set or the interface between two elements and is of one dimension less than the mesh. We solve the PDE problem in the weak sense and require that $\left.u\right|_{\Omega_{e}}$ is in the Sobolev space $H^{1}(\Omega_{e})$. In a CG formulation, we also require the solution to be continuous from one element to the other.

We formulate the solution to (4) in weak form: $find\ \vec{v}\in H^{1}(\Omega)\ such\ that$

[TABLE]

where $a(\cdot,\cdot)$ is a symmetric bilinear form, $l(\cdot)$ is a linear form and $H^{1}(\Omega)$ is formally defined as

[TABLE]

We want to solve this problem numerically and therefore consider solutions in a finite dimensional subspace $V_{N}\subset H^{1}(\Omega)$ and cast the problem: $find\ \vec{v}^{\delta}\in V_{N}\ such\ that$

[TABLE]

augmented with boundary conditions, (3). In the CG formulation, we also enforce the condition $V_{N}\subset C^{0}$.

We take a weighted sum of $N$ trial functions $\Phi_{n}(\vec{x})$ defined on $\Omega$ so we can represent $\vec{v}^{\delta}(\vec{x})=\sum_{n}\hat{v}_{n}\Phi_{n}(\vec{x})$. This transforms the problem to that of finding the coefficients $\hat{u}_{n}$ that define $\vec{v}^{\delta}(\vec{x})$ within an element. To obtain a unique choice of coefficients $\hat{u}_{n}$, we place a restriction on $R=\mathcal{L}\vec{v}^{\delta}$ that its $L^{2}$ inner product, with respect to the test functions $\Psi_{n}(\vec{x})$, is zero. In the Galerkin projection one chooses the test functions to be the same as the trial functions, i.e. $\Psi_{n}=\Phi_{n}$.

The contributions of each element in the domain must be taken into account to construct the global basis $\Phi_{n}$. A parametric mapping $\mathcal{X}_{e}:\Omega_{e}\to\mathcal{E}$ exists from each element $\Omega_{e}$ to a standard reference element $\mathcal{E}\subset{[-1,1]}^{d}$. This mapping is given by $\vec{x}=\vec{\mathcal{X}}_{e}\left(\vec{\xi}\right)$. It is important to distinguish $\vec{x}$, the physical coordinates, from the $\vec{\xi}$, the coordinates in the reference space.

We construct a local polynomial basis on the reference element to represent the solution. A one-dimensional order-$P$ basis is a set of polynomials $\Phi_{p}(\xi),0\leq p\leq P$, defined on the reference segment $-1\leq\xi_{1}\leq 1$. In two and three dimensional reference regions, a tensor basis is used, where the polynomial space is constructed as the tensor product of one dimensional bases on segments, quadrilaterals, and hexahedral reference regions.

In Nektar++, triangular, tetrahedral, prismatic, and pyramidal elements are created by collapsing one or more of the coordinate directions to create singular vertices. This allows it to support, for this work, easy-to-generate triangular meshes.

Finally, the discrete solution in a physical element $\Omega_{e}$ can be expressed as

[TABLE]

with $\hat{v}_{n}$ the coefficients computed by the Galerkin procedure. We restrict our solution space to

[TABLE]

where $\mathcal{P}_{P}(\Omega_{e})$ is the space of order $P$ polynomials on $\Omega_{e}$.

An assembly operator can then be designed to assemble the element contributions to the global solution. In the CG formulation, elemental contributions of neighbours are summed to enforce $C^{0}$-continuity. In the DG formulation, such mappings transfer flux values from the element interfaces into the global solution vector.

2.3 Critical Points

Points where $v=u=0$ are called critical points in the guiding field and singular points in the cross field. They are points where $\psi$ cannot be determined uniquely. Such points become the interior irregular nodes (valence $\neq 4$) in the quadrilateral decomposition.

The form of the boundary values usually implies that critical points will exist in the interior of the domain. Since there are interlaced zeros in $u$ and $v$ on the boundary, such as occur when the boundary is smooth and $\theta_{b}$ varies continuously, the fact that the interior field is smooth implies that zero contours of $u$ and $v$ must cross somewhere in the interior, thereby creating a critical point. Zeros in $u$ will exist along the boundary at any point where $\theta_{b}$ smoothly passes through $\theta_{b}=k\pi/8$, $k$ odd, and in $v$ for $k$ even. Therefore, one can infer the existence of interior irregular nodes in the final quadrilateral mesh from the boundary curves. A semi-circular section of a boundary, for instance will have four zeros in $u_{b}$ and five in $v_{b}$ and, in the absence of other nearby features, will create two interior critical points. See Fig. 6. We note that although the cross field directions are undefined at critical points, the field $\vec{v}$ is smooth because the solutions of the Laplace equation (4) are regular.

Irregular nodes in the mesh occur as recognizable critical point topologies in the $\vec{v}$ field [18]. They can be categorized analytically by the Poincaré index

[TABLE]

where

[TABLE]

When $\gamma$ encloses a critical point, $i_{\gamma}=\pm 1$, where $-1$ corresponds to a saddle point. Other critical points have index $+1$. If there is no critical point inside the contour $\gamma$ then the index is zero. The categorization of the critical point lets one determine the number of streamlines to emanate from the irregular node, and hence its valence in the mesh. When $i_{\gamma}=-1$ there are five separatices. When $i_{\gamma}=+1$ there are three [1, Lemma 5.1], [10]. It is possible to have an index of $+2$ at which six elements would meet. However such points are unstable and in practice will split into separate critical points of lower valence. For this reason, we do not consider such critical points here.

Using the definition (2) and relation (11), we similarly define the integral

[TABLE]

where $c$ is a (small) counter-clockwise circular contour centered on the critical point. Table 1 shows the values of $I_{c}$ for three separatrix topologies in the neighborhood of a critical point. It is equal to the sum of the jumps in $\psi$ around the contour. The valence $\cal V$ of the critical point is then

[TABLE]

Note, for consistency, that if there is no critical point inside the contour then $I_{c}=0$ and the point is regular with valence four.

If one generalizes the integral $I_{c}$ to be over a part of the circle, one can characterize boundary vertices, too [13]. Let

[TABLE]

so that $I_{c}=I\left(0,2\pi\right)$. At a boundary point, then, the number of quadrilaterals attached to the point (the boundary itself already being a streamline by construction) is

[TABLE]

where $\Delta\theta=\theta_{f}-\theta_{0}$. $I\left(\theta_{0},\theta_{f}\right)$ can be viewed as a correcting factor to $\Delta\theta$ towards a multiple of $\frac{\pi}{2}$.

2.4 Streamline Integration

The guiding field is also used to trace streamlines that will form the separatrices of the block decomposition. We can formulate the problem as finding the trajectory $\vec{x}(t)$ that satisfies

[TABLE]

where $t$ is the integration parameter and $\tilde{v}\left(\psi^{\prime}(\vec{x})\right)$ is the adjusted guiding field vector. Equation (2) gives a guiding field $\psi\in\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$, but crosses are invariant to rotations of $\pi/2$. It is therefore necessary to account for the branches of the cross that don’t lie within $\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$. These adjusted guiding directions can be obtained by adding a certain number $k=0,1,2,3$ of $\pi/2$ so that $\psi^{\prime}$ lies within the appropriate quadrant $\left[-\frac{\pi}{4}+k\frac{\pi}{2},\frac{\pi}{4}+k\frac{\pi}{2}\right]$. This adjustment operation is especially important when streamlines cross jump lines where $\psi$ abruptly changes by $\pm\pi/2$, but the streamline has to continue in the same overall direction.

The ordinary differential equation (16) requires initial conditions in the form of

[TABLE]

where $t_{0}$ is the initial integration time and $\vec{x}_{0}$ the location of the start of the streamline. In this case, the start of the streamline will either be an irregular node or a corner from which a non-boundary separatrix emerges. The ODE (16) can then be integrated in a traditional numerical manner as will be described below in Section 3.3.

3 The Mesh Generation Process

From the discussion above, we see that the decomposition of a domain into quadrilateral subdomains has four stages:

Computation of a guiding field. This includes the generation of a triangular mesh and solution of the Laplace problem (4). 2. 2.

Finding critical points of the guiding field and their valences, including the valences of corners at boundary vertices. 3. 3.

Integration of streamlines, from critical points and corners, to generate the interior boundaries of the quadrilateral decomposition. 4. 4.

Cutting of the domain into quadrilaterals using the separatrices as subdomain boundaries, then further subdivision into elements, as necessary.

This section presents details of the implementation of the process and walks through the four stages needed to generate a quadrilateral mesh. The procedure is illustrated using a geometry commonly used in the cross field literature [10, 11], that of a half disc as shown in Fig. 7.

3.1 Solution of the Field Equations

To solve the Laplace problem (4) we use the open source spectral element program Nektar++ [16, 17]. We first generate a triangular finite element mesh which is then made high order and curved by projecting interior nodes onto the curved boundaries [19]. This procedure is carried out in NekMesh [17], which also has the capability to optimize the high order mesh should some elements be of low quality or simply invalid [20]. Such a mesh is shown in Fig. 8 for the reference geometry.

For the purposes of generating the quadrilateral decomposition, we generate as coarse a triangular mesh as possible and then use high order polynomials to approximate the solution within the elements to get the desired accuracy. First, we want to take advantage of the spectral accuracy of the spectral element Laplace solver, where a fine mesh is not necessary. Second, as will be discussed in Sections 3.2 and 3.3, a small number of elements is desirable to simplify and speed up the critical point detection and streamline tracing. For the half circle geometry and triangle mesh shown in Fig. 8, it was sufficient to use third order polynomials (fourth order convergence with element size).

Once the curvilinear triangular mesh is generated, the Laplace problem (4) is solved. A special implementation of the Laplace solver has been implemented in the spectral element code Nektar++ [16] for the boundary value problem (4), where the boundary conditions are automatically computed using (3).

3.1.1 Choice of a Discretization

Nektar++ supports both continuous Galerkin and discontinuous Galerkin discretizations. The CG approximation is used when the boundary conditions (BCs) are continuous along all boundary curves. As can be seen in the discussion in Sec. 2, boundary conditions are continuous if the boundaries satisfy one of two conditions:

The curve at each point is $C^{1}$-continuous, i.e. it is a smooth curve; or

• 2.

The boundary curve is only $C^{0}$-continuous, i.e. a corner, and the angle is a multiple of $\pi/2$.

In this paper, we employ CG whenever possible and switch to DG only when the geometry requires discontinuous BCs. For example, the reference geometry of Fig. 7 contains only smooth curves and $\pi/2$ corners. We therefore use a CG formulation for it. The solution for $\vec{v}$ on the half disc geometry is shown in Fig. 6. As expected, boundary conditions and the solution are smooth, so the solution benefits from rapid convergence of the SEM.

If neither smoothness condition is satisfied at a boundary corner point, the boundary condition (3) has a jump discontinuity there. In that case, we use a DG formulation, which can account for discontinuous BCs. Unlike the CG approximation used in traditional low order solvers, no ad hoc smoothing of a corner’s BCs is required, leading to a discretization-consistent solution.

As an example of a geometry that does not satisfy the smoothness conditions, we present the polygon geometry shown in Fig. 9, which shows the solution and $\psi$ fields.

We note two consequences of using a DG formulation to handle corners with angles that are not multiples of $\pi/2$. First, jump lines can originate from corners. This fact is important, for it shows that one cannot a priori and ad hoc determine the valence of a boundary corner solely from its angle. Fig. 10 (left) shows how the jump line originating from the top corner
A in the geometry ends in the sole critical point of the domain. Although we have observed that curvature produces critical points, see Sec. 2.3, this straight-sided example shows that a corner whose angle is not a multiple of $\pi/2$ can also generate interior critical points, again, see Fig. 10 (left).

Second, observe that discontinuous BCs are naturally enforced in the DG formulation. No ad hoc smoothing of the boundary is required, as shown in Fig. 10 (right), because the Dirichlet BCs can be enforced either through fluxes when one triangle element used to compute the guiding field shares both edges of the corner (see lower corner,
B ), or exactly for all other corners (see upper corner,
C ).

When a DG method is used to solve the Laplace problems for the guiding fields, the fields are discontinuous between elements. This is especially true in domains where discontinuous BCs are present, as is clear in Fig. 10 (right). Enough resolution (either through h- or p-refinement) needs to be provided so that the jumps at the element interfaces are small. Too low resolution could result in zero isocontours located within solution jumps between elements. The current implementation would not be able to detect such isocontours, nor could it consequently detect critical points. See Sec. 3.2 for more details.

Fortunately, critical points should naturally be located in regions of smooth fields where jumps between elements are small. Discontinuities between elements are large only near boundaries where discontinuous BCs appear. The field on the inside of the domain is mostly smooth, due to the properties of the Laplace equation, resulting in very small jumps between elements. If strong curvature is present, we observe that critical points can appear near curved boundaries where BCs are continuous and the field is smooth. If no (or little curvature) is present, critical points are controlled by discontinuous BCs and are seen far away from the boundaries, on the inside of the domain, again where the field is smooth. In all cases, critical points are not expected to naturally occur in areas where large jumps between elements exist due to nearby discontinuous BCs.

Subsequently, the integration of streamlines is unaffected by the discontinuous nature of the discretization. By using a multi step high order integrator, as explained later in Sec. 3.3, streamline direction is not strongly affected by abrupt changes of direction between elements. Streamlines remain smooth even when gradients in the guiding field are large, e.g. when crossing a large jump between elements.

3.2 Critical Point Detection and Valence Evaluation

Once the solution to (4) is computed, zeros in the field $\vec{v}$ are found and their valences are evaluated. Unlike the usual cross field approach, we never, in fact, generate crosses. As seen in Sec. 2 the analysis can be performed on the $\vec{v}$ field, and $\psi$ is computed only when necessary, i.e. for computing the valence of vertices and tracing streamlines via (15) and (16). All other operations can be accomplished by operating on $\vec{v}$.

Critical points inside the domain are located at $\vec{v}=\vec{0}$. To find those, we first search all elements and flag each that contains a $\vec{v}=\vec{0}$ point at any quadrature point in the element. If we find at least one value of each $u\leq 0$, $u\geq 0$, $v\leq 0$ and $v\geq 0$, a critical point is either inside this or a nearby element. In each flagged element, we then use Newton’s method to approach the location of $\vec{v}=\vec{0}$. To do so efficiently, we perform the critical point search in parametric space. If a critical point is found outside the reference element, we dismiss it and assume that it will be found through parametric search in a neighbor element. In Fig. 6, we can clearly see the two symmetrically located critical points where the black and white contours cross, consistent with the irregular nodes found by others, e.g. [10], for this geometry.

After completing the search for all critical points, a valence must be computed for each. Referring to Sec. 2.3 and (13), we notice that the integral $I_{c}$ can be evaluated solely by looking for the presence of jump lines. In a counter-clockwise manner, a positive jump (i.e. from $-\frac{\pi}{4}$ to $\frac{\pi}{4}$) leads to a negative integral $I_{c}=-1$, and vice versa. The sign of the jump can itself be determined solely by values of $\vec{v}$, without the need to compute $\psi$ or the cross field.

A jump must satisfy two conditions: $u$ must be negative and $v$ must change sign. A positive increase of $v$ indicates a positive jump, and vice versa. To detect jumps, and therefore compute the valence of a critical point, we step counter-clockwise over a sequence of uniformly distributed quadrature points located on a circle of small radius $c$, centered on the critical point. Following this logic, we determine that the critical points seen in Fig. 6 both have a valence of three.

The valences of all boundary vertices are computed in a similar fashion, this time from expression (15). The angle at the corner, determined from $\theta_{0}$ and $\theta_{f}$, is computed from the geometry. To compute the integral $I_{c}$ we compute $\Delta\psi=\psi_{f}-\psi_{0}$ from the values of $\psi$ at the boundaries from the boundary conditions. To $\Delta\psi$, we add a jump contribution, if present, in the same way as for the computation of critical point valences. Because both corners in Fig. 6 have an angle of $\frac{\pi}{2}$, a valence of one can be trivially determined.

It is important to note that the boundary vertex operation does not require an a priori determination of the corner valence based on angle only. The corner valences can be determined through equation (15) which does not only depend on the known BCs but also on the existence of jump lines in the computed $\psi$ field in the neighborhood of the vertex, like that seen in Fig. 10.

The boundary vertex valence gives the number of quads located at the associated corner. Due to the definition, this valence could be zero. Indeed, $\Delta\theta$ and $\Delta\psi$ could cancel out if the corner is sufficiently sharp. Physically, this indicates the presence of a degenerate quadrilateral block (i.e. a triangular block) where all streamlines would converge towards the degenerate corner. This topology requires ad hoc manipulation where the streamlines are used later to construct the quadrilateral subdivision of the domain. This manipulation will be explained in the next section.

3.3 Streamline Integration and Manipulation

Streamlines are traced throughout the domain after the valences of all irregular nodes and vertices are determined. The first step is to find the initial direction of each streamline. Because $\vec{v}=\vec{0}$ at critical points, $\psi$ is undefined. Therefore, we choose to evaluate one of the streamline angles from $\vec{v}$ at a small distance $c$ of the critical point. This distance is currently set empirically based on the size of the elements.

The initial direction is refined iteratively from an initial guess using Algorithm 1, which is inspired by bisection. For the first streamline, we can take any initial guess for the direction. For subsequent streamlines, we take an initial guess at angles multiples of $\frac{2\pi}{\mathcal{V}}$. The size of the tolerance, $\epsilon$ is not critical; we use $10^{-9}$.

The angle search in Alg. 1 can also be applied to boundary vertices, including those where the boundary conditions might be discontinuous and $\psi$ is ambiguous. An initial guess for all streamlines may be obtained at angles multiples of $\frac{\Delta\theta}{\mathcal{V}}$.

After an initial direction is obtained, streamlines are synchronously advanced throughout the domain. This part of the procedure is the most computationally expensive. Each new streamline point requires the search for the element that contains it. Next, the inverse map $\vec{\xi}=\mathcal{X}^{-1}_{e}(\vec{x})$ is used to transform to parametric coordinates. Finally the high order interpolation of $\vec{v}$ is computed via eq. (8).

A multi stage integration like Runge-Kutta, where intermediate evaluations are used, would prove too expensive to integrate the streamlines. Instead, we use a $4^{th}$ order multi step Adams-Bashforth integrator. The current implementation reduces the order of the algorithm when too few steps are available. While advancing the streamlines, it is important to validate the direction $\alpha$ obtained at the latest point. When a streamline crosses a jump line, $\psi$ abruptly rotates $\pm\pi/2$ and $\alpha$ must be adjusted accordingly.

As streamlines are advanced throughout the domain, they may meet. To anticipate two streamlines meeting, the front points of each are compared at each step. If the distance between the front points of two streamlines is less than the step size and they are advancing in opposite directions, the streamlines are assumed to meet. When the absolute value of the difference in the latest $\alpha$ values equals $\pi$ when rounded to the nearest $\pi/2$, we consider the streamlines to be advancing in opposite directions. Each of two meeting streamlines are then advanced to the starting point of the other, keeping the same number of points in each.

A merged streamline is created by weight-averaging a pair of streamline points. A simple linear weight function could be used, but that will change the angles computed with Alg. 1 at the starting points, both irregular nodes and boundary vertices. This not desirable as the PDE problem has already optimized these angles to be as equally distributed as possible. Instead we use trigonometric weight functions of the form $W_{0}(x)=\cos^{2}(x)$ and $W_{1}(x)=\sin^{2}(x)$, where $x\in[0,1]$ is the parametrization along the streamline. Because $\frac{dW_{0}}{dx}$ and $\frac{dW_{1}}{dx}$ are zero at $x=0$ and $x=1$, angles are preserved at the initial and final points of the merged streamline.

3.3.1 Aggressive Merging

What has been described above constitutes expected or normal merging. But the user may also want to reduce the number of separatices and thereby simplify a valid block decomposition. One approach is to perform aggressive merging of the separatrices. For aggressive merging, a larger distance threshold is carefully chosen to make streamlines that would normally just miss each other merge anyway. Trigonometric weight functions are especially useful when performing aggressive merging of two reasonably distant streamlines. An example of aggressive merging with trigonometric weight functions is presented in Fig. 11 where the left figure corresponds to normal merging and the right one to aggressive merging. In this case, a distance threshold of 5 times the step size was used for aggressive merging. This example shows that the separatrix graph is greatly simplified and that the curving of some streamlines helps preserve the original direction at each end. Note, importantly, that irregular nodes obtained from the critical points are not modified a posteriori; only separatrices are.

3.4 Quad Mesh Generation

NekMesh relies on the OpenCASCADE platform [21] as its CAD engine for mesh generation and for the projection to high order boundary representations. OpenCASCADE also includes tools for CAD manipulation. We make use of these capabilities for spline and wire creation and for face splitting.

The separatrices computed from the guiding field are transformed into OpenCASCADE edges represented by interpolating splines. These edges can be joined into a set of wires, with each wire consisting of topologically connected edges. The original CAD file is loaded again, consisting of a single two dimensional face, representing the domain $\Omega$. This face is then iteratively split by each wire and a set of quadrilateral faces is obtained. These faces are topologically connected, meaning that the future mesh will be conforming. At this point, each face can be trivially meshed with a single quadrilateral and projected to high order. Fig. 12 (left) illustrates the quadrilateral decomposition and coarse quadrilateral mesh obtained for the half circle after computing the location of irregular nodes and tracing the separatrices.

We use the midpoint division approach [22] to split a triangular block that forms when a corner valence is zero. We insert an artificial 3-valence irregular node which we connect to each of the three sides of the triangle. One separatrix can then be physically integrated away from the degenerate corner and throughout the domain. The other two branches of the 3-valence irregular node are defined as straight lines at $\pm\frac{2\pi}{3}$ from the physical separatrix. Fig. 13 shows what the midpoint division of a degenerate quadrilateral might look like.

If desired, we make use of an isoparametric splitting approach described in [23] to further subdivide the quadrilateral subdomains into smaller elements. This technique, initially developed for boundary layer mesh division, guarantees that subdivided elements have the same quality as the macro element from which they are created. In the current context, it also ensures that angles are preserved. Each row of quads can be individually split using the isoparametric approach. The right-hand side of Fig. 12 shows what the mesh on the left looks like after a conformal splitting of each row of quads.

4 Examples

We present four examples of meshes generated for geometries taken from, or inspired by, the literature. The meshes shown in Fig. 12 are typical of those already generated by low order approaches [1, 10, 11]. We use the additional examples to illustrate the application of the high resolution approach on multiply connected geometry, one for which low order approaches have been shown to have difficulty, a polygon with both acute and obtuse angles, and a NACA 0012 airfoil.

The geometry for Example I shown in Fig. 14 is a rectangle with two quarter circle holes. This geometry has also been meshed, for example, in reference [10].

Some low order implementations have produced spurious asymmetric irregular nodes, see Fig. 14 [24], for example. The high order guiding field approach, however, produces the minimum number of irregular nodes (two 3-valence irregular nodes) while preserving the symmetry of the geometry without the need to coalesce spurious irregular nodes that can result from low order detection algorithms. Fig. 15 shows the coarse quadrilateral mesh obtained on this geometry using the guiding field approach.

Example II, is the nautilus. The domain and guiding fields shown in Fig. 16 give the decomposition shown in Fig. 17 when agressive merging is applied. This domain is one that low order methods generate spurious spiral separatices that end in a limit cycle, therefore creating an invalid block decomposition [1, 11]. To date, no approach using cross fields has been presented that generates a valid separatrix graph. The use of high order integration of an accurate guiding field allows NekMesh to generate separatrices that do not form a spiral. In fact, the irregular nodes found by NekMesh on this geometry are identical to those one would obtain for a simple disc. Fig. 17 shows the coarse quadrilateral mesh obtained on this nautilus geometry. This decomposition was generated using aggressive merging with a distance threshold of value 5 times the step size. If we were to use normal merging, the block decomposition would have looked more complicated, as shown in Fig. 11 (left), reducing the flexibility for users to split the mesh to their liking.

We note that the nautilus has a similar (mirrored) guiding field to that obtained on the half disc in Fig. 6. The nautilus itself consists of two half discs of different size attached by their chord. The irregular node pattern (four 3-valence nodes) is therefore expected. We also note one major difference with the irregular nodes obtained with previous techniques, e.g. references [1, 11]: The upper left irregular node is located further outside than in previous works. This seems to lead to streamlines successfully escaping a limit cycle.

Example III consists of a polygon whose corners include acute and obtuse angles. We use this example to demonstrate the use of the DG discretization for the solution of the guiding field. The rest of the meshing procedure remains the same and is in fact unaffected by the type of discretization. Fig. 9 has already shown the solution $\vec{v}$ and computed $\psi$ field for this geometry, with $u=0$ and $v=0$ isocontours shown in white and black, respectively. As mentioned earlier in Sec. 3.3, the only consideration when using the DG formulation is that the approximate solutions are accurate enough to not affect the streamline integration.

The Example III geometry features a sharp corner (far right, Fig. 18) whose valence is evaluated as zero. That creates a degenerate quadrilateral that must be split into three valid quads as explained in Section 3.3. Fig. 18 shows the coarse quadrilateral mesh obtained on the polygon geometry. It contains only two irregular nodes, one created ad hoc for the degenerate corner and the other detected as a critical point in the guiding field.

The coarse meshes can be split further according to the user’s preference using the isoparametric splitting described in Section 3.4. We show in Fig. 19 examples of two split meshes for the quadrilateral decompositions obtained above.

Finally, Example IV represents a geometry more relevant for CFD applications, that of a NACA 0012 profile in a rectangular domain. Since the trailing edge angle is not a multiple of $\pi/2$, a DG discretization was used to compute the guiding field. With this example, we demonstrate the flexibility of the isoparametric splitting module in NekMesh to define different distributions of elements. This allows the user to easily obtain elements of the required size, e.g. high aspect ratio elements in the boundary layer of the airfoil. Fig. 20 shows the coarse and the split quadrilateral meshes obtained on this NACA 0012 geometry.

5 Limitations and Extensions

The method that we have presented is inspired by cross field methods. Among the possibilities for computing what we are calling the guiding field, we solve a Laplace problem to smoothly propagate the boundary constraints to the interior. We differ from previous work in two ways: First, we use high resolution spectral element methods to solve for the guiding field; second, we do not compute the location of the critical points and separatrices from crosses, but instead from the continuous guiding field. As a result, the method has advantages over the low order techniques where numerical accuracy is the issue, but retains any disadvantages stemming from the original formulation of the problem. As a result there are both limitations and areas for further improvement.

The method has the same limitations due to the formulation as do cross field methods. One such limitation consists in the inability of cross fields to generate singular/critical points (and therefore irregular nodes) in domains where crosses are easily aligned with all nearby boundaries. Reference [14] documents a modified geometry of the nautilus (of Figs. 11, 16, 17) with a hole, for which a cross field is unable to generate a valid decomposition. Fig. 21 indeed shows that the field is aligned everywhere and that no critical point appears in the domain. Without critical points, however, a limit cycle is unavoidable and no valid quadrilateral decomposition can be obtained.

Also, the use of the high order approximations mitigates, but does not eliminate numerical errors. Therefore, although the nautilus geometry is easily decomposed with the high resolution approximation, it is still possible that spiral streamlines are generated if the accuracy is not high enough.

The limitations also point the way towards possible improvements to both the formulation and solution of the problem. Modifications to the guiding field formulation might be used to affect the number and location of critical points and even the existence of a valid quadrilateral decomposition. Modifications to the algorithm could be made to control errors.

The guiding field could be adjusted by adding a forcing term to (4), thus solving a Poisson rather than a Laplace problem, which would allow the guiding to be modified field in several ways. For example, critical points could be added [13, 12] so that a valid decomposition can be obtained for problems like the holed nautilus. (The degenerate triangle in Fig. 13 is also an example of placing an irregular point in the decomposition.) For the nautilus with hole, a hand sketch shows that adding one 3- and one 5-valent irregular nodes is sufficient. Existing critical points might also be moved away from boundaries to avoid thin blocks in the decomposition using the idea of control functions used in elliptic grid generation [25]. Finally, a forcing term would allow one to generate a guiding field that partially aligns with a metric field (incl. from a three dimensional surface), such as was demonstrated in [14].

Performance improvements can also be made. Streamline integration can be performed in reference space [26], akin to the location and analysis of the critical points. This would not only be computationally cheaper, it would also allow one to more easily handle geometries with large differences in scales. Where small scales are present, small elements are typically automatically generated for the background mesh to represent the geometry boundaries accurately. By performing all of the analysis of the guiding field in reference space, different scales would be automatically handled and better performance might be achieved. Streamline integration with a symplectic integrator might also be implemented to further mitigate the appearance of numerically created spirals. Finally, adaptivity and error control could be introduced to further mitigate numerical errors. Adaptive spectral element methods [27] could be used to automatically ensure that the guiding fields errors are acceptably small [28]. Adaptive ODE time integration methods could further enhance the accuracy of the streamline integration.

An extension of this methodology to volume block decompositions would represent a breakthrough in cross field based techniques. It would, however, require a reformulation of the problem based on frame (three dimensional crosses) arithmetic, which is currently lacking. Promising work in reference [29] points to a nine variable problem to solve.

6 Summary and Conclusions

We have developed a field guided method to generate quadrilateral meshes for general two dimensional domains. Inspired by cross field mesh generation methods, the procedure consists of four steps: The first is to compute the guiding field on an existing triangular mesh. For that, we use a spectral element method to compute a high resolution approximation to the two Laplace problems for the guiding field. We use either a continuous or discontinuous Galerkin method depending on the angles formed along the physical boundaries. If the angles are not multiples of $\pi/2$, the DG formulation can handle the discontinuous boundary conditions in a discretization consistent manner, without the need for ad hoc smoothing required by the CG formulation. Critical points in the field and their valences are then found by exploiting the local high order polynomial approximations within each element. Streamlines are computed in the field, again, using the high order solutions rather than crosses defined only at element corners. Finally, the domain can be cut into elements of the desired size, or used as is. The procedure is implemented in the open source program NekMesh, and uses the companion open source Nektar++ to compute the guiding fields.

An advantage of the approach is that it can generate meshes with naturally curved quadrilateral elements that do not need to be curved a posteriori to eliminate invalid elements. It avoids low order errors that make it difficult to locate irregular nodes and accurately trace separatrices.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675008. This work was supported by a grant from the Simons Foundation (#426393, David Kopriva). We thank Ms. Bing Yuan for the plots in Fig. 2. The first two authors would also like to thank Prof. Gustaaf Jacobs of the San Diego State University for his hospitality.