A Convenient Inclusion of Inhomogeneous Boundary Conditions in Minimal Residual Methods

1.1 MINRES Methods for Handling Inhomogeneous (Essential) Boundary Conditions

The possibility to handle inhomogeneous boundary conditions when solving PDEs is often mentioned as an advantage of minimal residual (MINRES) discretizations (e.g. [3, Chapter 12]). In most cases, however, it is not so clear how this can lead to satisfactory results.

Considering linear elliptic PDEs of second order, one option is to write the boundary value problem in an ultra-weak first-order system variational formulation, which renders all boundary conditions natural. The resulting “practical” MINRES method [16, §§ 4.4 & 5] yields a quasi-best approximation from the trial space to all variables with respect to L 2 -norms.

If one is interested in approximation with respect to other norms, then one can resort to first-order system weak or mild-weak variational formulations, or to the standard second-order variational formulation. In those cases, Dirichlet or Dirichlet and Neumann boundary conditions are essential ones. In the case that homogeneous versions of those boundary conditions lead to a well-posed variational problem, inhomogeneous ones can be appended by enforcing them weakly by introducing additional test spaces of functions defined on the boundary ∂ Ω of the domain Ω on which the PDE is posed. Then the combined variational formulation of the PDE and the boundary conditions can be shown to be well-posed. Indeed, for some Hilbert spaces 𝑋, Y 1 and Y 2 , let F : X → Y 1 ′ and T : X → Y 2 ′ be bounded, and with X 0 := ker ⁡ T , let F : X 0 → Y 1 ′ be boundedly invertible. In applications, 𝐹 and 𝑇 are operators corresponding to a weak imposition of the PDE and the boundary condition(s), and so 𝑋 and Y 1 are spaces of functions on Ω, and Y 2 is a (product of) space(s) of functions on ∂ Ω . Then, assuming 𝑇 is surjective, G := ( F , T ) : X → Y ′ := Y 1 ′ × Y 2 ′ is boundedly invertible [12, Lemma 2.7].

Given a finite-dimensional “trial” subspace X δ ⊂ X (“𝛿” stands for “discrete”), and a given forcing term and essential boundary data, the solution u ∈ X of the well-posed problem can be approximated by the residual minimizer u δ ∈ X δ in the norm of Y ′ . A complication is that some coordinate spaces of Y ′ will be negative – or, in particular for function spaces on ∂ Ω , fractional-order Sobolev spaces, all with norms that cannot be evaluated. Several solutions for this problem have been proposed, e.g. [4] for the case of negative-order Sobolev norms on Ω, and [18] for fractional-order Sobolev norms on ∂ Ω . The resulting MINRES methods, however, are not guaranteed to produce approximate solutions that are quasi-best.

Viewing all negative- or fractional-order Sobolev norms as dual norms, so writing e.g. H 1 2 ⁢ ( ∂ Ω ) as H − 1 2 ⁢ ( ∂ Ω ) ′ , our approach in [16] is to discretize dual norms by replacing the involved suprema over infinite-dimensional spaces by a supremum over a finite-dimensional (product) test space Y δ ⊂ Y . Under an inf-sup condition on the pair ( X δ , Y δ ) , the resulting MINRES approximation u δ ∈ X δ was shown to be quasi-best. The evaluation of the so-discretized dual norms is implemented by introducing the Riesz lift of the residual, viewed as an element of Y δ ′ , as an extra variable λ δ , which results in a saddle-point system for ( λ δ , u δ ) ∈ Y δ × X δ .

In the case that one or more components of 𝑌 are fractional Sobolev norms, a remaining issue is the evaluation of the scalar product between functions from Y δ . Without compromising quasi-optimality, it can be solved by replacing this scalar product on Y δ by an equivalent scalar product defined in terms of the inverse of an optimal preconditioner of linear complexity for the corresponding stiffness matrix. This construction permits then to eliminate the extra variable from the system, after which a symmetric positive definite system on X δ remains, with a system matrix that can be applied in linear complexity. Optimal preconditioners of linear complexity for fractional Sobolev spaces are available, e.g. BPX or multigrid for positive orders and [11, 19] for negative orders.

Finally, by applying an optimal preconditioner of linear complexity on the trial space X δ , the symmetric positive definite system can be solved iteratively within a tolerance of the order of the discretization error in linear complexity.

1.2 Current Paper

