# A locally calculable P^3-pressure using a P^4-velocity for incompressible Stokes equations

This paper will suggest a new finite element method to find a P^4-velocity and a P^3-pressure solving incompressible Stokes equations at low cost. The method solves first the decoupled equation for a P^4-velocity. Then, using the calculated velocity, a locally calculable P^3-pressure will be defined component-wisely. The resulting P^3-pressure is analyzed to have the optimal order of convergence. Since the pressure is calculated by local computation only, the chief time cost of the new method is on solving the decoupled equation for the P^4-velocity. Besides, the method overcomes the problem of singular vertices or corners.

## Authors

• 2 publications
04/12/2021

### Successive finite element methods for Stokes equations

This paper will suggest a new finite element method to find a P^4-veloci...
06/21/2020

### A stabilizer-free pressure-robust finite element method for the Stokes equations

In this paper, we introduce a new finite element method for solving the ...
12/06/2017

### Projection Method for Solving Stokes Flow

Various methods for numerically solving Stokes Flow, where a small Reyno...
05/28/2020

### Fully implicit and accurate treatment of jump conditions for two-phase incompressible Navier-Stokes equation

We present a numerical method for two-phase incompressible Navier-Stokes...
08/12/2020

### On Fourier analysis of polynomial multigrid for arbitrary multi-stage cycles

The Fourier analysis of the p-multigrid acceleration technique is consid...
09/13/2019

### Error Analysis of Supremizer Pressure Recovery for POD based Reduced Order Models of the time-dependent Navier-Stokes Equations

For incompressible flow models, the pressure term serves as a Lagrange m...
02/08/2022

### A stabilized formulation for the solution of the incompressible unsteady Stokes equations in the frequency domain

A stabilized finite element method is introduced for the simulation of t...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

High order finite element methods for incompressible Stokes equations have been developed well in 2 dimensional domain and analyzed along to the inf-sup condition [4, 6, 8]

. They, however, endure their large degrees of freedom and have to avoid singular vertices or corners.

Recently, we have found a so called sting function causing the discrete Stokes equation to be singular on the presence of exactly singular vertices in the Scott-Vogelius finite element method. Even on nearly singular vertices, the pressure solution is easy to be spoiled. To fix the problem, a new error analysis was devised in a successive way and restored the ruined pressure by simple post-process [7].

In this paper, employing the precedented new error analysis, we will suggest a new finite element method to find a -velocity and a -pressure solving incompressible Stokes equations at low cost.

The method will solve first the decoupled equation for a divergence-free -velocity. Then, using the calculated velocity, we will define a -pressure by exploiting locally calculable components in the Falk-Neilan and Scott-Vogelius finite element methods [4, 6, 8]. The resulting -pressure is analyzed to have the optimal order of convergence.

Since the pressure is calculated by local computation only, the chief time cost of the new method is on solving the decoupled equation for the -velocity. Besides, the method overcomes the problem arising from the singular vertices or corners.

The suggested -pressure consists of several components, each of whom needs others to be defined in a successive way. In the overall paper, the characteristics of sting functions depicted in Figure 1-(a) will also play a key role as in the precedented work in [7].

The paper is organized as follows. In the next section, the detail on finding a -velocity will be offered. Based on the orthogonality of sting and non-sting functions, we will introduce an orthogonal decomposition of the space of -pressures in Section 3. The non-sting and constant components will be defined in Section 4. Then, Section 5 will be devoted to defining 3 sting components for regular vertices, singular non-corner vertices, singular corners, respectively. It will be done on clustering sting functions vertex-wisely. Finally, after summing up all the components to define a -pressure in Section 6, numerical tests will be presented in the last section.

Throughout the paper, for a set , standard notations for Sobolev spaces are employed and is the space of all whose integrals over vanish. We will use , and for the norm, seminorm for and inner product, respectively. If , it may be omitted in the subscript. Denoting by , the space of all polynomials of degree less than or equals , will mean that coincides with a polynomial in on .

## 2 Velocity from the decoupled equation

Let be a simply connected polygonal domain in . In this paper, we will approximate a pair of velocity and pressure which satisfies an incompressible Stokes equation:

 (∇u,∇v)+(p,divv)+(q,divu)=(f,v) for all (v,q)∈[H10(Ω)]2×L20(Ω), (1)

for a body force .

Given a family of shape-regular triangulations of , define as the following space of piecewise polynomials:

 Pkh(Ω)={vh∈L2(Ω) : vh∣∣K∈Pk for all triangles K∈Th},k≥0.

Let be a space of -Argyris triangle elements [2, 3] such that

 Σh,0=P5h(Ω)∩H20(Ω), (2)

where

 H20(Ω)={ϕ∈H2(Ω) : ϕ,ϕx,ϕy∈H10(Ω)}.

