A spectral element method for meshes with skinny elements
Aaron Yeiser Advisor Alex Townsend





A spectral element method for meshes with skinny elements
    Aaron Yeiser
 MIT, Cambridge, MA 02139, US () 
 Advisor: Alex Townsend 
 Department of Mathematics, Cornell University, Ithaca, NY 14853.
    December 30, 2023
=============================================================================================================================================




  When numerically solving partial differential equations (PDEs), the first step is often to discretize the geometry using a mesh and to solve a corresponding discretization of the PDE.
  Standard finite and spectral element methods require that the underlying mesh has no skinny elements for numerical stability.
  Here, we develop a novel spectral element method that is numerically stable on meshes that contain skinny elements, while also allowing for high degree polynomials on each element.
  Our method is particularly useful for PDEs for which anisotropic mesh elements are beneficial and we demonstrate it with a Navier–Stokes simulation.
  Code for our method can be found at <https://github.com/ay2718/spectral-pde-solver>.



spectral element method, skinny elements, ultraspherical spectral method



15A03, 26C15




§ INTRODUCTION

A central idea for numerically solving partial differential equations (PDEs) on complex geometries is to represent the domain by a mesh, i.e., a collection of nonoverlapping triangles, quadrilaterals, or polygons. The mesh is then used by a spectral element method to represent the solution on each element by a high degree polynomial <cit.>. Engineering experience suggests that meshes containing skinny or skewed elements are undesirable and should be avoided, if possible <cit.>. However, from a mathematical point-of-view, long and thin elements can better resolve anisotropic solutions <cit.> that appear in computational fluid dynamics and problems involving boundary layers.

In this paper, we develop a spectral element method that is numerically stable on meshes containing elements with large aspect ratios and also allows for high degree polynomials on each element. In particular, the numerical stability of our method is observed to be independent of the properties of the mesh and polynomial degrees used. This has several computational advantages:  (1) A practical hp-refinement strategy does not need to take into consideration the quality of the final mesh, (2) Anisotropic solutions can be represented with fewer degrees of freedom, and (3) Low-quality meshes can be used without concern for numerical stability.

Over the last two decades, there has been considerable research on mesh generation algorithms that avoid unnecessarily skinny elements <cit.>. While these mesh generation algorithms are robust and effective, they can be the dominating computational cost of numerically solving a PDE. By allowing low-quality meshes, our method alleviates this current computational burden on mesh generation algorithms. This is especially beneficial for time-dependent PDE in which an adaptive mesh is desired, for example, to resolve the propagation of a moving wave front.

There are a handful of mesh generation algorithms that have been recently designed for anisotropic solutions that generate meshes with skinny elements <cit.>. It is therefore timely to have a general-purpose element method that allows for an accurate solution on such meshes. For time-dependent
PDEs, such as the Navier–Stokes equations, skinny elements can relax the Courant–Friedrichs–Lewy time-step restrictions, leading to significant computational savings in long-time incompressible fluid simulations.

A standard element method constructs a global discretization matrix with a potentially large condition number when given a mesh containing skinny elements. This is because the global discretization matrix is assembled element-by-element from a local discretization of the PDE on each element. Since the elements may have disparate areas, the
resulting discretization matrix can be poorly scaled. One can improve the condition number of the global discretization matrix by scaling the local discretization matrices by the element's area.  Since we are using a spectral element method, we can go a little further and scale the local discretization matrices by the Jacobian of a bilinear transformation that maps the element to a reference domain (see Section <ref>).

The solutions of uniformly elliptic PDEs with Dirichlet boundary conditions defined on a convex skinny element are not sensitive to perturbations in the Dirichlet data or forcing term.
Since we know it is theoretically possible to solve elliptic PDEs on meshes with skinny elements, it is particularly disappointing to employ a numerically unstable method (see Section <ref>).

The structure of the paper is as follows:
In Section <ref>, we define some useful terminology and provide the theory to show that certain elliptic PDEs on skinny elements are well-conditioned.
Section <ref> describes the process of using our spectral element method on a single square element, and Section <ref> extends that process to convex quadrilateral domains.
In Section <ref>, we define a method for joining together two quadrilateral domains along an edge, and Section <ref> extends that method to joining an entire quadrilateral mesh.
Key optimizations to our method are described in Section <ref>.
Finally we present a Navier–Stokes solver using our spectral element method in Section <ref>.




§ UNIFORMLY ELLIPTIC PDES ON SKINNY ELEMENTS

It is helpful to review some background material on skinniness of elements and the condition number of PDEs on skinny elements.
Given an element 𝒯 of any convex 2D shape, let r_in denote the inradius, i.e., the radius of largest circle that can be inscribed in 𝒯, and
r_out denote the min-containment radius, i.e.,  the radius of smallest circle that contains 𝒯. We define the skinniness of 𝒯
to be

    s( 𝒯) =  r_in/r_out.

For finite element methods, it is widely considered that a skinny triangle with a large angle, i.e., θ≈π, is more detrimental to numerical stability than those triangles with only acute angles. This is based on the theoretical results on interpolation and discretization errors in <cit.>.  The spectral element method that we develop is numerically stable on meshes containing any type of skinny quadrilateral.  Therefore, it is also numerically stable on meshes containing any kind of skinny triangle, because a triangle can be constructed out of three non-degenerate quadrilaterals by dissecting along the medians, as in Figure <ref>.



 §.§ Apriori bounds