A disadvantage of the method from [16] is that its implementation is rather demanding, in particular because of the incorporation of the preconditioner(s) for fractional Sobolev norms on the boundary. In the current paper, we introduce an alternative approach, which can be expected to be computationally more costly, but that can be very easily implemented with a finite element package like NGSolve [17].

As in our approach from [16], given a trial space X δ ⊂ X , we select a test space Y δ ⊂ Y such that ( X δ , Y δ ) is (uniformly) inf-sup stable, replace the norm on Y ′ for the residual by the discretized norm on Y δ ′ and introduce its Riesz lift in Y δ as an extra variable resulting in a saddle-point problem.

We proceed differently to deal with the problem that also the scalar product on 𝑌 might not be evaluable. Recalling that G : X → Y ′ and thus G ′ : Y → X ′ is boundedly invertible, we equip 𝑌 with the equivalent norm ∥ G ′ ⋅ ∥ X ′ , known as the optimal test norm. The resulting dual norm on Y ′ then equals ∥ G − 1 ⋅ ∥ X so that the (unfeasible) exact residual minimization in this norm would yield the best approximation from X δ with respect to ∥ ⋅ ∥ X . Since also the new scalar product ⟨ G ′ ⋅ , G ′ ⋅ ⟩ X ′ on 𝑌 cannot be evaluated, for some finite-dimensional X ̂ δ ⊂ X , we replace it by a discretized version ⟨ G ′ ⋅ , G ′ ⋅ ⟩ X ̂ δ ⁣ ′ . Assuming that also ( Y δ , X ̂ δ ) satisfies a (uniform) inf-sup condition, the resulting MINRES approximation is quasi-best. To make this approximation computable, we introduce the Riesz lift θ δ ∈ X ̂ δ of G ′ ⁢ λ δ , viewed as an element of X ̂ δ ⁣ ′ , as a third variable, and so saddle the system once more.

The resulting 3 × 3 block system for the triple ( θ δ , λ δ , u δ ) ∈ X ̂ δ × Y δ × X δ has one block that corresponds to the stiffness matrix of X ̂ δ with respect to the “easy” scalar product ⟨ ⋅ , ⋅ ⟩ X , whereas all four other non-zero blocks correspond to system matrices of the bilinear form ( G ⋅ ) ( ⋅ ) with respect to X δ × Y δ or X ̂ δ × Y δ , or transposes of those blocks. So the “difficult” scalar product ⟨ ⋅ , ⋅ ⟩ Y vanished completely from the system, and the implementation is simple.

For suitable Y δ and X ̂ δ , we show that ∥ θ δ ∥ X is an efficient and asymptotically reliable a posteriori estimator for ∥ u − u δ ∥ X . We will use local 𝑋-norms of θ δ to drive an adaptive finite element method.

1.3 A Nonintrusive Approach

Although it is not subject of this paper, for completeness and comparison, in our abstract framework, we recall an alternative classical approach to deal with inhomogeneous essential boundary data (e.g. [10]).

Recall the setting of T : X → Y 2 ′ being bounded and surjective, F : X 0 := ker ⁡ T → Y 1 ′ being boundedly invertible and G := ( F , T ) : X → Y ′ := Y 1 ′ × Y 2 ′ being boundedly invertible. We will use that 𝑇 has a bounded right inverse E : Y 2 ′ → X (e.g. E := g ↦ G − 1 ⁢ ( 0 , g ) ).

Given ( f , g ) ∈ Y 1 ′ × Y 2 ′ , consider the problem G ⁢ u = ( f , g ) . Let X δ be a finite-dimensional subspace of 𝑋, and X 0 δ := X δ ∩ X 0 . Suppose that one has a method available that for g = 0 (i.e., “homogeneous essential boundary data”) produces a u δ ∈ X δ that is quasi-best, i.e., for some constant C > 0 , ∥ u − u δ ∥ X ≤ C ⁢ inf w δ ∈ X 0 δ ∥ u − w δ ∥ X . Using this method, one can approximate the solution 𝑢 for general 𝑔 as follows. First, construct g δ ∈ ran ⁡ T | X δ , i.e., g δ = T ⁢ v δ for some v δ ∈ X δ , that approximates 𝑔; and second, approximate with aforementioned method the solution z ∈ X 0 of G ⁢ z = ( f − F ⁢ v δ , 0 ) by z δ ∈ X 0 δ . Then