The degrees of freedom of are the Hessians, gradients and values of at vertices and its normal derivatives at midpoints of edges.

Define a divergence-free space as

 Zh,0={(ϕh,y,−ϕh,x)  : ϕh∈Σh,0}⊂[P4h(Ω)∩H10(Ω)]2.

We note that

 dimΣh,0=dimZh,0=6Vin+Ein+Vbdy−Vcnr,

where , and are the numbers of interior vertices, interior edges, boundary vertices and corners, respectively [4].

Then, we can solve satisfying the following decoupled equation:

 (∇uh,∇vh)=(f,vh) for all vh∈Zh,0. (3)
###### Theorem 2.1.

Let and satisfy (1), (3

), respectively. Then we estimate

 |u−uh|1≤Ch4|u|5, (4)

if , where is a constant independent of .

###### Proof.

Since satisfies (1), we have . Thus, there exists a stream function such that [5]

 u=(ϕy,−ϕx). (5)

Let be a projection of such that the Hessians, gradients and values of at vertices and its normal derivatives at midpoints of edges vanish. Then, since , by Bramble-Hilbert lemma, we have

 |ϕ−Πhϕ|2≤Ch4|ϕ|6. (6)

If we denote , then from (5) and (6), we estimate

 |u−Πhu|1≤Ch4|u|5. (7)

We note that for all . Thus, from (1) and (3), we deduce

 (∇u−∇uh,∇vh)=0 for % all vh∈Zh,0.

It is written in the form:

 (∇Πhu−∇uh,∇vh)=(∇Πhu−∇u,∇vh) for all vh∈Zh,0. (8)

Then we can establish (4) from (7) and (8) with . ∎

## 3 Orthogonal decomposition of P3

For a triangle and an integer , define

 Pk(K)={q∈L2(Ω): q∣∣K∈Pk, q=0 on Ω∖K}.

In the remaining of the paper, we will use the following notations:

• : a generic constant which does not depend on of ,

• : the union of all triangles in sharing a vertex as in Figure 4,

• : the counterclockwise

rotation of a vector

,

• : the area or length of a set ,

• : the average of a function over .

We assume the following on to exclude pathological meshes.

###### Assumption 3.1.

Every triangles in has at most one corner point of .

### 3.1 sting function

Let be a vertex of a triangle . Then there exists a unique function satisfying the following quadrature rule:

 ∫KsVK q dxdy=|K|100q(V) for all q∈P3, (9)

since the both sides of (9) are linear functionals on . If is a reference triangle with vertices and , we have

 sˆVˆK(x,y)=285y3−6310y2+95y−110, (10)

as depicted in Figure 1-(a). Given a vertex of , we note that

 sVK=sˆVˆK∘F−1 for an affine transformation F:ˆK⟶K. (11)

Thus, the values of are inherited from those of as

 (12)

If for a scalar , we call it a sting function of on , named after the shape of its graph as in Figure 1-(a).

For a triangle , define a subspace of as

 S(K)=, (13)

where are 3 vertices of . From (12), it is easy to prove that

 dimS(K)=3. (14)

### 3.2 non-sting function

For a triangle , let

 B(K)={v∈[P4(K)]2 : v=0 on % ∂K}, (15)

and define a subspace of as

 N(K)={divv : v∈B(K)}. (16)

If , we will call it a non-sting function on . By definition in (15), (16), every non-sting function has the following properties:

 q(V)=0 for every vertex V of K,∫Kq dxdy=0. (17)

An example of its graph is depicted in Figure 1-(b).

Then, the following orthogonality is clear from (13), (17) and the quadrature rule in (9),

 N(K)⊥S(K). (18)

The fact induces the following, with an aid of Lemma 3.1 below,

 dimN(K)=6. (19)

If and , then .

###### Proof.

Since on and on , there exists such that

 (ϕy,−ϕx)=v on K,ϕ=0 on ∂K.

Let be 3 infinite lines containing 3 line segments of , respectively. Then vanish on . It implies that vanishes on any line which passes 3 points in . Thus we have and on . ∎

The above lemma tells that is a norm of . Furthermore, we can show that

 |v|1,K≤C∥divv∥0,K for all v∈B(K). (20)

Actually, in (16) is the space of all function satisfying (17) [6].

### 3.3 orthogonal decomposition

Let be a constant function of value 1 on . Then, we can decompose as in the following lemma. We will notate

 A⊥⨁B for A⨁B, if A⊥B.
###### Lemma 3.2.
 P3(K)=N(K)⊥⨁(<1K>⨁S(K)). (21)
###### Proof.

From (17), we have . Thus, by (14), (18), (19), it is enough to prove (21) that

 1K∉S(K).

For the vertices of , let

 sK1=54(sV1K+sV2K+sV3K)∈S(K). (22)