Consider the second-order elliptic differential operator operating on domain Ω.

    ℒ u = ∑_i,j=1^d a_ij(x) ∂^2 u/∂ x_i∂ x_j + ∑_i=1^d b_i(x)∂ u/∂ x_i+ c(x) u,

where the coefficients a_ij, b_i, and c are continuous functions on Ω. We say that ℒ is uniformly elliptic if for some constant θ_0 > 0,

    ∑_i,j=1^d a_ij(x) ξ_iξ_j ≥θ_0 ∑_i=1^dξ_i^2,       x∈Ω,     ξ∈^d.




Consider the uniformly elliptic partial differential operator

    ℒ u = ∑_i, j = 1^n a_ij u_x_i x_j + ∑_i = 1^n b_i u_x_i,

with coefficients a_ij and b_i continuous.  Let ℒ u = f on domain Ω with boundary ∂Ω.
If f ≥ 0 on Ω, then the maximum value of u occurs on ∂Ω.
If f ≤ 0 on Ω, then the minimum value of u occurs on ∂Ω.



For a proof, see Theorem 1 in Section 6.5.1 of <cit.>.




Consider the uniformly elliptic partial differential operator

    ℒ u = ∑_i, j = 1^n a_ij u_x_i x_j + ∑_i = 1^n b_i u_x_i + cu,

with coefficients a_ij, b_i, and c continuous.  Let ℒ u = f on domain Ω with boundary ∂Ω, and let c < 0 on Ω.
If f ≥ 0 on Ω, then the maximum positive value of u occurs on ∂Ω.
If f ≤ 0 on Ω, then the minimum negative value of u occurs on ∂Ω.



For a proof, see Theorem 2 in Section 6.5.1 of <cit.>.


Theorem <ref> uses Lemmas <ref> and <ref> to show that perturbations to the boundary conditions of a uniformly elliptic PDE with c ≤ 0 will not be amplified in solution on the interior of the domain.



Let ℒ u = ∑_i, j = 1^n a_ij u_x_i x_j + ∑_i = 1^n b_i u_x_i + cu, with c ≤ 0.  Let u be the solution to ℒ u = f on domain Ω with boundary ∂Ω, and u(∂Ω) = g.  If v is the solution to ℒ u = f with boundary conditions u(∂Ω) = g + ϵ, then max |v - u| ≤max |ϵ|.



Let us write v = u + δ.  Since L is a linear operator, we know that ℒ v = ℒ u + ℒ u δ = f.
Therefore, ℒδ = 0, with δ(∂Ω) = ϵ.  By Lemma <ref>, we know that if δ achieves a nonnegative maximum on the interior of Ω, then δ is constant.  Similarly, if δ achieves a nonpositive minimum on the interior of Ω, then δ is constant.  If max |δ| > max |ϵ|, then δ will achieve a nonpositive minimum or a nonnegative maximum on the interior of Ω.  However, Lemma <ref> states that δ must then be constant, so therefore max |δ| = max |ϵ|, hence a contradiction.
Therefore, max |δ| ≤max |ϵ|.


Theorem <ref> provides a weak bound for the effect of perturbations to the right hand side on uniformly elliptic PDEs with c ≤ 0.