Furthermore, following [1, § 3.3], let us consider the case that there exists a uniformly bounded projector J δ : X → X δ that preserves “homogeneous essential boundary conditions”, i.e., T ⁢ J δ ⁢ ( Id − E ⁢ T ) = 0 . Then

and so

Finally, if g δ := P δ ⁢ g for some uniformly bounded projector P δ , then

and by combining (1.1), (1.2) and (1.3), we conclude that z δ + v δ is a quasi-best approximation from X δ to 𝑢.

An a posteriori estimator of the error in u − ( z δ + v δ ) must control the error in both z − z δ in 𝑋 and in g − g δ in Y 2 ′ . Reliable error estimators for the latter term have been developed in [14] (for Dirichlet data) and [5] (for Neumann data), but require additional smoothness of 𝑔 beyond being in Y 2 ′ (Dirichlet and Neumann data are required to be in H 1 - or L 2 -spaces on the boundary).

1.4 Outline

In Section 2, we recall the principle of a MINRES discretization and recall approaches to deal with the situation when the residual is measured in a dual norm ∥ ⋅ ∥ Y ′ and when possibly also the norm on 𝑌 cannot be evaluated. In Section 3, we present an alternative solution for the second problem by equipping 𝑌 with the optimal test norm. In addition, we discuss a posteriori error estimation. The method from Section 3 requires two uniform inf-sup conditions to be valid. In Section 4, we discuss them for four different well-posed variational formulations of Poisson’s problem with (generally inhomogeneous) Dirichlet and/or Neumann boundary conditions. Numerical results are presented in Section 5.

2.1 Well-Posed Operator Equation

For some Hilbert spaces 𝑋 and 𝑌, for convenience over ℝ, an operator G ∈ L ⁢ is ⁢ ( X , Y ′ ) and an f ∈ Y ′ , we consider the equation

With the notation G ∈ L ⁢ is ⁢ ( X , Y ′ ) , we mean that 𝐺 is a boundedly invertible linear operator X → Y ′ , i.e., G ∈ L ⁢ ( X , Y ′ ) and G − 1 ∈ L ⁢ ( Y ′ , X ) . In the following, 𝑢 will always denote the solution of (2.1).

For any closed, in applications finite-dimensional subspace { 0 } ⊊ X δ ⊊ X from a family { X δ } δ of such subspaces^[1], consider the minimal residual or least squares approximation

2.2 Discretizing the Norm on Y ′ and Saddle-Point Formulation

Unless 𝑌 is such that the Riesz map R Y : Y ′ → Y can be efficiently evaluated, i.e., 𝑌 is a (Cartesian product of) L 2 -space(s), the minimizer u δ from (2.2) is not computable.

Therefore, for a closed, in applications finite-dimensional subspace { 0 } ⊊ Y δ = Y δ ⁢ ( X δ ) ⊊ Y with

we replace above u δ by

Assuming inf δ γ δ > 0 , the following theorem shows that u δ is a quasi-best approximation to 𝑢 from X δ .

The solution u ∈ X of G ⁢ u = f is equal to argmin w ∈ X 1 2 ⁢ ∥ G ⁢ w − f ∥ Y ′ 2 , and so it solves the corresponding Euler–Lagrange equations

Introducing λ := R Y ⁢ ( f − G ⁢ u ) , the pair ( λ , u ) ∈ Y × X is the unique solution of

Completely analogously, the minimal residual approximation u δ from (2.4) can be computed as the second component of the pair ( λ δ , u δ ) ∈ Y δ × X δ that solves

2.3 Changing the Norm on Y δ

In several applications, not only ∥ ⋅ ∥ Y ′ but also ⟨ ⋅ , ⋅ ⟩ Y cannot be efficiently evaluated. This occurs for example when 𝑌 is a Cartesian product of spaces with one or more of them being fractional Sobolev spaces. Therefore, let ⟨ ⋅ , ⋅ ⟩ Y δ be an (efficiently evaluable) scalar product on Y δ , whose corresponding norm ∥ ⋅ ∥ Y δ satisfies, for some 0 < m δ ≤ M δ < ∞ ,

Now we replace u δ from (2.4) by

being the second component of the pair ( λ δ , u δ ) ∈ Y δ × X δ that solves

