CORE Metadata, citation and similar papers at core.ac.

uk
Provided by Brunel University Research Archive

TR/24 February 1974

LOCAL MESH REFINEMENT WITH

FINITE ELEMENTS FOR ELLIPTIC PROBLEMS

by

J.A. Gregory and J.R. Whiteman

 1.

1. Introduction

When a finite element technique is used for the solution of

a two dimensional elliptic boundary value problem, P, the region

of definition Ω of the problem is divided in a number of non-overlapping

elements. In this paper the region Ω is rectangular, the elements are

all rectangles and we consider a technique of local mesh refinement.

A weak formulation of the problem P is constructed and it is the

solution of this weak problem, the generalized solution of P lying in

a larger space, W, rather than the strong solution of P, which is

sought. In the Galerkin technique an approximation U(x,y) to the

generalized solution u (x,y) is constructed from a finite-dimensional

space Sh, which is usually a subspace of W, where h is a parameter

designating the size of the elements. The key step in the successful

application of the Galerkin method is the construction of Sh which

generally consists of functions which are piecewise polynomial over Ω.

In each element the approximating function is derived from an interpola-

tion function which interpolates the values of u, and frequently also

certain derivatives of u, at nodes in the element.

Let the interpolant in each element have the form

~ k
u (x, y ) = ∑ u i φ i ( x, y ) , (1.1)
i=1

where the Φi. are the cardinal basis functions of the interpolation

with respect to ui , the function and derivative values at the

nodal points. The approximating function then has in each element

the form

k
U ( x, y ) = ∑ U i φi ( x, y ) , . (1.2)
i=1

where now the Ui. are the corresponding values of U and certain

of its derivatives at the nodal points of the element. The unknown

 2.

values of Ui are determined by the Galerkin method. The functions

φ i considered over the totality of elements generate the finite

dimensional space Sh .

The global piecewise polynomial approximating function must

satisfy certain continuity requirements across interelement

boundaries in order that Sh be a subspace of W. This is the

conforming condition. For Poisson type problems the conforming

condition is that Sh ⊂ Co ( Ω ), whilst for biharmonic problems it is

that Sh ⊂ C1 ( Ω ); see e.g. Zlamal [10].

When a standard rectangular mesh is taken together with bilinear

interpolation to the function values at the four corners of each

element, the global approximating function is in Co and the conforming

condition for Poisson problems is satisfied. In this case the Lagrange

basis functions in each element are the linear pyramid functions of

the finite element method. However, if the mesh is refined locally

about some point 0, as for example in Figure 1, mid-side nodes are

introduced and special interpolants must be derived for use in these

five-node rectangle elements so that the global approximating function

may again be in C o . Such interpolants are derived in Section 2.

Figure 1
 3.

For biharmonic problems, when a standard rectangular mesh

is used, together with a bicubic interpolant in each element to

∂u ∂u ∂ 2u
the values of u, , and at each of the four corners,
∂x ∂y ∂x∂y

the global approximating function will be in C1 . Again, if the mesh

is refined locally, special interpolants are necessary in the

five-node elements for the global approximating function to he in C1 .

These are derived in Section 3.

In Section 4 a Galerkin procedure is described briefly for

a problem involving Laplace's equation, and in Section 5 numerical

results are given.

2. Co approximating functions

We consider the unit square with vertices at (0,0), (1,0),

(1,1) and (0,1). The linear interpolant to the values U ( 0 , 0 ) and

U(1,0) along [0,1 ] can be written as

U(x,0) = (1-x) U(0,0) + x U(1,0) . (2.1)

By using tensor products we see immediately that the bilinear

interpolant to U ( 0 , 0 ) , U ( 1 , 0 ) , U(1,1) and U ( 0 , 1 ) over the square is

U(x,y) = (1-x)(l-y)U(0,0)+x(l-y)U(l,0)+xyU(l,l)+(l-x)yU(0,1),

4
= ∑ Ui φ i (x, y) . (2.2)
i=1

The φ i. are the basis functions referred to in Section 1, and

use of (2.2) in each element of a regular mesh produces a C°

approximating function.

When the mesh is refined locally by successive halving

of the mesh length, as in Figure 1, mid-side nodes are introduced

 4.

and situations as in Figure 2 arise. An obvious approach for