Then, by (12), vanishes at all vertices of . By (9), it means that

 (23)

Assume . Then, and from (23). It contradicts to

 ∫KsK1 dxdy=54∫KsV1K+sV2K+sV3K dxdy=380|K|≠|K|=∫K1K dxdy.

Let’s define the following subspaces of :

 Nh=⨁K∈ThN(K),Ch=⨁K∈Th<1K>,Sh=⨁K∈ThS(K). (24)

Then by Lemma 3.2, we have

 P3h(Ω)=Nh⊥⨁(Ch⨁Sh). (25)

### 3.4 decomposition of Πhp

For satisfying (1), let

be a Hermite interpolation of

such that

 ∇Πhp(V)=∇p(V),Πhp(V)=p(V),Πhp(G)=p(G), (26)

at all vertices and gravity centers of triangles in . Then, if , we have

 ∥p−Πhp∥0≤Ch4|p|4. (27)

By (25), we can decompose into

 Πhp=ΠhpN+ΠhpC+ΠhpS, (28)

for

 ΠhpN∈Nh,ΠhpC∈Ch,ΠhpS∈Sh,

called the non-sting, constant and sting components of , respectively. We will approximate them component-wisely in next two sections, exploiting the following equation for ,

 (Πhp,divv)=(f,v)−(∇u,∇v)−(p−Πhp,divv) for all v∈[H10(Ω)]2. (29)

## 4 Non-sting and constant components

Define the following spaces:

 Vh,0={vh∈[P4h(Ω)∩H10(Ω)]2 : ∇vh is continuous at all vertices in Th},Qh,0={qh∈P3h(Ω)∩L20(Ω) : qh is continuous at all vertices in Th},Qh,00={qh∈Qh,0 : qh % vanishes at all corners of ∂Ω}. (30)

Then, for in (1), there exists a unique satisfying

 (31)

since the following Stokes complex is exact with in (2) [4]:

 0⟶Σh,0curl−−→Vh,0div−−→Qh,00⟶0.

We note that in (31) is the previously calculated in (3) and is the unique solution of the following equation,

 (rh,divvh)=(f,vh)−(∇uh,∇vh) for all vh∈Vh,0. (32)

Similarly to (28), decompose into

 rh=rNh+rCh+rSh % for rNh∈Nh, rCh∈Ch, rSh∈Sh. (33)

### 4.1 non-sting component pNh

For a triangle , we note . Thus, by (24), (33) and orthogonality in Lemma 3.2, we can rewrite (32) into

 (rNh,divvh)=(f,vh)−(∇uh,∇vh) for all% vh∈B(K). (34)

Then by definition of in (16), is determined by (34). Thus we can calculate locally on each triangle in .

###### Lemma 4.1.

Define in (33), then for in (28), we estimate

 (35)
###### Proof.

From (29) and Lemma 3.2, satisfies that

 (36)

If we set , from (34) and (36), we have

 (37)

Then (35) comes from (20), (37) and the definition of in (16). ∎

### 4.2 piecewise constant component pCh

Define a subspace of as

 Vh,00={vh∈Vh,0 : ∇vh vanishes at all vertices in Th}. (38)

Then, by quadrature rule in (9) and definition of in (38), we can reduce (32) into

 (rCh,divvh)=(f,vh)−(∇uh,∇vh)−(pNh,divvh) for all vh∈Vh,00, (39)

where defined in Lemma 4.1.

For each triangle , let be a constant such that

 CK=rCh∣∣K. (40)

If are two adjacent triangles sharing an edge, there exists a test function such that

 the support of vh is in K1∪K2,∫Kjdivvh dxdy=(−1)j+1,j=1,2. (41)

Then from (39)-(41), we have

 CK1−CK2=(f,vh)−(∇uh,∇vh)−(pNh,divvh). (42)

Fix a triangle . Then for each triangle , we can choose a telescoping sequence of triangles:

 K0,K1,K2,⋯,Kn=K, (43)

so that share an edge for . Then, can be obtained from

 CK0−CK1,CK1−CK2,⋯,CKn−1−CKn,

which are calculated locally as in (42).

The knowledge of for all gives us from

 0=∫ΩrCh−\mathpzcm(rCh) dxdy=∑K∈ThCK|K|−\mathpzcm(rCh)|Ω|=∑K∈Th(CK−CK0)|K|+(CK0−\mathpzcm(rCh))|Ω|.

Now, we can obtain , since for each ,

 rCh∣∣K−\mathpzcm(rCh)=CK−CK0+CK0−\mathpzcm(rCh).
###### Remark 4.2.

The difference does not depend on the choice of telescoping sequence in (43) from the existence of satisfying (39).

###### Lemma 4.3.