Assuming inf δ γ δ > 0 and sup δ M δ / m δ < ∞ , the next result shows that this u δ is a quasi-best approximation to 𝑢 from X δ .

Concerning the selection of ∥ ⋅ ∥ Y δ , let K δ = K δ ′ ∈ L ⁢ is ⁢ ( Y δ ′ , Y δ ) be an operator whose application can be computed efficiently. Such an operator is called a preconditioner for A δ ∈ L ⁢ is ⁢ ( Y δ , Y δ ′ ) defined by ( A δ ⁢ v ) ⁢ ( v ~ ) := ⟨ v , v ~ ⟩ Y . Setting

the corresponding norm ∥ ⋅ ∥ Y δ satisfies (2.6) with m δ = λ min ⁢ ( K δ ⁢ A δ ) and M δ = λ max ⁢ ( K δ ⁢ A δ ) . By substituting above choice for ⟨ ⋅ , ⋅ ⟩ Y δ in (2.7) and by subsequently eliminating λ δ from the saddle-point system, one infers that u δ ∈ X δ can be computed as the solution of the symmetric positive definite system

3.1 Optimal Test Norm

In applications, the solution proposed in Section 2.3 to circumvent the evaluation of a “difficult” scalar product ⟨ ⋅ , ⋅ ⟩ Y by means of the introduction of a preconditioner may require quite some efforts concerning coding. This holds true when 𝑌 is a Cartesian product of spaces with one or more of them being fractional Sobolev spaces on a manifold. In the current section, we give an alternative for this approach by replacing the canonical norm on 𝑌 by the so-called optimal test norm.

As a consequence of G ∈ L ⁢ is ⁢ ( X , Y ′ ) , it holds that G ′ ∈ L ⁢ is ⁢ ( Y , X ′ ) . From here on, we replace the canonical norm on 𝑌 by the equivalent norm ∥ G ′ ⋅ ∥ X ′ and, correspondingly, equip Y ′ with the resulting dual norm

We conclude that, with respect to these new norms on 𝑌 and Y ′ , G ′ : Y → X ′ and G : X → Y ′ are isometries. For this reason, ∥ G ′ ⋅ ∥ X ′ is known as the optimal test norm on 𝑌 (see [2, 23, 8, 6]).

3.2 Discretizing the Norm on Y ′ and Saddle-Point Formulation

Also with the new norm on Y ′ , the minimizer in (2.2) cannot be computed. Following Section 2.2, writing this norm in dual form sup 0 ≠ v ∈ Y | ⋅ ( v ) | ∥ G ′ ⁢ v ∥ X ′ , we discretize it by supremizing over Y δ ∖ { 0 } and rewrite the resulting least-squares problem in saddle-point form.

So we consider

being the second component of ( λ δ , u δ ) ∈ Y δ × X δ that solves

Theorem 2.1 applies, where now | | | ⋅ | | | X reads as the canonical norm ∥ ⋅ ∥ X , i.e.,

with the definition of γ δ reading as

3.3 Discretizing the Norm on X ′ and Saddling the System Once More

Unless 𝑋 is such that the Riesz map R X : X ′ → X can be efficiently evaluated, i.e., 𝑋 is a (Cartesian product of) L 2 -space(s), as in Section 2.3, we are in the situation that the norm on Y δ , here ∥ G ′ ⋅ ∥ X ′ , cannot be evaluated so that (3.1) does not correspond to an implementable method. Similar to Section 2.3, we will therefore replace this norm by an equivalent one.

Let { 0 } ⊊ X ̂ δ = X ̂ δ ⁢ ( Y δ ) ⊊ X be a finite-dimensional subspace for which^[2]

Setting, on X ̂ δ ′ ⊃ X ′ , ∥ ⋅ ∥ X ̂ δ ′ := sup 0 ≠ w ∈ X ̂ δ | ⋅ ( w ) | ∥ w ∥ X and denoting with ⟨ ⋅ , ⋅ ⟩ X ̂ δ ′ the corresponding bilinear form, we have

being the counterpart of (2.6). As in Section 2.3, we now replace (3.1) by the problem of finding ( λ δ , u δ ) ∈ Y δ × X δ that solves

Recalling that, with our current norm ∥ G ′ ⋅ ∥ X ′ on 𝑌, the norm | | | ⋅ | | | X reads as ∥ ⋅ ∥ X , an application of Proposition 2.2 shows the following.