Figure 2

1
dealing with the mid-side node ( , 0 ) is to take ( 2 . 2 ) in element 1,
2
1
which thus gives the value at ( ,0) as the linear interpolant "between
2
(0,0) and (1,0). This can then be used directly when (2.2), suitably
scaled;, is applied in elements 2 and 3. Thus the unknown value at
1
( ,0) is not introduced in element 1. However, the effect of this scheme
2
is to spread the domain of influence of the coarse mesh into the region
of fine mesh. The effect of the refinement is therefore reduced;
see e.g. Wait and Mitchell [ 5 ] where this procedure is adopted. In order
to avoid this we choose the alternative scheme given below.
A suitably scaled form of (2.2) is used in each of the elements
1
2 and 3. This function in element 2 interpolates U(0,0) and U( ,0),
2
1
whilst at the same time it is linear on L12 = {(x,y); 0 ≦ x ≦ , y=0}.
2
1
Similarly in 3 the function interpolates U ( ,0) and U ( 1 , 0 ) and is linear
2
1 1
on L13 ≡ {(x,y); ≦ x ≦ 1, y =0} . As there is a node at ( ,0) in element 1,
2 2
 5.
the interpolant (2.2) has to be adapted so that in 1 it will
1
interpolate the value U ( ,0) as well as the values at the four
2
corners of the element, whilst being linear on L12 and L 13 . The
piecewise polynomial function over the union of the elements 1,2 and 3
will then again be in Co .

The technique is to consider separately the two rectangles

1 (2.3)
R1 ≡ {(x,y) ; 0 ≦ x ≦ , 0 ≦ y ≦ 1}
2

and

1 (2.4)
R2 = {(x,y); ≦ x ≦ 1, 0≦y ≦1} .
2

In R1 the interpolating function is

1 1
U(x,y) = ( 1 - 2 x ) ( 1 - y ) U ( 0 , 0 ) + 2x (1-y)U( ,0) + 2xy U( , 1) (2.5)
2 2

+ (1-2x)yU(0,l).
However, the point ( 2 , 1 ) is not a node of the element 1, and
1
so the value U ( ,1) is eliminated using the continuity of the
2

approximating function across { ( x , y ) ; 0 ≦ x ≦ 1, y = 1 } by the

substitution of

1 1
U( , l) = (U(0,1)+ U(1,1)).
2 2

Thus for (x,y) ε R1

1
U(x,y) = (1-2x)(1-y)U(0,0) + 2x(1-y)U( ,0)+ xy U( 1,1)+ y(1-x)U(0,1 ). (2.6)
2

A similar technique is adopted in R2 so that in element 1

U(x,y) = y(1-x)U(0, 1) + xy U (1,1)

⎧ 1 1
⎪ ( 1 − 2x) (1 − y)U(0,0) + 2x( 1 − y) U)( 2 , 0 ), , 0 ≤ x < 2 ,
+ ⎨ (2.7)
1 1
⎪ 2(1 − x) (1 − y)U( , 0) + (2x − 1) (1 − y) U (1,0) , < x < 1.
⎩ 2 2
 6.
The function U(x,y) in ( 2 . 7 ) is bilinear in R1 and R2 and
1
continuous across x= , 0 ≦ y ≦ 1 . Use of trial functions of the
2
type (2.7) in the five node elements together with bilinear trial
functions in the standard elements will ensure that the resulting
global approximating function is in Co .

Note that, in terms of "+ functions" commonly used in splines,

it is possible to write

U(x,y) = y ( 1 - x ) U ( 0 , l ) + x y U ( 1 , 1 )
1
+ ( 1 - 2 x ) ( 1 - y ) U ( 0 , 0 ) + 2x(l-y)U( ,0)
2
1
+ (1-y) (2x-1)+ [U(0 ,0 )- 2U( ,0) + U ( 1 , 0 ) ] , (2.8)
2

where

⎧ 1
⎪ ( 2x − 1 ) , x>
2
,
(2x − 1) + = ⎨ 1
⎪ 0 , x< .
⎩ 2
3. C1 Approximating Functions

The notation

∂U(x, y) ∂U(x, y) ∂ 2 U(x, y)