Define for in (33), then for in (28), we estimate

 ∥ΠhpC−\mathpzcm(ΠhpC)−pCh∥0≤C(|u−uh|1+∥p−Πhp∥0). (44)
###### Proof.

By (29) and quadrature rule in (9), satisfies for all ,

 (45)

If we set , from (39) and (45), it satisfies

 (46)

for all .

We can establish Lemma 4.4 below that tells the existence of a nontrivial such that

 β∥eh∥0|vh|1≤(eh,divvh) for β>0 regardless of h. (47)

Then (44) comes from (46), (47) and Lemma 4.1. ∎

###### Lemma 4.4.

satisfies the inf-sup condition.

###### Proof.

Given , there exists a nontrivial such that [1]

 β∥ch∥0|w|1≤(ch,divw), (48)

for a constant regardless of .

For each triangle in , define so that

 ∇zK=∇w at all % vertices of K,∫EzK dℓ=0 for each edge E of K,zK=0 at all 3 vertices and midpoints of 3 medians of % K. (49)

Then, for a reference triangle and an affine map , we have

 |zK|1,K≤C|zK∘F|1,^K≤C|w∘F|1,^K≤C|w|1,K. (50)

If we define by for all , then belongs to , since derivatives of along to edges are continuous. We note that . If so, from the second and third conditions in (49).

Thus, from (49) and (50), satisfies

 w−z∈Vh,00∖{0},|z|1≤C|w|1,(1,divz)K=0 for all K. (51)

Then, the following comes from (48) and (51), which completes the proof:

For each triangle , define as in (22) so that vanishes at all vertices of . If we modify in (33) by

 rCh′=∑K∈ThrCh∣∣K(1K−sK1),rSh′=rSh+∑K∈ThrCh∣∣KsK1,

we have

 rh=rNh+rCh′+rSh′. (52)

We note is continuous at every vertex, since so is and vanish at all vertices in . By (32) and (52), the vertex-continuous sting function satisfies the following system:

 (rSh′,divvh)=(f,vh)−(∇uh,∇vh)−(rNh+rCh′,divvh) for all vh∈Vh,0. (53)

If has no singular corner, we could revise (53) slightly with the Falk-Neilan finite element space [4], to obtain a sting component in approximating of in (28).

To avoid a global system such as (53), we will seek locally calculable sting components vertex-wisely in the next section by unleashing the vertex-continuity.

## 5 Sting component

For each vertex , let be the space of all sting functions of , that is,

 S(V)=,

where are all triangles in sharing . An example of a function in is represented in Figure 2.

Then we note that

 Sh=⨁K∈ThS(K)=⨁V:vertexS(V).

Thus the sting component of in (28) can be decomposed by vertex into

 ΠhpS=∑V:vertexΠhpV for ΠhpV∈S(V). (54)

In this section, we will define to approximate for regular vertex , singular non-corner vertex , singular corner , respectively in order: .

### 5.1 two test functions for two adjacent triangles

Let be two adjacent triangles sharing an edge and a vertex as in Figure 3. Denote other 3 vertices and a unit tangent vector by so that

 ¯¯¯¯¯¯¯¯¯¯¯¯¯VW0=K1∩K2,τ=−−−−→VW0|¯¯¯¯¯¯¯¯¯¯¯¯¯VW0|,Wj∈Kj∖{V,W0},j=1,2.

Then, there exists a function such that [7]

 ∂w∂τ(V)=1, ∂w∂τ(W0)=0, ∫¯¯¯¯¯¯¯¯¯¯¯¯VW0w dℓ=0,  the support of w is K1∪K2. (55)

Assuming are counterclockwisely numbered with respect to , by simple calculation, we have

 ∇w∣∣K1(V)=|¯¯¯¯¯¯¯¯¯¯¯¯¯VW0|2|K1| −−−−→VW1⊥,∇w∣∣K2(V)=−|¯¯¯¯¯¯¯¯¯¯¯¯¯VW0|2|K2| −−−−→VW2⊥. (56)

For a vector , denote . Then from (55) and (56), vanishes at all vertices in except

 divwξh∣∣K1(V)=|¯¯¯¯¯¯¯¯¯¯¯¯¯VW0|2|K1| −−−−→VW1⊥⋅ξ,divwξh∣∣K2(V)=−|¯¯¯¯¯¯¯¯¯¯¯¯¯VW0|2|K2| −−−−→VW2⊥⋅ξ. (57)

Let be a sting function on . It is represented with some constants as

Then, by (57) and quadrature rule of sting functions in (9), we have that

It can be written in simpler form:

 (qh,divwτh)=ℓ1ℓ0sinθ1200α1+ℓ0ℓ2sinθ2200β1=1100(|K1|α1+|K2|β1),(qh,divwτ⊥h)=ℓ1ℓ0cosθ1200α1−ℓ0ℓ2cosθ2200β1, (58)

where