To turn (3.5) into an equivalent evaluable system, we introduce an additional variable. With the Riesz map R X ̂ δ : X ̂ δ ′ → X ̂ δ , let θ δ := R X ̂ δ ⁢ G ′ ⁢ λ δ , i.e.,

Then

By substituting the latter relation into (3.5) and by adding the preceding equation which defined θ δ , we arrive at the equivalent problem of finding ( θ δ , λ δ , u δ ) ∈ X ̂ δ × Y δ × X δ that satisfies

for all ( θ ~ δ , λ ~ δ , u ~ δ ) ∈ X ̂ δ × Y δ × X δ , or, in equivalent block form,

Notice that, in comparison to (2.5), the “difficult” scalar product ⟨ ⋅ , ⋅ ⟩ Y completely disappeared from the system and that, apart from the usually “easy” scalar product ⟨ ⋅ , ⋅ ⟩ X , only the bilinear form ( G ⋅ ) ( ⋅ ) has to be implemented. To satisfy the inf-sup conditions γ δ > 0 and μ δ > 0 , the auxiliary spaces Y δ and X ̂ δ have to be sufficiently large (in any case, dim X ̂ δ ≥ dim Y δ ≥ dim X δ is needed), which makes solving (3.6) computationally relatively expensive. On the other hand, the implementation of the method is quite simple, whereas Proposition 3.1 shows that, for “large” Y δ and X ̂ δ , the obtained solution will be close to the best approximation from X δ . Indeed, notice that, for Y δ = Y and X ̂ δ = X , u δ is the best approximation to 𝑢 from X δ with respect to ∥ ⋅ ∥ X , whereas the second line in (3.6) shows that then θ δ = u − u δ .

3.4 A Posteriori Error Estimation

As is well known, γ δ > 0 is guaranteed by existence of a Π δ ∈ L ⁢ ( Y , Y δ ) with

where γ δ ≥ ∥ G ′ ⁢ Π δ ⁢ G ′ − 1 ∥ L ⁢ ( X ′ , X ′ ) − 1 , and conversely, γ δ > 0 guarantees existence of such a “Fortin interpolator”, being even a projector onto Y δ , with ∥ G ′ ⁢ Π δ ⁢ G ′ − 1 ∥ L ⁢ ( X ′ , X ′ ) = 1 / γ δ (e.g. [9, Lemma 26.9] or [20, Proposition 5.1]).

In the following, let Π δ be a Fortin interpolator with ∥ G ′ ⁢ Π δ ⁢ G ′ − 1 ∥ L ⁢ ( X ′ , X ′ ) ≂ 1 / γ δ (so that inf δ γ δ > 0 is equivalent to sup δ ∥ G ′ ⁢ Π δ ⁢ G ′ − 1 ∥ L ⁢ ( X ′ , X ′ ) < ∞ ).

The derivation of the a posteriori error estimator in this subsection is similar to that in [7] specialized to the use of the optimal test norm on 𝑌. Modifications were needed because of the replacement on Y δ of the optimal test norm ∥ G ′ ⋅ ∥ X ′ by ∥ G ′ ⋅ ∥ X ̂ δ ′ and the introduction of the extra variable θ δ . In [7], the 𝑌-norm of the approximate Riesz lift λ δ ∈ Y δ of the residual f − G ⁢ u δ was used as the a posteriori error estimator, whereas we use the 𝑋-norm of the approximate error θ δ ∈ X ̂ δ . Local 𝑋-norms of θ δ ∈ X ̂ δ will turn out to be effective for driving an adaptive finite element method.

4.1 Model Elliptic Second-Order Boundary Value Problem

On a bounded Lipschitz domain Ω ⊂ R d , where d ≥ 2 , and closed Γ D , Γ N ⊂ ∂ Ω , with

consider the following elliptic second-order boundary value problem:

where B ∈ L ⁢ ( H 1 ⁢ ( Ω ) , L 2 ⁢ ( Ω ) ) and A ⁢ ( ⋅ ) = A ⁢ ( ⋅ ) ⊤ ∈ L ∞ ⁢ ( Ω ) d × d satisfies ξ ⊤ ⁢ A ⁢ ( ⋅ ) ⁢ ξ ≂ ∥ ξ ∥ 2 ( ξ ∈ R d ). We assume that the matrix 𝐴 and the first-order operator 𝐵 are such that^[3]