U1,0 (x, y) ≡ , U 0,1 (x, y) ≡ , U1 , 1 (x, y) =
∂x ∂y ∂x∂y

is adopted, and with this notation the cubic interpolant to

U(0,0), U(1,0), U1,0(0,0), U1 , 0(0,1) over [0,l] can be written

U(x,0) = φ 1(x)U(0,0) + φ 2(X)U1,0(0,0) + φ 3(X)U(1,0)+ φ 4 (x)U1,0(1,0) , (3.1)

φ1 (t) = ( t − 1) 2 (2t + 1), ⎫

⎪
φ 2 (t) = ( t − 1) 2 t , ⎪⎪
where ⎬ (3.2)
φ 3 (t) ≡ φ1 (1 − t) = t 2 (−2t + 3) ,⎪
⎪
φ 4 (t) ≡ − φ 2 (1 − t) = t 2 (t − 1) . ⎪⎭
 7.

The φ 's are the cardinal basis functions for Hermite

interpolation and it is noted that, if φ '(t) = d φ (t)/dt,

φ 1 (0) = φ 3(1). = φ 2(0) = φ 4’ (1) = 1 ,

φ 2 (0) = φ 3(0) = φ 4(0) = 0 ,

φ 1 (0) = φ 3(0) = φ 4(0) = 0 ,

φ ’1(0) = φ ’3(0) = φ '4(0) = 0 ,

φ ’1(1) = φ ’2(1) = φ ’3(1) = 0 .

Taking tensor products we obtain the bicubic interpolant to

Z ≡ {U(xi,yi), U1,0 (xi, yi ),U0,1(xi,y.), U1,1(xi,yi)} , (3.3)

at respectively each of the four points (xi yj) =( 0 , 0 ) ,(1,0),(1,1),(0,1)

over the unit square as

