A nonconforming primal hybrid finite element method for the two-dimensional vector Laplacian

We introduce a nonconforming hybrid finite element method for the two-dimensional vector Laplacian, based on a primal variational principle for which conforming methods are known to be inconsistent. Consistency is ensured using penalty terms similar to those used to stabilize hybridizable discontinuous Galerkin (HDG) methods, with a carefully chosen penalty parameter due to Brenner, Li, and Sung [Math. Comp., 76 (2007), pp. 573-595]. Our method accommodates elements of arbitrarily high order and, like HDG methods, it may be implemented efficiently using static condensation. The lowest-order case recovers the $P_1$-nonconforming method of Brenner, Cui, Li, and Sung [Numer. Math., 109 (2008), pp. 509-533], and we show that higher-order convergence is achieved under appropriate regularity assumptions. The analysis makes novel use of a family of weighted Sobolev spaces, due to Kondrat'ev, for domains admitting corner singularities.

It is well established that conforming finite element methods for (1) have severe difficulties.For instance, a finite element approximation that is both div-and curl-conforming will also be H 1 -conforming, but when Ω is non-convex, H 1 (Ω) of edges, partitioned into interior edges E • h and boundary edges E ∂ h .Denote the broken L 2 inner products (•, •) •⟩ e .The method of [6] is based on the variational problem where γ| e := γ e > 0 is a penalty parameter on each e ∈ E h , to be detailed further in Section 2, and • is the jump in both tangential and normal components across an interior edge.(See also Brenner, Li, and Sung [9,10,11,12] for related work on curl-curl source problems and eigenproblems arising in Maxwell's equations.)In [6], u h and v h are linear vector fields continuous at the midpoint of each e ∈ E • h (i.e., both components live in the P 1 -nonconforming space of Crouzeix and Raviart [17]) whose tangential components vanish at the midpoint of each e ∈ E ∂ h .Brenner and Sung [13] later developed a quadratic nonconforming element for this problem and conjectured that it could be generalized to higher degree, as well as to dimension three.The two-dimensional conjecture was subsequently proved by Mirebeau [27], who also gave a counterexample to the three-dimensional case.However, for k > 1, the order-k elements are not simply P k vector fields: they are enriched by additional vector fields up to degree 2k − 1 that are gradients of harmonic polynomials.
In this paper, we present a three-field primal hybridization of (3) in the following form: Find (u h , p h , ûh ) ∈ V h × Q h × Vh such that, for all (v h , q h , vh ) where ph := p h + γ(u h − ûh ) and ⟨•, •⟩ ∂T h := K∈T h ⟨•, •⟩ ∂K .With appropriately chosen finite element spaces, as detailed in Section 2, this method has the following properties: • The lowest-order case is a hybridization of the method of Brenner et al. [6].
• Arbitrarily high order may be obtained using standard polynomial finite elements.The more exotic Brenner-Sung-Mirebeau spaces and projections play a crucial role in the analysis but are not needed for implementation.• As with HDG methods [14], the hybrid formulation enables efficient local assembly and static condensation, where u h and p h may be eliminated to solve a smaller global system involving only the approximate trace ûh .In addition to these contributions, we also present a novel error analysis using weighted Sobolev spaces, cf.Costabel and Dauge [16].This approach allows us to obtain error estimates on domains admitting corner singularities, without imposing the mesh-grading conditions on T h required by [6].
The paper is organized as follows.In Section 2, we describe the method and discuss its fundamental properties.Next, in Section 3, we present the error analysis of the method.Finally, in Section 4, we present the results of numerical experiments, which demonstrate and confirm the analytically obtained convergence results.impose inter-element continuity and boundary conditions, we define Vh := vh ∈ e∈E h P 2k−1 (e) 2 : vh × n = 0 on E ∂ h , whose elements are single-valued on edges.The extra variables p h and ûh , and their role in the variational problem, may be understood as follows.From (4b), we see that p h acts as a Lagrange multiplier, constraining the degree ≤ k − 1 moments of u h and ûh to agree on E h .Consequently, u h satisfies weak inter-element continuity and boundary conditions, and ûh may be seen as an approximate trace of u.Next, on each K ∈ T h , taking the inner product of the strong form (2a) with v ∈ H(div; K) ∩ H(curl; K) and integrating by parts implies that the solution to (1) satisfies

Comparing with (4a) and writing ⟨p
it follows that ph • n| ∂K and ph × n| ∂K can be seen as approximating −∇ • u| ∂K and ∇ × u| ∂K , respectively.Finally, (4c) shows that ûh also acts as a Lagrange multiplier, constraining ph • n and ph × n to be single-valued on interior edges and ph • n = 0 on boundary edges.The latter may be seen as an approximation of the natural boundary condition (2c).
To ensure convergence of the method for solutions with corner singularities, the penalty γ must be chosen carefully.Here, we recall the penalty used by Brenner et al. [6], which is the same one that we will use.Denote the corners of Ω by c 1 , . . ., c L , and let r ℓ (x) := |x − c ℓ | be the distance from x ∈ Ω to each corner.Given a multi-exponent λ = (λ 1 , . . ., λ L ), we denote r λ := L ℓ=1 r λ ℓ ℓ .Now, at each corner c ℓ , with interior angle ω ℓ , choose a parameter µ ℓ such that and µ := (µ 1 , . . ., µ L ).For each e ∈ E h , whose midpoint is denoted m e , we then define Finally, the penalty parameter on e is taken to be where |e| is the length of e.This ensures that γ e ∼ 1/|e| away from corners, while being appropriately weakened near corners to allow convergence to singular solutions, as we will see in Section 3.

2.2.
The Brenner-Sung-Mirebeau element and projection.We now recall the nonconforming finite element developed in Brenner and Sung [13] and Mirebeau [27], which we call the Brenner-Sung-Mirebeau (BSM) element.While it is not used to implement the method described above, this element and its associated projection play an important role in the numerical analysis of the method-and will also make clear why we have taken polynomial spaces of degrees 2k − 1 and k − 1.
Definition 2.1.Given a positive integer k, define the Brenner-Sung-Mirebeau (BSM) space on a triangle K ⊂ R 2 to be where H j (K) is the space of homogeneous harmonic polynomials of degree j on K. (By harmonic, we mean having vanishing Laplacian.) We immediately see that 2 , with equality if and only if k = 1.Indeed, since dim H j (K) = 2 for each j, it follows that dim BSM k (K) = k(k + 5).Brenner and Sung [13] conjectured, and Mirebeau [27] proved, that an element v h ∈ BSM k (K) is uniquely determined by the k(k + 5) degrees of freedom Moreover, the canonical interpolation using these degrees of freedom naturally defines a projection for all q h and w h as above.Letting P h : L 2 (K) → P k−1 (K) be the L 2 -orthogonal projection for scalar fields, we obtain the following commuting-projection property; the proof is basically identical to that in Brenner and Sung [13].
Proof.For all ϕ h ∈ P k−1 (K), integrating by parts using the divergence theorem gives with q h = ϕ h n and w h = ∇ϕ h .This proves the first equality; the proof of the second is essentially the same, using Green's theorem instead of the divergence theorem.□ The solution to (1) satisfies the hypotheses of this lemma on each K ∈ T h , as a result of the regularity theory discussed in Section 3.1.

2.3.
Equivalence to reduced methods with jump terms.We next show that the three-field hybrid method described in Section 2.1 may be reduced to a two-field or one-field method with jump terms.The coupling introduced by the jump terms prevents static condensation, so we generally prefer the three-field formulation for implementation.However, these reduced formulations will be useful analytically, and will help in relating our method to that of Brenner et al. [6]. 1irst, we introduce notation and definitions for the average and jump of a vector field across interior edges.Suppose e ∈ E • h is an interior edge shared by two triangles, K + and K − , and let n ± denote the unit normal to e pointing outward from K ± .If a vector field w takes values w ± on the K ± sides of e, we define the average and jump of w at e to be where w ⊗ n := wn ⊤ is the outer product.It is straightforward to see that the i-th row of w e is the transpose of the usual scalar jump w i e = w + i n + + w − i n − for i = 1, 2. This definition of • encodes the jump in both tangential and normal directions, without requiring a global orientation of the edges.It is then easily verified that the ⟨•, •⟩ ∂T h inner product of vector fields (just as for scalar fields) may be expanded as (7) ⟨w, where the inner product of the matrix-valued jumps is taken in the Frobenius sense.Since functions are single-valued on boundary edges, we leave average and jump undefined on Taking this as the test function in (4b) and applying the identity (7) gives There are no interior jump-jump or tangential boundary terms, since this choice of q h has q h = 0 on E • h and q h × n = 0 on E ∂ h .Similarly, vh = 0 on E • h and vh × n = 0 on E ∂ h for all vh ∈ Vh , so (4c) may be rewritten as (9) 2 The terms involving p h vanish by (8), leaving Finally, substituting these equalities into (9) gives Note that { {w} } e = 0 can be rewritten as w + = −w − .Since the outer normals satisfy n + = −n − , it follows that w + ×n + = w − ×n − and w + •n + = w − •n − , i.e., the tangential and normal components of w agree on both sides of e.Thus, Lemma 2.3 says that the tangential and normal components of p h and u h − ûh are single-valued, with normal components vanishing on boundary edges.In particular, the same is therefore true of ph , as previously remarked in Section 2.1.
Using Lemma 2.3 and the identity (7), observe that the edge terms in (4a) reduce to Similarly, the edge terms in (4b) reduce to for all q h ∈ Q h .This allows us to eliminate ûh and the equation (4c) from the variational problem.A two-field reduced formulation is defined as follows.Let and define the bilinear forms a h : We then consider the problem: This resembles a standard two-field hybrid method in saddle-point form, where Qh is the space of Lagrange multipliers.Compare with the nonconforming hybrid method of Raviart and Thomas [31] for the scalar Poisson equation.
Finally, we may reduce even further to a one-field formulation on consisting of degree-(2k − 1) vector fields whose degree ≤ k − 1 moments are continuous on E • h and have vanishing tangential component on E ∂ h .We then consider the problem: Find This is precisely (3), modulo a constant factor of 1 2 for the penalty on interior edges, and the lowest-order case k = 1 recovers the method of Brenner et al. [6].
We have thus shown the equivalence of the three-field, two-field, and one-field formulations, which we now state as a lemma.
Lemma 2.5.The following are equivalent: and ûh is as in (ii).

2.4.
Existence/uniqueness and static condensation.The problem (1) is well-posed if and only if α is not an eigenvalue of the vector Laplacian on H(div; Ω) ∩ H(curl; Ω).In particular, the bilinear form and if the complement of Ω is connected (e.g., Ω is simply connected), then it is also coercive for α = 0 by Friedrichs's inequality (cf.Monk [28]).
We are now ready to prove our main result on existence and uniqueness for the hybrid method.
Theorem 2.6.For the problems (4), (10), and (11), existence and uniqueness of solutions holds-or fails to hold-simultaneously for all three.In particular, all three are uniquely solvable if α > 0, and if the complement of Ω is connected, then this is also true for α = 0.
Proof.By Lemma 2.5, unique solvability of ( 4) is equivalent to that of (10), since ûh is uniquely determined by u h , so it suffices to show equivalence of ( 10) and (11).Using classic saddle-point theory surjective and the restriction of a h (•, •) to its kernel is an isomorphism.The isomorphism-on-thekernel condition is precisely the unique solvability of (11), so it remains to show that the surjectivity condition holds.
In fact, we will show something slightly stronger, which is surjectivity of the map B h : □ Next, we discuss the static condensation of the hybrid method, which eliminates the spaces V h and Q h from (4) to obtain a smaller global variational problem on Vh alone.We take a similar approach to that used for HDG methods in Cockburn, Gopalakrishnan, and Lazarov [14].Observe that, given ûh and f , (4a To separate the influence of ûh and f , we define two local solvers: Lemma 2.7.If α ≥ 0, then the local solvers are well-defined, i.e., (12) is uniquely solvable.
Proof.First, we show that K∈T h a K (•, •) is coercive.This is obvious when α > 0; when α = 0, a K (u h , u h ) = 0 implies that u| ∂K = 0, so Friedrichs's inequality implies u| K = 0. Finally, the surjectivity of already been shown in the proof of Theorem 2.6.□ Assuming the local solvers are well-defined-which always holds for α ≥ 0, by Lemma 2.7-we now define Pû h := Pû h + γ(Uû h − ûh ) and Pf := Pf + γUf .Substituting into (4c) and rearranging gives the condensed problem: Find ûh ∈ Vh such that, for all vh ∈ Vh , Since the local solvers may be computed element-by-element, in parallel if desired, the condensation from (4) to ( 13) is efficient to implement.The condensed bilinear form âh (û h , vh ) := −⟨ Pû h , vh ⟩ ∂T h on the left-hand side of (13) has the following useful symmetric expression.
Lemma 2.8.For all ûh , vh ∈ Vh , Proof.We begin by writing For the first term, (12a For the second term, (12b) implies ⟨Pû h , Uv h − vh ⟩ ∂T h = 0, so which completes the proof.□ Theorem 2.9.Assuming the local solvers are well-defined, (u h , p h , ûh ) ∈ V h × Q h × Vh is a solution of (4) if and only if ûh is a solution of (13) with u h = Uû h + Uf and p h = Pû h + Pf .Consequently, (4) is uniquely solvable if and only if (13) is.In particular, âh (•, •) is symmetric positive-definite if α > 0, and if the complement of Ω is connected, then this is also true for α = 0.
Proof.The equivalence of ( 4) and ( 13) has already been demonstrated in the discussion above.When α ≥ 0, Lemma 2.7 states that the local solvers are well-defined, and Lemma 2.8 implies that âh (•, •) is symmetric positive-semidefinite.Furthermore, if âh (û h , ûh ) = 0, then Uû h = ûh on E h , so Uû h ∈ H(div; Ω) ∩ H(curl; Ω) with a(Uû h , Uû h ) = 0. Hence, as in the proof of Theorem 2.6, âh (•, •) is positive-definite whenever a(•, •) is.□ Remark 2.10.These results tell us that static condensation from ( 4) to (13) does not merely reduce the size of the global system.It also makes the system more amenable to efficient global solvers, such as the conjugate gradient method in the case where ( 13) is positive-definite.

Regularity and error analysis
3.1.Weighted Sobolev spaces and regularity.Costabel and Dauge [16] characterize the regularity of solutions to Maxwell's equations in two dimensions (as well as in three) using a family of weighted Sobolev spaces due to Kondrat'ev [23].We now recall these spaces and give corresponding regularity results for the problem (1), combining the approach used in [16] with that of Brenner et al. [6,Section 2].For detailed treatments of Kondrat'ev spaces and elliptic regularity in domains with corners, we refer the reader to Nazarov and Plamenevsky [29] and Kozlov, Maz'ya, and Rossmann [24].As in Section 2.1, let r ℓ (x) denote the distance from x ∈ Ω to a corner c ℓ and r λ := L ℓ=1 r λ ℓ ℓ for a multi-exponent λ = (λ 1 , . . ., λ L ).Given a nonnegative integer m, define the weighted Sobolev space for all |β| ≤ m , where β is a multi-index, equipped with the natural norm defined by This space also has the following equivalent characterization: If Ω = Ω 0 ∪ L ℓ=1 Ω ℓ , where Ω 0 contains none of the corners and Ω ℓ contains only corner c ℓ , then From the definitions, we immediately obtain the continuous inclusion , which may be interpolated to obtain fractional-order spaces.That is, if s ≥ 0, then V s λ (Ω) may be defined by complex interpolation between V V s+ϵ λ+ϵ (Ω) ⊂ V s λ (Ω), for all s ≥ 0 and ϵ > 0. Additionally, the continuous inclusions extend in the obvious way from nonnegative integer m to real s ≥ 0. Remark 3.1.Schneider [33] uses an alternative notation for Kondrat'ev spaces, For example, the inclusion K m+1 p,a (Ω) ⊂ K m p,a (Ω) in the notation of [33] gives V m+1 λ+1 (Ω) ⊂ V m λ (Ω) in our notation, since p = 2 and a = m − λ = (m + 1) − (λ + 1).Fractional Kondrat'ev spaces are denoted in [33] by K s p,a (Ω), and similarly we have V s λ (Ω) = K s 2,s−λ (Ω).Finally, we note that an intrinsic treatment of fractional weighted Sobolev spaces may be found in Dauge [18,Appendix A].
Suppose now that u ∈ H(div; Ω) ∩ H(curl; Ω) satisfies (1).We recall that ∇ • u ∈ H1 (Ω), since it can be seen as the solution to the Dirichlet problem Assuming Ω is simply connected, we can express u in terms of its Helmholtz decomposition u = ∇ϕ + ∇ × ψ, where ϕ ∈ H1 (Ω) and ψ ∈ H 1 (Ω) solve with homogeneous Dirichlet and Neumann boundary conditions, respectively.For uniqueness, we take Ω ψ = 0. To determine the regularity of ϕ and ψ, we follow Chapter 2 of Nazarov and Plamenevsky [29], which characterizes the regularity of solutions to Dirichlet and Neumann problems in plane domains with corner points; similar results are also found in Kozlov et al. [24, §6.6.1-6.6.2].
) is well-posed, then we have the stability estimate Proof.Pick ϵ > 0 such that s < µ ℓ − ϵ for all ℓ.Then this follows by the continuous inclusions □ Finally, we note that this also implies the following well-known unweighted Sobolev regularity result, cf.Assous, Ciarlet, and Sonnendrücker [2].
is well-posed, then we have the stability estimate ∥u∥ 2s ≲ ∥f ∥ Ω .
2 , and we have the continuous inclusions In particular, since ω ℓ < 2π for all ℓ, we may take s > 1 4 in Corollary 3.4 to conclude that u ∈ H σ (Ω) 2 with σ > 1 2 .3.2.Preliminary estimates.We now establish two weighted Sobolev norm approximation results that will be useful in the subsequent error analysis; compare Lemmas 5.2 and 5.3 in [9].
For the remainder of the paper, we assume that T h is shape-regular, but we make no additional assumptions about quasi-uniformity or grading.Let h K denote the diameter of K ∈ T h and h := max K∈T h h K .We denote the weighted Sobolev norm on V s λ (Ω)| K by ∥•∥ s,λ,K (with distances taken to the corners of Ω, not those of K) and the ordinary Sobolev seminorm on for all K ∈ T h and e ⊂ ∂K.
Proof.If K does not have any of the corners c ℓ as a vertex, then v| K ∈ H s+1 (K) 2 , so the trace inequality with scaling and Bramble-Hilbert lemma imply . By shape regularity, we have Φ µ (e) = r 1−µ (m e ) ∼ r 1−µ (x) for all x ∈ K, and therefore On the other hand, if K has c ℓ as a vertex, then the inclusions Hence, the trace inequality with scaling and Bramble-Hilbert give , and therefore where the last inequality is due to the continuity of the inclusion for all K ∈ T h and e ⊂ ∂K.
Thus, for all K ∈ T h and e ⊂ ∂K, the trace inequality with scaling and Bramble-Hilbert lemma give By shape regularity, Φ µ (e) −1 = r µ−1 (m e ) ≲ r µ−1 (x) for all x ∈ K, and therefore which completes the proof.□ 3.3.Error estimates.We now estimate the error u − u h , where u satisfies (1) and u h satisfies (11).The argument follows a similar general outline to that in Brenner et al. [6], but the details differ in several important respects-especially in the use of weighted Sobolev regularity hypotheses and higher-order polynomial approximation, and in the absence of mesh-grading assumptions.
As in [6], we will first estimate the error in the mesh-dependent energy norm If we extend a h (•, •) from Vh to H(div; Ω) ∩ H(curl; Ω) + Vh , then in the special case α = 1, this is precisely the norm associated to a h (•, •) considered as an inner product.
For arbitrary α, we immediately see that a h (•, •) is bounded with respect to ∥•∥ h .For α > 0, we have the coercivity condition a h (v, v) ≥ min(1, α)∥v∥ 2 h .If the complement of Ω is connected, then we also have coercivity for α = 0, by the argument in the proof of Theorem 2.6.In general, for α ≤ 0, we have a Gårding inequality (which is actually an equality), h .This implies the following Strang-type abstract estimates, whose proofs are identical to those of Lemma 3.5 and Lemma 3.6 in Brenner et al. [10].Lemma 3.7.If α > 0, u is the solution to (1), and u h is the solution to (11), then and if the complement of Ω is connected, then this also holds for α = 0.If α ≤ 0, u satisfies (1), and u h satisfies (11), then We will proceed by estimating the two terms on the right-hand side of (15), which correspond to approximation error and consistency error, respectively.
Proof.The first inequality holds since the BSM projection maps u ∈ H σ (Ω) 2 with σ > 1 2 to Π h u ∈ Vh .Next, letting µ min := min ℓ µ ℓ , the continuous inclusion 2 , so polynomial approximation theory gives It remains to estimate the contributions from the penalty terms.Observe that, for e ∈ E • h , by the parallelogram identity.Therefore, so it suffices to estimate the contribution from each K ∈ T h and e ⊂ ∂K.By Lemma 3.5, (17), (18), and (19) and summing over K ∈ T h completes the proof.□ Lemma 3.9.Suppose u satisfies (1) and u h satisfies (11).
Proof.Subtracting ( 11) from ( 5) with v = v h = w h ∈ Vh , we get where we have denoted the normal and tangential jump components on e ∈ E • h by The condition w h ∈ Vh says that w h • n e and w h × n e are each L 2 -orthogonal to P k−1 (e) for e ∈ E • h , and that w h × n| e is L 2 -orthogonal to P k−1 (e) for e ∈ E ∂ h .Therefore, letting P h be projection onto either triangle where the last step uses the Cauchy-Schwarz inequality.Applying Lemma 3.6 with η = ∇ • u to the first term and the definition of the energy norm to the second, we conclude that and the result follows.□ Next, we use a duality argument to control the error in the L 2 norm.
Lemma 3.10.Suppose u satisfies (1) and u h satisfies (11).Suppose also that ∇•u, ∇×u ∈ V s µ−1 (Ω) with s ≤ k, and let t < µ min .If (1) is well-posed, then , and we have the stability estimate ( 22) Hence, Lemma 3.8 implies To express ∥u − u h ∥ 2 Ω in terms of z, we would like to take v = u − u h in ( 21), but we cannot do so since generally u h / ∈ H(div; Ω) ∩ H(curl; Ω).Instead, integrating by parts as in (5) gives which we will estimate term-by-term.
For the first term of ( 24), we write By the boundedness of a h (•, •) in the energy norm and ( 23), we have Next, by (20) with w h = Π h z ∈ Vh , we have By a similar argument to that used in Lemma 3.9, along with the fact that z where the last two lines use Lemma 3.6 with η = ∇ • u and (23).Similarly, Thus, we have estimated the first term of ( 24) by ( 25) For the remaining terms of ( 24), we use a similar argument to the one above to get where the last two lines use Lemma 3.6 with η = ∇ • z and (22).Likewise, Altogether, estimating (24) by combining ( 25), (26), and ( 27), we have Finally, we are ready to state the main energy and L 2 error estimates.
Theorem 3.11.Suppose u satisfies (1) , where s ≤ k, and let t < µ min .If α > 0, then the solution u h to (11) satisfies the error estimates and if the complement of Ω is connected, then this also holds for α = 0.If α < 0 and (1) is well-posed, then (11) is uniquely solvable for sufficiently small h, and the solution u h satisfies these same estimates.
Proof.If α > 0, or if α = 0 with Ω having connected complement, then the proof is fairly immediate.The energy estimate follows from the abstract estimate (15) in Lemma 3.7, together with Lemmas 3.8 and 3.9, and the L 2 estimate follows by Lemma 3.10.
When α < 0 is such that (1) is well-posed, we follow the approach in Brenner et al. [10, Theorem 4.5], which uses a technique for indefinite problems due to Schatz [32].Suppose that u h satisfies (11).From the abstract estimate ( 16) in Lemma 3.7, along with Lemmas 3.8, 3.9, and 3.10, we have where the constant C has been made explicit.Now, choose h * small enough that Ch t * < 1.It follows that, whenever h ≤ h * , we may subtract Ch t ∥u − u h ∥ h from both sides of (28) to obtain In particular, when f = 0, well-posedness of (1) gives the unique solution u = 0, and ∥u h ∥ h ≲ 0 implies that (11) has the unique solution u h = 0. Hence, (11) is uniquely solvable and satisfies the energy estimate whenever h ≤ h * , and the L 2 estimate follows by another application of Lemma 3.10.1. Convergence to a smooth solution on a square domain.
Corollary 3.12 (minimum-regularity case).If α > 0, u is the solution to (1), and u h is the solution to (11), then for all s < µ min , we have the error estimates and if the complement of Ω is connected, then this also holds for α = 0.If α < 0 and (1) is well-posed, then these estimates hold for sufficiently small h.
Proof.This is immediate from Corollary 3.3 and Theorem 3.11 with s = t.□

Numerical experiments
In this section, we present numerical experiments illustrating the convergence behavior of the method, showing how convergence is affected by the interior angles of Ω and by the regularity of the exact solution u, and relating these numerical results to the theoretical results of Section 3.For all numerical experiments, we take α = 1.
All computations have been carried out using the Firedrake finite element library [30] (version 0.13.0+4959.gac22e4c5),and a Firedrake component called Slate [19] was used to implement the local solvers for static condensation and postprocessing.Since all four corners of Ω are π/2, we have µ = (1, 1, 1, 1).Given N ∈ N, we construct a uniform triangle mesh by partitioning Ω uniformly into N × N squares, then dividing each into two triangles.
Table 1 shows the result of applying our method to the problem whose exact solution is (This is the same u that Brenner et al. [6] use for their numerical experiments on the square.)Since u is smooth, we observe convergence rates of k for the energy error and k + 1 for the L 2 error, consistent with Theorem 3.11.
Given N ∈ N, we construct a uniform triangle mesh of Ω by taking a uniform 2N × 2N mesh of the square − 1 2 , 1 2 2 , as in Section 4.1, and removing the first quadrant.Table 2. Convergence to the minimum-regularity singular harmonic on an L-shaped domain.
which is a harmonic vector field with ∇ • u = 0 and ∇ × u = 0. We observe that u ∈ V m m−1/3−µ (Ω) for all m, since the condition for ∂ β u to be in the appropriate weighted L 2 space in a δ-neighborhood of the origin is Table 2 shows the results of applying our method to this problem, where the inhomogeneous boundary conditions are imposed on ûh × n by interpolating u × n on E ∂ h .Since s = 1/3, we observe minimal convergence rates of approximately 1/3 for the energy error and 2/3 for the L 2 error for all k, consistent with Theorem 3.11.4.2.2.Minimum-regularity nonsingular vector field.The next example shows that even a nonsingular vector field may have minimum regularity, owing to the conditions ∇ • u, ∇ × u ∈ V s µ−1 (Ω).Given arbitrarily small ϵ > 0, consider the problem whose exact solution is where the inhomogeneous boundary conditions may also be dealt with as in Remark 4.1.By a similar calculation as for the singular harmonic vector field, we have u ∈ V m m−5/3−µ (Ω) 2 for all m, and thus u ∈ V 5/3+1 1−µ (Ω) 2 .However, we merely have µ−1 (Ω).Hence, even though u does not have a singularity at the origin, the regularity hypotheses of Theorem 3.11 hold merely with s = 1/3.
Table 3 shows the results of applying our method to this problem with ϵ = 0.001.As with the previous example, since s = 1/3, we observe minimal convergence rates of approximately 1/3 for the energy error and 2/3 for the L 2 error for all k, consistent with Theorem 3.11.

4.3.
Higher-regularity solutions on L-shaped domain.Finally, we present numerical results for convergence to solutions with higher regularity on the L-shaped domain, observing improved convergence for larger k.As in Section 4.2, we consider both a harmonic and a non-harmonic example-here, both having s = 7/3-on the same family of uniform meshes, where inhomogeneous boundary conditions for u × n on ∂Ω are handled in the same way.4.3.1.Higher-regularity harmonic vector field.Consider the problem whose exact solution is which is a harmonic vector field with ∇ • u = 0 and ∇ × u = 0.By a similar calculation to that in Section 4.2.1, we get u ∈ V m m−7/3−µ (Ω) 2 for all m.By interpolation, we see that u ∈ so the hypotheses of the error estimates in Section 3.3 hold with s = 7/3.Table 4 shows the results of applying our method to this problem.Since s ≤ 3, for k = 3 we observe the maximum convergence rates predicted by Theorem 3.11: roughly 7/3 for the energy error and 8/3 for the L 2 error.For k = 2, however, we also observe rates of approximately 7/3 for the energy error and 8/3 for the L 2 error.This is explained by the fact that u is the curl of a harmonic function, and BSM k (K) contains gradients (hence curls) of harmonic polynomials with degree ≤ 2k.In this special case, the condition s ≤ k in the approximation estimate Lemma 3.8 improves to s ≤ 2k − 1, while the consistency error in Lemma 3.9 vanishes due to ∇ • u = 0 and ∇ × u = 0.For k = 1, we observe the expected energy-norm convergence rate of 1, but the L 2 -norm convergence rate of 2 is better than the duality-based estimate of 4/3 in Theorem 3.11.We do not yet have a satisfying analytical explanation for this better-than-expected gap between the energy-norm and L 2 -norm rates when s > k; see further discussion in the next example and in Section 5.
Table 5 shows the results of applying our method to this problem with ϵ = 0.001.For all k, we observe a convergence rate of approximately min(k, 7/3) in the energy norm, consistent with Theorem 3.11.This also supports the argument that the improved energy error in Section 4.3.1, which is not observed here, was due to that exact solution being the curl of a harmonic function.For k = 3, we observe the expected L 2 -norm convergence rate of approximately 8/3.However, for k = 2 and k = 1, we observe better-than-expected rates of 8/3 (rather than 7/3) and 2 (rather than 4/3), respectively, similar to what we saw with the k = 1 case in Section 4.3.1.

Conclusion
We have presented a nonconforming primal hybrid finite element method for the two-dimensional vector Laplacian that extends the P 1 -nonconforming method of Brenner et al. [6] to arbitrary order k.The method uses only standard polynomial finite elements, although the more exotic BSM element and projection play a key role in the analysis, and the method may be implemented efficiently using static condensation.Using the weighted Sobolev spaces of Kondrat'ev for domains with corners, we have obtained error estimates that hold on general shape-regular meshes, without mesh-grading conditions.These estimates establish the convergence of the method, even for minimum-regularity solutions with corner singularities, and the convergence rate improves with k to the extent regularity allows.
Let us conclude with a brief discussion of one area where the numerical experiments in Section 4 suggest possible room for improvement.Dropping the hypothesis that s ≤ k, we may rewrite the estimates of Theorem 3.11 as ∥u − u h ∥ h ≲ h min(k,s) ∥u∥ s+1,1−µ + ∥∇ • u∥ s,µ−1 + ∥∇ × u∥ s,µ−1 , (29) ∥u − u h ∥ Ω ≲ h min(k,s)+t ∥u∥ s+1,1−µ + ∥∇ • u∥ s,µ−1 + ∥∇ × u∥ s,µ−1 .(30) From the numerical experiments, it appears that the energy estimate ( 29) is sharp, and in general one cannot relax the restriction t < µ min in the L 2 estimate (30).However, when s > k, it appears that a sharper estimate than (30) holds, which we now state as a conjecture.Establishing this would require some new analytical arguments, perhaps involving weighted-norm error estimates.Indeed, a duality estimate of the sort in Lemma 3.10 can only give an L 2 error estimate of the form (30), where the improvement over the energy rate is the same for all k, and the numerical experiments suggest that there is no way to sharpen this uniformly in k.
Finally, it is natural to ask how the two-dimensional method presented in this paper might be generalized to three dimensions.In contrast with some other approaches that are limited to dimension two-such as methods that use the Hodge decomposition to transform vector problems into scalar problems [7,25,26,8]-the variational form of the hybrid method (4) extends naturally to the three-dimensional case.The main challenge is to choose suitable finite element spaces and a suitable penalty, and here there are two obstacles to overcome.First, when k > 1, we do not yet know a three-dimensional version of the BSM element and commuting projection, which would be needed to make the analysis work.(A naive extension to three dimensions fails to satisfy unisolvence when k = 2, as shown by Mirebeau [27].)Consequently, it is not clear what polynomial degrees would be needed for the finite element spaces V h , Q h , and Vh .Second, the weighted Sobolev analysis in three-dimensional domains becomes more complicated, since singularities can form along boundary edges, as well as at corners where edges meet, cf.Costabel and Dauge [16].Consequently, a penalty γ would need to be carefully constructed, likely involving the distances both to edges and to corners with suitable exponents.

4. 1 .
Smooth solution on square domain.We begin by considering the square domain Ω = 0,

δ 0 r − 1 / 3 −µ 1 +|β| r 2 / 3 − 1 −|β| 2 r dr = δ 0 r 2 ( 1 / 3 −µ 1 )
−1 dr < ∞, which holds since µ 1 < 1/3.(Compare Costabel and Dauge [16, Theorem 6.1].)By interpolation, we getu ∈ V 1/3+1 1−µ (Ω)2 , so the hypotheses of the error estimates in Section 3.3 hold with s = 1/3.Remark 4.1.Although u does not satisfy the homogeneous boundary condition u × n = 0 on all of ∂Ω, it does satisfy this condition on the boundary edges θ = π/2 and θ = 2π adjacent to the reentrant corner.Thus, taking ϕ to be a smooth cutoff function supported in a small neighborhood of the origin, we may write u = uϕ + u(1 − ϕ), where uϕ satisfies the homogeneous boundary condition and u(1 − ϕ) is a smooth extension of the inhomogeneous boundary condition.It follows that u and uϕ have the same regularity, and standard arguments may be used to extend the numerical properties of the method from the homogeneous boundary value problem with exact solution uϕ to the inhomogeneous boundary value problem with exact solution u.
since it can be seen as the zero-mean solution to a Neumann problem.See Costabel and Dauge [15, Theorem 1.2] and similar arguments in Brenner et al. [6, Section 2].Using this, we may now obtain a minimum weighted Sobolev regularity result for u itself.Theorem 3.2.If u satisfies (1), then u ∈ V 2 2−2µ+ϵ (Ω) 2 for all ϵ > 0. Furthermore, if (1) is well-posed, then we have the stability estimate ∥u∥ 2,2−2µ+ϵ ≲ ∥f ∥ Ω .Proof.As in Brenner et al. [6, Section 2], it is sufficient to establish regularity and stability for Ω simply connected, since the general case follows by a partition of unity argument.

Table 3 .
Convergence to a minimum-regularity nonsingular solution on an L-shaped domain.4.2.1.Minimum-regularity singular harmonic vector field.In polar coordinates (r, θ), we first consider the problem whose exact solution is