4.2 Well-Posed Variational Formulations

In [21, 16], the following variational formulations i–iv of (4.1) have been shown to be well-posed. From the formulations given, one easily derives the expressions for 𝐺, 𝑢, 𝑓, 𝑋 and 𝑌. Implicitly we will assume that the data 𝑔, h D and h N are such that 𝑓 is in the dual of 𝑌. In the case of homogeneous essential boundary conditions, formulations i–iii can be simplified by incorporating such conditions in the definition of the domain 𝑋, an option we do not consider here.

1. (Second-order weak formulation) Find u ∈ H 1 ⁢ ( Ω ) such that

   ∫ Ω A ⁢ ∇ u ⋅ ∇ v 1 + B ⁢ u ⁢ v 1 ⁢ d ⁢ x + ∫ Γ D u ⁢ v 2 ⁢ d s = g ⁢ ( v 1 ) + ∫ Γ N h N ⁢ v 1 ⁢ d s + ∫ Γ D h D ⁢ v 2 ⁢ d s

   for all ( v 1 , v 2 ) ∈ H 0 , Γ D 1 ⁢ ( Ω ) × H ̃ − 1 2 ⁢ ( Γ D ) .

2. (First-order mild formulation) Find ( p ⃗ , u ) ∈ H ⁢ ( div ; Ω ) × H 1 ⁢ ( Ω ) such that

   ∫ Ω ( p ⃗ − A ⁢ ∇ u ) ⋅ v ⃗ 1 + ( B ⁢ u − div ⁡ p ⃗ ) ⁢ v 2 ⁢ d ⁢ x + ∫ Γ D u ⁢ v 3 ⁢ d s + ∫ Γ N p ⃗ ⋅ n ⃗ ⁢ v 4 ⁢ d s = ∫ Ω g ⁢ v 2 ⁢ d x + ∫ Γ D h D ⁢ v 3 ⁢ d s + ∫ Γ N h N ⁢ v 4 ⁢ d s

   for all ( v ⃗ 1 , v 2 , v 3 , v 4 ) ∈ L 2 ⁢ ( Ω ) d × L 2 ⁢ ( Ω ) × H ̃ − 1 2 ⁢ ( Γ D ) × H 00 1 2 ⁢ ( Γ N ) .

3. (First-order mild-weak formulation) Find ( p ⃗ , u ) ∈ L 2 ⁢ ( Ω ) d × H 1 ⁢ ( Ω ) such that

   ∫ Ω ( p ⃗ − A ⁢ ∇ u ) ⋅ v ⃗ 1 + p ⃗ ⋅ ∇ v 2 + B ⁢ u ⁢ v 2 ⁢ d ⁢ x + ∫ Γ D u ⁢ v 3 ⁢ d s = g ⁢ ( v 2 ) + ∫ Γ N h N ⁢ v 2 ⁢ d s + ∫ Γ D h D ⁢ v 3 ⁢ d s

   for all ( v ⃗ 1 , v 2 , v 3 ) ∈ L 2 ⁢ ( Ω ) d × H 0 , Γ D 1 ⁢ ( Ω ) × H ̃ − 1 2 ⁢ ( Γ D ) .

4. (First-order ultra-weak formulation) Assuming B ⁢ w = b ⃗ ⋅ ∇ w + c ⁢ w for some b ⃗ ∈ L ∞ ⁢ ( Ω ) d and c ∈ L ∞ ⁢ ( Ω ) , find ( p ⃗ , u ) ∈ L 2 ⁢ ( Ω ) d × L 2 ⁢ ( Ω ) such that

   ∫ Ω A − 1 ⁢ p ⃗ ⋅ v ⃗ 1 + u ⁢ div ⁡ v ⃗ 1 + ( b ⃗ ⋅ A − 1 ⁢ p ⃗ + c ⁢ u ) ⁢ v 2 + p ⃗ ⋅ ∇ v 2 ⁢ d ⁢ x = ∫ Γ D h D ⁢ v ⃗ 1 ⋅ n ⃗ ⁢ d s + g ⁢ ( v 2 ) + ∫ Γ N h N ⁢ v 2 ⁢ d s

   for all ( v ⃗ 1 , v 2 ) ∈ H 0 , Γ N ⁢ ( div ; Ω ) × H 0 , Γ D 1 ⁢ ( Ω ) .