U(x,y) = φ 1 (x) [ φ 1(y)U(0,0) + φ 2(y)U0.1(0,0) + φ 3(y)U(0,l)+ φ 4(y)U0,1,(0,1)

+ φ 2(x)[ φ 1(y)U1,0(0)+ φ 2(y)U1,1(0,0)+ φ 3(y)U1,0(0,1)+ φ 4(y)U1,1(0,1)]

+ φ 3(x) [ φ 1(y)U(1,0) + φ 2(y)U0,1 (1,0)+ φ 3 (y)U(l,l)+4(y)U0,1 (1,1)]

+ φ 4(x)[ φ 1(y)U 1, 0(1,0)+ φ 2(y)U1,1(1,0) + φ 3(y)U1,0(1,1) + φ 4(y)U1,1(1,1)],

(3.4)

where the φ 's are as in ( 3 . 2 ) . Use of (3.4) as the trial function in.

each element of a standard rectangular mesh, together with the specifying

of Z as in (3.3) at each node, will produce a C1 approximating function.

However, we wish to refine the mesh as in Figure 1, whilst retaining C1

continuity in the global approximating function. Referring again to the

situation as in Figure 2, a special trial function is thus needed in elements

such as 1. Following Section 2 we split the element 1 into the two rectangles

R1 and R2 of (2.3) and (2.4). The bicubic interpolant to the vales of Z at the
 6.
four vertices of R1 is, from (3.4),

u(x.y) = φ 1(2x)[ φ 1(y)U(0,0) + φ 2(y)U0,1 (0,0) + φ 3 (y)U(0,1) + φ 4(y)U0,1 (0,1)]

1
+ φ 2(2x)[ φ 1(y)U1, 0(0,0) + φ 2(y)U1,1 (0,0) + φ 3(y)U1,0(0,1) + φ 4(y)U1,1 (0.1)]
2
1 1 1 1
+ φ 3(2x)[ φ 1(y)U( ,0) + φ 2(y)U0,1( ,0) + φ 3(y)U( ,1) + φ 4(y)U0,1 ( ,1)]
2 2 2 2
1 1 1 1 1
+ φ 4(2x)[ φ 1(y)U1,0( ,0) + φ 2(y)U1.1 ( ,0)+ φ 3(y)U1.0 ( ,1) + φ 4(y)U1,1 ( ,1)].
2 2 2 2 2

(3.5)
1
The right hand side of (3.5) involves values of Z at the point ( ,1),
2
which is not a node of the discretization. In order that these values
may be eliminated interpolants having the form (3.1) are used on

{(x,y) ;0 ≦ x ≦ 1, y = 1 }, so that
1 1 1 1 1
U( ,1) = φ 1( )U(0,1) + φ 2( )U1,0 (0,1) + φ 3( )U(1,1) + φ 4( )U1, 0(1,1), (3.6)
2 2 2 2 2
1 1 1 1 1
U0,1 ( ,1)= φ 1 ( ) U0,1(0,1) + φ 2( )U1,1 (0,1) + φ 3( )U0,1 (1, 1) + φ 4 ( )U1,1(1,1), (3.7)
2 2 2 2 2
1 1 1 1 1
U1,0( ,1)= φ 1( )U(0,1) + φ 2( )U1,0(0,1) + φ 3( )U(1,1) + φ 4( )U1,0 (1,1), (3.8)
2 2 2 2 2
1 1 1 1 1 (3.9)
U1,1( ,1)= φ 1( )U0,1(0,1) + φ 2( )U 1,1 (0,1) + φ 3( )U0,1 (1,1) + φ 4( )U1,1 (1,1).
2 2 2 2 2

But

1 1 1 1 1 1
φ 1( ) = φ 3( ) = , φ 2( ) = - φ 4( ) - ,
2 2 2 2 2 8

1 3 1 1 1 1 3
φ ’ 1( ) = - . φ ’ 2 ( ) = φ 4( ) = - , φ 3( ) = .
2 2 2 2 4 2 2
Substitution of these values into (3 .6 ) - ( 3 . 9 ) and the subsequent
1
substitution of the resulting expressions for Z at ( ,1) in (3. 5)
2
 9.

leads to
1 1 1 1
U(x,y) = φ 1(y)[Φ1(2x)U(0,0) + φ 2(2x)U1,0(0,0) + φ 3(2x)U( ,0) + φ 4(x)U1,0 ( , 0)]
2 2 2 2
1 1
+ φ 2 (y)[ φ 1 (2x)U l, 0 (0,0)+ 2 (2x)U 1,1 (0.0)+ φ 3 (2x)U 1,0 ( ,0)
2 2
1 1
+ φ 4 (2x) Ul,1 ( ,0)]
2 2
1 3
+ φ 3(y)[{ φ 1(2x) + φ 3 + 3(2x)- φ 4(2x)} U(0,1)
2 4
1 1 1
+ { φ 2(2x) + φ 3(2x) - φ 4 (2x)}U1 , 0(0,1)]
2 8 8
1 3 1
+ { φ 3(2x)+ φ 4(2x)}U(1,1)- { φ 3(2x)+ φ 4(2x)}U1, 0(1,1)]
2 4 8
1 3
+ φ 4(y)[{ φ 1(2x)+ φ 3(2x)- φ 4(2x)}U0,1 (0,1)
2 4
1 1 1
+ { φ 2(2x) + φ 3(2x) + φ 4(2x)}U 1,1 (0,1)
2 8 8
1 3 1
+ { φ 3(2x)+ φ 4(2x)}U0 , 1(l,1)- { φ 3(2x)+ φ 4(2x)}U1,1 (1,1)]
2 4 8

(3.10)

1
for 0 ≦ x ≦ , 0 ≦ y ≦ 1.
2
An expression of a similar form to that in (3.10) is
1
obtained for the interpolant in R 2 ≡ {(x.y); ≦ x ≦ 1,0 ≦ y ≦1}.
2
This taken together with the interpolant ( 3 - 1 0 ) produces in element
1 an interpolant which is C1 and which is cubic on the sides, the top
and each of the halves of the bottom of the element. Incorporation of
this into the space of piecewise bicubic functions will ensure that the
resulting global approximating function is in C 1 .
 10.

4. Galerkin Method For Model Problem

We consider the problem in which u(x,y) satisfies

⎫
⎪
− Δ [u(x, y) ] = 0 , (x, y) ε Ω , ⎪
⎪
u(x, y) = 500 , (x, y) ε BC, ⎪
⎪⎪
u(x, y) = 0 , (x, y) ε EO , ⎬ (4.1)
∂u(x, y) ⎪
= 0 , (x, y) ε OB U CD ,⎪
∂y ⎪
⎪
∂u(x, y) ⎪
= 0 , (x, y) ε DE , ⎪⎭
∂x