Let Ω∈ℝ^d be a bounded domain and ℒ a second-order uniformly elliptic differential operator with c≤ 0. Let u be a weak solution of

    {ℒ u    = f   on  Ω,
    
    u    = g   on  ∂Ω,
    .

and û be a weak solution of

    {ℒû   = f + Δ f   on  Ω,
    û   = g + Δ g   on  ∂Ω,
    .

where u,û∈𝒞^2(Ω)∩𝒞^0(∂Ω) and ∂Ω is the closure of Ω. Then,

    sup_x∈Ω |u(x)-û(x)| ≤sup_x∈∂Ω |Δ g(x)| + (e^(β+1) diag(Ω)-1)sup_x∈Ω|Δ f(x)/θ(x)|,

where β = max_i=1^d sup_x∈Ω |b_i(x)/θ(x)|.


Since u satisfies ℒu = f on Ω and u = g on ∂Ω, Δ u = û - u satisfies ℒ(Δ u) = Δ f on Ω and Δ u = Δ g on ∂Ω.
One can apply the bound in Theorem 3.7 of <cit.>. 


Theorem <ref> shows that solutions to uniformly elliptic PDEs with c≤ 0 are not sensitive to perturbations in the Dirichlet boundary data and Theorem <ref> shows the sensitivity of the solution to perturbations in the right–hand side does not depend on the aspect ratio of the domain Ω.
However, this bound is extremely weak.  We provide a better bound for the screened Poisson equation with Theorem <ref>.




Let ℒ u = u_xx + u_yy - k^2 u, and let u be the solution to ℒ u = f on domain Ω with boundary ∂Ω, and u|_∂Ω = g.  Let r be the radius of the smallest circle completely containing Ω.  If s is the solution to ℒ s = f + ϵ with boundary conditions s|_∂Ω = g, then max |s - u| ≤max |ϵ| r^2/4.


First, consider the smallest disk ω with radius r that completely contains region Ω.  Let this disk have center (a, b).
For (x, y) ∈ω and c < 0, the solution t = c(x - a)^2 + (y - b)^2 - r^2 is nonnegative, max t = -r^2 c, and t_xx + t_yy = 4c.

Next, consider solution w such that w_xx + w_yy = 4c for c < 0 and w = 0 on ∂Ω.
The solution w - t satisfies (w - t)_xx + (w - t)_yy = 0 and w - t ≤ 0 on ∂Ω, since t ≥ 0 for all (x, y) ∈Ω.
By Lemma <ref>, (w - t) ≤ 0 for all (x, y) ∈Ω.
Therefore, max w ≤ -r^2 c.

Finally, we will consider the solution v such that v_xx + v_yy = ϵ with v = 0 on ∂Ω and max |ϵ| = c.  By Lemma <ref>, v is bounded above by w_+ satisfying ∇^2 w_+ = -c and w_+ = 0 on ∂Ω, and v is bounded below by w_- satisfying ∇^2 w_- = c and w_- = 0 on ∂Ω.
We can write ϵ = ϵ_+ + ϵ_-, where ϵ_+ = max(ϵ, 0), and ϵ_- = min(ϵ, 0).

Let δ_+ be the solution to ℒδ_+ = ϵ_- with zero Dirichlet boundary conditions and let δ_- be the solution to ℒδ_- = ϵ_+ with the same boundary conditions.  By Lemma <ref>, δ_+ ≥ 0 on Ω, and δ_- ≤ 0.

We can write that ∇^2 δ_+ = k^2 δ_+ + ϵ_-.  Since δ_+ ≥ 0, min(k^2 δ_+ + ϵ_-) ≥min(ϵ_-), so δ_+ is bounded above by the solution of ∇^2 v = ϵ_- with zero Dirichlet boundary conditions by Lemma <ref>.  Similarly, δ_- is bounded below by the solution of ∇^2 v = ϵ_+ with zero Dirichlet boundary conditions.  Therefore, s = u + δ_+ + δ_-.
Thus, max |s - u| = max(maxδ_+,  max -δ_-) ≤max|ϵ|r^2/4.


Theorem <ref> proves that any change to the boundary conditions of the screened Poisson equation will not be amplified in the solution on any domain.
Theorem <ref> gives a weak bound for solving all second–order uniformly elliptic PDEs with c ≤ 0, and Theorem <ref> shows that the amplitude of a perturbation in the right hand side of the screened Poisson equation is bounded by a function of the radius of the domain.  Therefore, solving the screened Poisson equation on any domain is a well-posed problem.

It seems unfortunate to have an unstable numerical algorithm for a well-conditioned problem. Therefore, we ask ourselves the question: Can we create a numerically stable algorithm for solving PDEs on meshes with skinny elements?


 §.§ Stability and conditioning


Our method represents any linear differential operator as a matrix operator.  One method of assessing the numerical stability of the method is to examine the condition number of the matrix.
In order to demonstrate that the condition number is bounded, let us define a skinny triangle with coordinates (0, 0), (1, 1 + ϵ), (2, 2 - ϵ).  We can now define a skinny quadrilateral with vertices (0, 0), (0.5, 0.5 + 0.5ϵ), (1, 1), (1, 1 - 0.5ϵ).
When we construct a differential operator matrix on that skinny quadrilateral, we still need a preconditioner.  In fact, we found that row scaling so that each row has a supremum norm of 1 is fairly close to optimal.
Figure <ref> shows the condition number κ of the normalized matrix for solving Poisson's equation as ϵ approaches zero.
Even as the quadrilateral becomes very skinny, the condition number of the matrix is bounded from above.
The process of creating the matrix is shown in sections <ref> and <ref>.



Figures <ref> and <ref> show that our method gives extremely accurate results even on meshes with skinny triangles.  Figure <ref> shows that our method performs equally well on meshes with or without skinny elements.










§ SPECTRAL METHODS

In order to develop a useful spectral method, one must choose the proper set of basis functions.
Although wavelets and complex exponentials are tempting possibilities, certain sets of polynomials
promise simplicity and easy computability. Although polynomial interpolation is
numerically unstable using equispaced points, it is numerically stable using Chebyshev
points, with the nth set of Chebyshev points x_j = cosn-1-j/(n-1)π
 <cit.>.

Chebyshev polynomials are an efficient and numerically stable basis for
polynomial approximation <cit.>.  We define the nth degree Chebyshev
polynomial as T_n(x) = cos(narccos(x)), for x ∈ [-1, 1] <cit.>.



 §.§ The ultraspherical spectral method

Spectral methods are especially useful for solving linear differential equations.
The goal of using spectral methods is to approximate an infinitely-dimensional
linear differential operator with a finite-dimensional sparse matrix operator.
We can represent a function as an infinite vector of Chebyshev coefficients,
and we can truncate that vector at an appropriate size to achieve a sufficiently
accurate approximation of a smooth function.
If we attempted to create a differential operator D_λ that mapped a
function in the Chebyshev basis to its λ^th derivative in the Chebyshev basis,
then that operator would be dense, as most of its matrix elements would be nonzero.
Current spectral methods map Chebyshev polynomials to Chebyshev polynomials,
giving dense, ill-conditioned differentiation matrices.
In contrast, our method constructs D_λ to map a function in Chebyshev coefficients
to its λth derivative in the ultraspherical basis of parameter λ.
For every positive real λ, ultraspherical polynomials of parameter λ
are orthogonal on [-1, 1] with respect to the weight function (1 - x^2)^λ - 1/2
<cit.>.  Here, we will ultraspherical polynomials with positive integer parameter.

Using ultraspherical coefficients of parameter λ, the matrix λth
differentiation matrix D_λ is

    D_λ = 2^λ-1(λ-1)!
    [ 0   ⋯   0^λ times                 λ;                                                 λ + 1;                                                                   λ + 2;                                                                                         ⋱ ].

These operators are diagonal, but combining different-order
derivative coefficients leads to a problem:  different-order D_λ transform functions
in the Chebyshev basis to ultraspherical bases of different parameters.
The solution is to create basis conversion operators, denoted as S_λ,
which convert between ultraspherical bases of different orders.
Specifically, S_0 converts functions from the Chebyshev basis to the
ultraspherical basis of parameter 1, and S_λ converts functions from the
ultraspherical basis of parameter λ to the ultraspherical basis of parameter λ + 1.
For example, we can represent the differential operator ux + 2x as
the matrix operator (D_1 + 2S_0)u.  The ultraspherical conversion operators take the form

    S_0 = 1/2[  2  0 -1;     1  0 -1;        1  0  ⋱;           1  ⋱;              ⋱ ]
    ,   
    S_λ =
    [        1        0 -λ/λ + 2;           λ/λ + 1        0 -λ/λ + 3;                    λ/λ + 2        0        ⋱;                             λ/λ + 3        ⋱;                                            ⋱ ]
    ,   λ > 0.

Since D_λ and S_λ are sparse for all valid λ,
the final differential operator will be sparse and banded.
Let L be some n-order differential operator such that
L = c_0 + c_1D_1 + c_2D_2 + … + c_nD_n where c_0, …, c_n are constants
and D_n is equivalent to d^n/dx^n.
If we wish to solve the differential equation Lu = f, we create the
matrix operator 𝐋.  We also create vectors u⃗ and f⃗
as vectors of Chebyshev coefficients approximating u(x) and f(x).
Now, we can write
𝐋 = c_0S_n-1⋯ S_0 + c_1S_n-1⋯ S_1D_1 + ⋯ + c_n-1D_n-1 + c_nD_n,
and we can write the matrix equation 𝐋u⃗ = S_n-1⋯ S_0 f⃗.

This method is not limited to equations with constant coefficients.
It is possible to define multiplication matrices of the form M_λ[f], where M_λ[f] will multiply a vector in the ultraspherical basis of parameter λ by function f <cit.>.
The bandwidth of multiplication matrices increases with the degree of f, so high speed computation requires f with low polynomial degree.

In order to find a particular solution to the differential equation Lu = f, we need boundary conditions.
An nth order differential equation requires n distinct boundary conditions.
Since the highest ordered coefficients of f are often very close to zero,
we can simply replace the last n rows of 𝐋 and the last n elements of
f⃗ with boundary condition rows.
To set u(x) = a, we create the boundary condition row of [T_0(x), T_1(x), T_2(x), …, T_n-1(x)],
and the corresponding value of f⃗ is set to a.
If we use the boundary condition u'(x) = a, we create the boundary condition row
of [T'_0(x), T'_1(x), T'_2(x), …, T'_n-1(x)], and the corresponding value of f⃗ is set to a.
We now have a sparse, well-conditioned differential operator,
assuming that the problem is well-posed <cit.>.















































































 §.§ Discrete differential operators in two dimensions

In two dimensions, we can write any polynomial p as its bivariate Chebyshev expansion

    p(x, y) = ∑_i, j = 0^n a_ij T_i(y)T_j(x),      x, y ∈ [-1, 1].

To form a matrix operator, we can stack the columns of the coefficient matrix to form the vector

    u = [a_00 ,  a_10 ,  a_20 , … , a_n-1,0 ,  a_01 ,  a_11 , … , a_n-1,1 , … , a_n-1,n-1]^T.

We use Kronecker products (denoted as `⊗') in order to construct differential operators in two dimensions.
For example ux is represented as D_1 ⊗ I and uy is represented as I ⊗ D_1.
To put together differential operators, we can use the identity
(A ⊗ B) (C ⊗ D) = AC ⊗ BD,
so to represent an operator like ^3x ∂ y^2, we use D_1 ⊗ D_2.
We can also add Kronecker products to construct more complicated operators; for instance,

    ^2x^2 - ^2x ∂ y + 2x  ⟺  
    D_2 ⊗ S_1S_0 - S_1D_1 ⊗ S_1D_1 + 2S_1D_1 ⊗ S_1S_0.

A derivation of this two-dimensional sparse operator construction can be found in <cit.>.

Boundary conditions in two dimensions are slightly more complicated than in one dimension.
Consider solving partial differential equations on a square with resolution n × n.
We now remove rows corresponding to a right-hand side term of degree n-2 or n-1.
These 4n-4 rows are replaced with boundary conditions corresponding to the 4n-4 boundary points.

If we want to set Dirichlet conditions at the point (a, b), we create the boundary row

    [T_0(a), T_1(a), …, T_n-1(a)] ⊗ [T_0(b), T_1(b), …, T_n-1(b)].


In order to apply Neumann conditions, S_0^-1D_1u will give the first derivative of Chebyshev vector u in Chebyshev coefficients.  To find the value the derivative of u at point a, we use [T_0(a), T_1(a), …, T_n-1(a)]S_0^-1D_1u.
Thus, a boundary condition row representing ∂ u/∂ x at point (a, b) would be

    [T_0(a), T_1(a), …, T_n-1(a)]S_0^-1D_1 ⊗ [T_0(b), T_1(b), …, T_n-1(b)]

and a boundary condition row representing ∂ u/∂ y at point (a, b) would be

    [T_0(a), T_1(a), …, T_n-1(a)] ⊗ [T_0(b), T_1(b), …, T_n-1(b)]S_0^-1D_1.




§ SPECTRAL METHODS ON CONVEX QUADRILATERALS

The method of constructing matrices described above works only on square domains. We now adapt existing spectral methods to solve differential equations on quadrilaterals, and eventually, quadrilateral meshes.

We use the following bilinear map from the square [-1, 1]^2 to the quadrilateral with vertices (x_1, y_1)… (x_4, y_4) in counterclockwise order.

    x = a_1 + b_1r + c_1s + d_1rs,     
    y = a_2 + b_2r + c_2s + d_2rs,

where (x, y) refers to coordinates on the quadrilateral and (r, s) refers to coordinates on the square.

    a_1    = 1/4x_1 + x_2 + x_3 + x_4,       b_1 = 1/4x_1 - x_2 - x_3 + x_4,
    
    c_1    = 1/4x_1 + x_2 - x_3 - x_4,          d_1 = 1/4x_1 - x_2 + x_3 - x_4,

and a_2, b_2, c_2, and d_2 are similarly defined with y_1, y_2, y_3, and y_4.
The transformation from the quadrilateral to the square is more complicated, but it will not be needed.

When the quadrilateral domain is mapped onto the square, the differential equations on the quadrilateral are also distorted.  Through the use of the chain rule, we obtain

    u_x    = u_r r_x + u_s s_x,
    
    u_y    = u_r r_y + u_s s_y,
    
    u_xx   = u_rr(r_x^2) + 2u_rs(r_x s_x) + u_ss(s_x^2) + u_r r_xx + u_s s_xx,
    
    u_xy   = u_rr(r_x r_y) + u_rs(r_x s_y + s_x r_y) + u_ss(s_x s_y) + u_rr_xy + u_s s_xy,
    
    u_yy   = u_rr(r_y^2) + 2u_rs(r_y s_y) + u_ss(s_y^2) + u_r r_yy + u_s s_yy.

We still need to evaluate derivatives r_x, r_y, s_x, s_y, r_xx, r_xy, r_yy, s_xx, s_xy, and s_yy.
Using the Inverse Function Theorem, we can invert the Jacobian of the transformation to obtain

    [ r_x r_y; s_x s_y ]
    =
    [ x_r x_s; y_r y_s ]^-1
    =
    1/(x, y)[  y_s -x_s; -y_r  x_r ],
        (x, y) = x_r y_s - x_s y_r.

Using these equations, we can also find second derivatives r_xx, r_xy, r_yy, s_xx, s_xy, and s_yy.  These derivatives can be represented as a bivariate cubic polynomial divided by (r, s)^3.

In order to keep our differential operators sparse, we simply multiply the equations and right-hand side by (r, s)^3.
That way, we only multiply the operator by multiplication matrices with polynomial degree < 3, rather than by a dense approximation of a rational function, and our differential operators remain sparse and banded.



§ INTERFACE CONDITIONS

In order to actually do something useful with quadrilateral domains, it is necessary to somehow stitch them together into a larger mesh.
Here, we will first consider the case of stitching two quadrilaterals together along an edge.
We found that the most efficient way to stitch together multiple domains was to use the Schur Complement Method <cit.>.
Essentially, we first solve for the values along the interface between the two domains, and then we use that interface to find the bulk solution inside of each domain.  The advantage to using the Schur Complement Method is that once we have solved for the interface solution, we can compute the solution on each element in parallel.

In order to use the Schur Complement method on two quadrilaterals, we must have a linear system of the form

    [  A_11     0  A_1Γ;     0  A_22  A_2Γ; A_Γ 1 A_Γ 2  A_ΓΓ ][ u_1; u_2; u_Γ ]
    =
    [ f_1; f_2; f_Γ ].


Here, u_1 is the solution on the first quadrilateral, u_2 is the solution on the second quadrilateral, and u_Γ is the solution on the interface between u_1 and u_2.
We will represent u_Γ as a vector of values at Chebyshev points.

We also have A_11 and A_22 for solving partial differential equations on quadrilateral elements.
Finally, A_1Γ and A_2Γ constrain u_Γ to u_1 and u_2, and A_Γ 1 and A_Γ 2 ensure that the solution across the interface is differentiable.  A_ΓΓ is zero, as it is not necessary here.

First, we need to set A_1Γ and A_2Γ.
As shown in Figure <ref>, if an edge of u_1 is connected to u_Γ, then we connect all of the nodes on that edge of u_1 to u_Γ except for the node on the counterclockwise corner of that edge.
This arrangement ensures that all nodes on the exterior boundary of the domain are assigned the proper boundary conditions.
It also avoids redundant boundary conditions when multiple quadrilaterals meet at a point on the interior of the mesh.


Now, we need to ensure that the solution is differentiable across the interface with A_Γ 1 and A_Γ 2.
In order to do so, we need to find the derivative perpendicular to the boundary.
Let us say that the boundary has endpoints (x_1, y_1) and (x_2, y_2).
Let us define

    α = Δ x/√(Δ x^2 + Δ y^2),     β = Δ y/√(Δ x^2 + Δ y^2).


Let z = β x - α y, so the derivative normal to the boundary is u_z = β u_x - α u_y. We want this derivative in terms of r and s. Using (<ref>), we have

    u_x = u_r r_x + u_s s_x,      u_y = u_r r_y + u_s s_y.

From (<ref>), we have

    u_x = 1/(r, s)u_r y_s - u_s y_r,      u_y = 1/(r, s)-u_r x_s + u_s x_r.

By differentiating (<ref>), we have

    u_x    = 1/(r, s)u_r (c_2 + d_2 r) - u_s (b_2 + d_2 s)
    
        u_y    = 1/(r, s)-u_r (c_1 + d_1 r) + u_s (b_1 + d_1 s).

Finally, we can generate boundary rows for u_r and u_s with (<ref>) and (<ref>).
Each row of A_Γ 1 represents the derivative of u_1 perpendicular to the interface at one of the n Chebyshev points on the interface,
and each row of A_Γ 2 represents the negative derivative of u_2 at each Chebyshev point on the boundary.

Now, we want to ensure continuity across the boundary at the two endpoints. We can simply replace the differentiability conditions corresponding to the endpoints of the interface with continuity conditions from (<ref>).



§ CONSTRUCTING DIFFERENTIAL OPERATORS ON MESHES

The key to designing a mesh solver is strict bookkeeping.
Every quadrilateral in the mesh is given a unique number, as is every vertex and every edge.
We can define our mesh with a list of coordinates of vertices, and a list of the vertices contained in each triangle.
It is important that all of the vertices are numbered in counterclockwise order.
Next, each edge must be given a unique global number.
Edges also have local definitions—an edge can be expressed as two vertices, or as
a quadrilateral number and a local edge number on the quadrilateral.
We generate arrays to convert between these global and local definitions.
Finally, we must construct an array that keeps track of whether vertices are on the
inside or boundary of the domain.  Each vertex is also assigned a single edge
adjacent to that vertex.  This “vertex list” is essential for constructing non-singular
boundary conditions.

Next, we define u_Γ such that elements n(k-1)+1 : nk of u_Γ correspond to the kth interior edge. Thus, columns n(k-1)+1 : nk of A_iΓ and rows n(k-1)+1 : nk of AΓ i correspond to interface conditions across the kth interior edge.

We now generate a matrix A_i for each quadrilateral.
If quadrilateral i and j share an edge with interior edge number k, then columns n(k-1)+1 : nk of A_iΓ
and columns n(k-1)+1 : nk of A_jΓ are set like in the previous section.
Next, rows n(k-1)+1 : nk of AΓ i and rows n(k-1)+1 : nk of AΓ j are set to ensure differentiability, with the endpoints set to continuity conditions.
Finally, we check the two endpoints of the edge on the vertex list.
If an endpoint is listed as an interior vertex and the vertex list marks edge k, then the continuity conditions are removed for that endpoint and replaced with differentiability conditions across the edge at that endpoint.
That way, we make sure that the solution is fully continuous and almost fully differentiable without singularities in the linear system.



§ OPTIMIZATION

When solving systems of equations, banded matrices are typically efficient to solve.
Our matrices for finding solutions on quadrilaterals are “almost banded,” meaning that most of the elements lie close to the main diagonal, but a few rows extend the whole length of the matrix.
These dense rows correspond to boundary conditions, as defined in Section <ref>.
Therefore, sparse LU decomposition and Gaussian elimination have large backfill—they are forced to set many zero matrix elements to nonzero values, resulting in a much more computationally intensive solve.
We found that we can use the Woodbury matrix identity to efficiently replace the rows lying outside of the bandwidth,
creating a banded matrix that can easily be stored as an LU decomposition.
The Woodbury matrix identity can compute the inverse a matrix given the inverse of another matrix and a rank-k correction.  It satisfies

    (A + UCV)^-1 = A^-1 - A^-1U(C^-1 + VA^-1U)^-1VA^-1,

where A is n-by-n, U is n-by-k, V is k-by-n, and C is k-by-k <cit.>.

In our case, we regard our matrix as A + UCV, where A is a matrix that is fast to solve and UCV is a low–rank correction.
We generate V by stacking every row of A + UCV that corresponds to a boundary condition, since these rows are the only dense rows.
We define U as a sparse binary matrix with exactly one `1' in each column and no more than one `1' in each row, and C is simply the identity matrix.  The replacement boundary condition rows are 1's in columns corresponding to the lowest-order coefficients of the solution.  Our matrix A is now sparse and banded, so fast sparse LU decomposition is possible.  Therefore, once the LU decomposition is completed, back-substitution can be evaluated for each individual right-hand side extremely rapidly, allowing a system of variables with over 60,000 degrees of freedom to be solved in under a second.

Another place for optimization is the matrix 
    Σ = A_ΓΓ
     - ∑_j A_Γ j A_jj^-1 A_j Γ.

The matrix Σ is typically quite large, so it is important to be efficiently solvable.
Consider A_Γ j A_jj^-1 A_j Γ.
This matrix has nonzero columns only in columns of A_j Γ
that contain nonzero values, and it has nonzero rows only in rows of A_Γ j
that contain nonzero values.
Fortunately, most rows of A_Γ j and columns of A_j Γ are zero.
In order to optimize the bandwidth of Σ, we must optimize the placement of the nonzero elements of A_Γ j and A_j Γ.

In general, if there are m different boundaries between quadrilaterals, and n is the discretization size,
then u_Γ is composed of m blocks of n elements each, with each block
corresponding to a boundary between two quadrilaterals.  Each boundary is assigned a number, such that
elements n(a-1)+1  :  na of u_Γ correspond to boundary a.
If boundary a is contained by quadrilateral j, then A_Γ j has nonzero rows
n(a-1)+1  :  na and A_j Γ has nonzero columns n(a-1)+1  :  na.
Thus, if two quadrilaterals contain boundaries a and b, then Σ has a bandwidth of
at least (|a-b| + 1)n.  Therefore, in order to minimize the bandwidth of Σ, we
minimize the maximum |a - b| for any two boundaries contained by the same quadrilateral.
The theoretical minimum is (√(N)), where N is the number of elements, and
standard techniques can efficiently get pretty close to this minimum.  Therefore,
computing the LU decomposition of Σ has a time complexity of (N^2n^3)
and a space requirement of (N^1.5n^2).  The time complexity of back-substitution
is (N^1.5n^2).

Using the Woodbury identity for all of the matrices A_jj and computing A_jj^-1 A_j Γ
has a time complexity of (N n^4).  Back substitution for all A_jj has a time complexity of (N n^3), and the storage space for the matrices is also (N n^3).
In order to solve partial differential equations on a mesh, there is an (N^2 n^3 + N n^4) precomputation time, and for each right hand side, the solve time is (N^1.5 n^2 + N n^3).
It may be possible to use an iterative method to multiply a vector by Σ^-1, reducing the factors of (N^2)
and (N^1.5) to (N).



§ THE NAVIER–STOKES EQUATIONS



We decided to demonstrate our method with the direct numerical simulation of the Navier–Stokes equations in two dimensions at moderately high Reynolds numbers.
For incompressible flows, the Navier–Stokes equation are

    𝐮t + (𝐮·∇) 𝐮 + ∇ p = ∇^2 u,     ∇·𝐮 = 0,

where 𝐮 is the velocity vector field and p is the internal pressure field.
Note that these equations set density and viscosity to 1, so the Reynolds number can be varied solely by changing the velocity of the fluid or the dimensions of the domain.
Since these equations are coupled and nonlinear, we use the multistep method described in <cit.>.  We can use a first-order projection method to linearize and decouple the Navier–Stokes equations.  With no-slip boundary conditions, we have

    ∇^2 𝐮^n+1/2 - 𝐮^n+1/2/ t = (𝐮^n ·∇) 𝐮^n - 𝐮^n/ t,           𝐮^n+1/2 = 0,   on ∂Ω
       ∇^2 p^n+1 = 𝐮^n+1/2/ t,           p^n+1𝐧 = 0,   on ∂Ω
       𝐮^n+1 = 𝐮^n+1/2 -  t ∇ p^n+1,

where 𝐧 is the normal vector to the boundary <cit.>.

Navier–Stokes simulations are extremely useful in the form of wind tunnel simulations, where a test object is placed in a steady flow and analyzed.
For an incompressible flow, different boundary conditions are applied at the inlet, outlet, walls, and test object.
Typically, the test object will have no-slip boundary conditions, and the walls will have free-slip boundary conditions.  The inlet has a constant fluid velocity and zero pressure gradient, and the outlet has a constant pressure and zero velocity gradient.

Figure <ref> shows a wind tunnel simulation of an object in a moving fluid.  A tilted ellipse was chosen to generate asymmetries in the flow, leading to vortex shedding earlier in the simulation.  



§ ACKNOWLEDGMENTS

We thank Pavel Etingof, Slava Gerovitch, and Tanya Khovanova for their work on the MIT PRIMES program, which allowed us to collaborate together.  In March 2017, the first author was awarded the second prize at the Regeneron Science Talent search worth $175,000 based on this work that will cover the tuition fees at MIT starting in August 2017. We are sincerely humbled by Regeneron's very generous support. We thank Grady Wright for giving us advice on the implementation of the Navier–Stokes simulation.  We also thank Dan Fortunato and Heather Wilber.  This work was supported by National Science Foundation grant No. 1645445.

3
Apel_99_01 T. Apel, Anisotropic finite elements: local estimates and applications, vol. 3, Stuttgart, Teubner, 1999.
Babuska_76_01 I. Babuška and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal., 13 (1976), pp. 214–226.
evans L. C. Evans, Partial differential equations, American Mathematical Society, 2010.
Gilbarg_15_01 D. Gilbarg and N. S. Trudinger, Elliptic partial differential equations of second order, Springer-Verlag Berlin Heidelberg, reprint of 1998 edition, 2015.
higham N. Higham, Accuracy and stability of numerical algorithms, Second edition, SIAM, 2002.
Mathew T. Mathew, Domain decomposition methods for the numerical solution of partial differential equations, Springer Science & Business Media, 2008.
nseqs J. G. Liu, Projection methods I: Convergence and numerical boundary layers, SIAM J. Numer. Anal., 32 (1995), pp. 1017–1057.
NIST:DLMF F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Digital Library of Mathematical Functions, <http://dlmf.nist.gov>
Olver_13_01 S. Olver and A. Townsend, A fast and well-conditioned spectral method, SIAM Review, 55 (2013), pp. 462–489.
Patera_84_01 A. T. Patera, A spectral element method for fluid dynamics: laminar flow in a channel expansion, J. Comput. Phys., 54 (1984), pp. 468–488.
Ruppert_95_01 J. Ruppert, A Delaunay refinement algorithm for quality 2-dimensional mesh generation., J. Algor., 18 (1995), pp. 548–585.
Schoen_08_01 J. Schoen, Robust, guaranteed-quality anisotropic mesh generation, PhD Dissertation, Department of Electrical Engineering and Computer Sciences, University of California at Berkeley, 2008.
Shewchuk_02_01 J. R. Shewchuk, What is a good linear finite element? interpolation, conditioning, anisotropy, and quality measures, University of California at Berkeley, 73, 2002.
Shewchuk_02_02 J. R. Shewchuk, Delaunay refinement algorithms for triangular mesh generation., Computational Geometry, 22 (2002), pp. 21–74.
Solin_03_01 P. Šolín, K. Segeth, and I. Doležel, Higher-order finite element methods, Chapman & Hall/CRC Press, 2003.
townsend_computing_2014 A. Townsend, Computing with functions in two dimensions,PhD Dissertation Department of Mathematics, Oxford University, 2014.
townsend_automatic_2015 A. Townsend and S. Olver, The automatic solution of partial differential equations using a global spectral method, J. Comput. Phys., 299 (2015), pp. 106–123.
trefethen L. N. Trefethen, Approximation theory and approximation practice, SIAM, 2013.




§ DETERMINANTS

Let (x_1, y_1), (x_2, y_2), and (x_3, y_3) be the vertices of an arbitrary triangle in counterclockwise order.
By the shoelace formula, the area of the triangle is

    A = 1/2[(x_1y_2 + x_2y_3 + x_3y_1) - (x_2y_1 + x_3y_2 + x_1y_3)].

We can construct three quadrilaterals from the triangle by partitioning the triangle along the line segments connecting the centroid of the triangle to the midpoints of the sides.  Without loss of generality, let our quadrilateral contain vertex 1 of the triangle, so it has the following vertices in counterclockwise order:

    (x_1, y_1),   (x_1 + x_2/2, y_1 + y_2/2),   (x_1 + x_2 + x_3/3, y_1 + y_2 + y_3/3),   (x_1 + x_3/2, y_1 + y_3/2).

Using equations (<ref>) and (<ref>), we have

    x = 1/24[(14 x_1 + 5 x_2 + 5 x_3) + (4 x_1 - 5 x_2 + x_3)r + (4 x_1 + x_2 - 5 x_3)s + (2 x_1 - x_2 - x_3)rs ],

and similarly for y.
We can use equations (<ref>) and (<ref>) to find that

    (r, s) = (b_1 c_2 - b_2 c_1) + (b_1 d_2 - b_2 d_1)r + (c_2 d_1 - c_1 d_2) s.

Plugging in values for b_1, c_1, d_1, b_2, c_2, and d_2, we have

    (r, s) = 1/576[   ( (4 x_1 - 5 x_2 + x_3) (4 y_1 + y_2 - 5 y_3) - (4 y_1 - 5 y_2 + y_3) (4 x_1 + x_2 - 5 x_3) )
    
     +    ((4 x_1 - 5 x_2 + x_3) (2 y_1 - y_2 - y_3) - (4 y_1 - 5 y_2 + y_3) (2 x_1 - x_2 - x_3) )r
    
     +    ((4 y_1 + y_2 - 5 y_3) (2 x_1 - x_2 - x_3) - (4 x_1 + x_2 - 5 x_3) (2 y_1 - y_2 - y_3) )s ].

We can note that

    (∑_i=1^n v_i x_i )(∑_i=1^n w_i y_i ) - (∑_i=1^n w_i x_i ) (∑_i=1^n v_i y_i ) = ∑_1 ≤ i < j ≤ n (x_i y_j - x_j y_i)(v_i w_j - v_j w_i),

and use this identity to simplify (<ref>) to obtain

    (r, s) = 1/96[(x_1y_2 + x_2y_3 + x_3y_1) - (x_2y_1 + x_3y_2 + x_1y_3)](4 + r + s).

Substituting in (<ref>) and letting A be the area of the triangle gives the desired result of

    (r, s) = A(4 + r + s)/48.