where Ω is the rectangular region OBCDEO of Figure 3, in which EO=OB=BC=0.5.

The problem(4.1) is derived using symmetry from a well known problem

in a rectangle containing a slit which has been much studied; see for

example Whiteman [ 6 ], [7 ] and Wait and Mitchell [ 5 ] . We define the

two disjoint parts of the boundary ∂Ω1 ≡ BC U EO , ∂Ω 2 ≡ OB U CD U DE

and let ∂Ω = ∂Ω1 U ∂Ω2 with Ω ≡ Ω U ∂Ω.

Let W12 (Ω) be the Sobolev space of functions which together with

their generalized derivatives of order one are in L2(Ω). The subspace

of functions in W12 (Ω) which satisfy a homogeneous boundary condition

on ∂Ω1 is written W12(Ω)∩ (∂Ω1)0 ; that is for v e W12(Ω)∩ (∂Ω1)0 ,

 11.
V ε W12 (Ω) and v = 0 on ∂Ω1 .

The weak problem corresponding to (4.1) is ;

find u ε f + W12 (Ω)such that

a(u,v) =0 ∀ v ε W12 (Ω)∩ (∂Ω1)0 , (4.2)

where f ε W12( Ω ) with f=500 on BC and f=0 on EO

The notation u ε f + W12 (Ω)means that u = f+v, where v ε W12(Ω) ∩ (∂Ω1 ) 0 .

In (4.2) the bilinear functional a (u,v) is defined as

⎛ ∂u ∂v ∂u ∂v ⎞
a (u, v) = ∫∫ ⎜⎜ + ⎟⎟ dx dy ∀ u, v ε W1
2 (Ω) .
Ω ⎝ ∂x ∂x ∂y ∂y ⎠

The energy norm | | v | | E is defined by

1
||v||E = ( a ( v . v ) 2 . (4.3)

The region Ω is discretized into rectangular elements of

generic length h so that there are m internal nodes, n nodes

on ∂Ω1 and p nodes on ∂Ω2 . The global approximating function

U ε Sh is written as
m+ p n
U(x, y) = h
∑ U i B i (x, y) + ∑ f j C j (x, y) , (4.4)
i=1 j=1

where the fhj. are the known values of u =U at the nodes of

∂Ω1 and the Bi and Cj are formed from the Φi of (1,2). The Bi

and Cj. are linear pyramid basis functions at respectively the

nodal points of Ω U ∂Ω2 and the nodal points of ∂Ω1

In the Galerkin procedure we seek U as in (4.4) such that

a(U , Bk) = 0 , k= 1, 2, ... , m+p. (4. 5)

 12.
Equations (4.5) are the normal equations for the solution of

the unknown coefficients U i . Bounds on the error in the Galerkin

solution U, having the form

||u-U||E ≤ Kh|u|2 , (4.6)

where

⎧ 2 2 2 ⎫
⎪ ∂ 2u ∂ 2u ∂ 2u ⎪
| u |2 = ⎨ + + ⎬ ,
2 ∂y 2 L (Ω ) ⎪
⎪ ∂x L (Ω ) ∂x∂y L (Ω )
⎩ 2 2 2 ⎭

are well known; see e.g. Ciarlet and Raviart [ 3 ] . For problems

containing boundary singularities, for example due to re-entrant

corners such as in the problem in the slit rectangle which is

equivalent to (4.1), the second derivatives of u are not in L2.

However , error bounds involving hπ/a for problems containing

re-entrant corners with interior angles a > π have been derived;

see Babuska and Aziz [l,p.274].

The effect of the singularity is to reduce the accuracy of the

Galerkin solution, particularly in the neighbourhood of the re-entrant

corner. It is shown by Babuska [2 ] that with "proper" refinement

of the elements around the corners the effect of the singularity can

be removed. For ease and automation of computation we now use the local

refinement scheme of Section 2, which is based on successive halving

of the mesh length to obtain accurate Galerkin solutions to the

problem (4. 1).

5. Numerical Results

In applying the Galerkin technique of Section 4 to the

problem (4.1) we use square elements and take the basis functions
 13.

B i (x,y) and C j (x,y) of ( 4 . 4 ) to be pyramids in the appropriate

elements. Thus in each square of a standard mesh of length h the

local trial function has the form (2.2). For the value h= 1/14 the results

obtained at three specific points are given in Figure 4. For comparison accurate

results obtained using a conformal transformation method (CTM) due

to Whiteman and Papamichael [ 9 ] are also included in Figure 4.

It is seen that, as expected, the greatest inaccuracies occur in

the neighbourhood of 0.

The refinement scheme of Figure 1 is now used with the trial

functions (2.4) in five node elements. We note that each level of

refinement introduces 8 new equations into the system (4,5). The

ordering of the nodes is chosen so that the inclusion of the extra

equations is performed automatically. This is achieved by ordering

peripherally about the singularity. In each case the matrix of coefficients

is symmetric and positive definite, and full advantage is taken of the band

structure. The three sample points P ≡ (0,1/14) and Q ≡( 1/14,0) near to 0

and R ≡ (-3/7,3/7) remote from 0, with the origin of co-ordinates at 0, are again

chosen and results at these points are given in Figure 4 for levels of

refinement ranging from 1 to 8 with the original mesh length h again

1
taken as /14 .

It is seen that with continued local mesh refinement the stage

has been reached where the Galerkin solutions are more accurate near

the singularity than they are at a point in Ω remote from 0. The effect

of the singularity on the numerical solution has thus been neutralized

by the refinement. The error at points on the coarse mesh is due to

the coarse mesh spacing.

An alternative technique to local mesh refinement is to include

in the space Sh singular functions having the form of the singularity,

as has been done by Fix [4] and Whiteman [8] . In order to do this one
 14.

has of course to know the form of the singularity. Further

the augmentation of Sh poses difficult computational problems.

We feel that local mesh refinement is a viable alternative.

Number of Value at U (x, y) at Number of

Levels of Linear
Refinement P≡(0,
1
/14) Q ≡(
1
/14,0) R≡(3/7,3/7) Equations

0 (h=1/14) 97.05 147.05 88.73 104

1 99.61 150.52 89.78 112

2 101.62 153.39 90.31 120

3 102.72 154.92 90.57 128

4 103.27 155.69 90.70 136

5 103.54 156.07 90.78 144

6 103.68 156.26 90.80 152

7 103.75 156.36 90.82 160

8 103.78 156.40 90.83 168

CTM [9]
103.77 156.48 91. 34 -
Results

Figure 4
 15.

References.

1. Babuska I. and Aziz A. K. , Foundations of the finite element

method, pp.5-359 of A.K.Aziz (ed.), The Mathematical
Foundations of the Finite Element Method with Applications
to Partial Differential Equations. Academic Press, New York, 1972.
2. Babuska I. , Finite element method for domains with corners.
Computing 6, 264-273, 1970.
3. Ciarlet P.G and Raviart P.A. , General Lagrange and Hermite
interpolation in one and two variables with applications
to partial differential equations. Numer.Math. 11, 232-256, 1968.
4. Fix G., Higher order Rayleigh-Ritz approximations.
J.Math.Mech. 18, 645-657, 1969.
5. Wait R. and Mitchell A.R. , Corner singularities in elliptic problems
by finite element methods. J.Comp. Phys.8, 45-52, 1971.
6. Whiteman J.R., Treatment of singularities in a harmonic mixed
boundary value problem by dual series methods.
Q.J.Mech.Appl.Math.21, 4l-50, 1968.
7. Whiteman J.R., Finite difference techniques for a harmonic mixed
boundary problem having a re-entrant boundary.
Proc.Roy.Soc.Lond.A.323, 271-276,1971.
8. Whiteman J.R. , Numerical solution of steady state diffusion problems
containing singularities, (to appear in The Finite Element
Method in Flow Problems, Gallagher, Oden, Taylor and Zienkiewicz(eds.)
Wiley, London).
9. Whiteman J.R. and Papamichael M. , Treatment of harmonic mixed boundary
problems by conformal transformation methods.
Z.angew. Math.Phys.23, 655-664, 1972.
10. Zlamal M. , Some recent advances in the mathematics of finite elements,
pp. 59-81 of J.R.Whiteman ( e d . ) , The Mathematics of Finite
Elements and Applications. Academic Press, London, 1973.