Math 5520 Lecture 18 Jeffrey Connors March 4, 2016 University of Connecticut

Math 5520 Lecture 18
Jeffrey Connors
University of Connecticut
March 4, 2016
Hermite cubic elements on triangles
There are four nodes zi with
I
z1 , z2 and z3 at the vertices
I
z4 at the center (z1 + z2 + z3 )/3
Lemma
Given a triangle T , a function v ∈ P3 (T ) is uniquely determined
by v (zi ) for i = 1, 2, 3, 4 and ∇v (zi ) for i = 1, 2, 3.
Proof. Let v (zi ) = 0 for i = 1, 2, 3, 4 and ∇v (zi ) = 0 for
i = 1, 2, 3. The restriction of v to an edge e of T may be
parameterized as a cubic function of a single variable. Call this
q(s) = v |e . Then
q(0) = q(1) = q 0 (0) = q 0 (1) = 0
implies that q(s) = 0 on e. This holds for all three edges, hence
we may assume that v is a cubic “bubble” function, which is
determined by the value at z4 , which is zero. Thus v = 0 on T .
Rectangular elements use tensor-product polynomials
We think of polynomials p ∈ Qk , where
p(x, y ) =
k X
k
X
αij x i y j .
i=0 j=0
Let τh = ∪E be the mesh; a union of mesh rectangles E . The
finite element spaces of order k are
S h = v ∈ C 0 (Ω) : v |E ∈ Qk (E ), ∀E ∈ τh .
I
Nodes are just a (k + 1) × (k + 1)-grid on each element.
I
Restriction to an edge is uniquely determined by k + 1 nodal
values per edge.
I
Continuity across edges follows immediately.
Shape functions for rectangles, case k = 1
We construct the basis functions φi for i = 1, 2, 3, 4 on a general
rectangle as follows:
φ1 (x, y ) = ψ1 (x)η1 (y ) φ2 (x, y ) = ψ2 (x)η1 (y )
φ3 (x, y ) = ψ2 (x)η2 (y ) φ4 (x, y ) = ψ1 (x)η2 (y )
where
ψ1 (x) =
η1 (y ) =
x2 −x
x2 −x1
y2 −y
y2 −y1
ψ2 (x) =
η2 (y ) =
x−x1
x2 −x1
y −y1
y2 −y1
The reference rectangle is [−1, 1] × [−1, 1]. The basis functions:
φ̂1 = 14 (1 − x̂)(1 − ŷ ) φ̂2 = 14 (1 + x̂)(1 − ŷ )
φ̂3 = 41 (1 + x̂)(1 + ŷ ) φ̂4 = 14 (1 − x̂)(1 + ŷ )
The reference mapping is easily handled for rectangles
Transformation:
x
=
y
=
1
1
(x2 − x1 )x̂ + (x2 + x1 )
2
2
1
1
(y2 − y1 )ŷ + (y2 + y1 )
2
2
The Jacobian is quite simple and |DF | = |E |/4 = |E |/|Ê |. It
follows
Z
Z
|E |
g (x, y ) dx dy =
g (x̂, ŷ ) d x̂ d ŷ .
4 Ê
E
The stiffness matrix uses
∂φi
∂x
∂φi
∂y
=
=
2 ∂ φ̂i
x2 − x1 ∂ x̂
2 ∂ φ̂i
y2 − y1 ∂ ŷ
General quadrilaterals have a wider range of applicability
I
For meshes on non-rectangular domains.
I
For Lagrangian methods where the grid deforms.
It turns out that defining the shape functions is non-trivial! Direct
usage of bilinear functions on quadrilaterals destroys continuity
across edges. The remedy:
1. Define a mapping F from the reference element Ê to
quadrilateral E .
2. Use bilinear basis functions φ̂i on Ê .
3. Invert F to get basis functions on E :
φi (x, y ) = φ̂i (x̂, ŷ ) = φ̂i (F −1 (x, y ))
The reference mapping is itself bilinear, but the basis
functions on E will not be
We take F to have the following form:
x
= a1 + a2 x̂ + a3 ŷ + a4 x̂ ŷ
y
= b1 + b2 x̂ + b3 ŷ + b4 x̂ ŷ
Now we require the mapping to transform corners of Ê to corners
of E , which requires one to solve two linear systems of the form
M~a = ~x ,
M ~b = ~y .
It is easy to find M −1 , so an explicit formula may be used for this
part. However, note that F −1 will not be bilinear, but rational.
This means that φ(x, y ) = φ̂i (F −1 (x, y )) is rational on E .
These rational basis functions provide continuity on the
general quadrilaterals
Let us choose two nodes on Ê , say ẑi and ẑj .
1. Let ê be the edge of Ê containing ẑi and ẑj .
2. Let P̂ be any point on ê.
3. Write P̂ = α1 ẑi + α2 ẑj with α1 + α2 = 1 and 0 ≤ αk ≤ 1, for
k = 1, 2.
4. It holds that P = F (P̂) satisfies P = α1 zi + α2 zj .
5. F −1 |e acts linearly as a map from e onto ê.
6. Now let φ(zi ) = φ(zj ) = 0. Then
φ̂(ẑi ) = φ̂(ẑj ) = 0
⇒ φ̂(P̂) = 0, ∀P̂ ∈ ê
⇒ φ(P) = 0, ∀P ∈ e.
It is convenient to use the reference mapping to calculate
matrix or load vector entries
The Jacobian of the reference mapping F : Ê → E is
a2 + a4 ŷ a3 + a4 x̂
J =
.
b2 + b4 ŷ b3 + b4 x̂
Load vector entries look like:
Z
Z
f (z)φi (z) dz =
f (F (ẑ))φ̂i (ẑ)|det(J )| d ẑ.
E
Ê
Mass matrix entries look like:
Z
Z
a0 (z)φj (z) φi (z) dz =
a0 (F (ẑ))φ̂j (ẑ)φ̂i (ẑ)|det(J )| d ẑ.
E
Ê
Stiffness matrix entries look like:
Z
A(z)∇φj (z) · ∇φi (z) dz =
E
Z
ˆ φ̂j (ẑ) · J −T ∇
ˆ φ̂i (ẑ) |det(J )| d ẑ.
A(F (ẑ)) J −T ∇
Ê
Brief aside: the next homework requires Richardson
extrapolation
If one knows the error on two grids, the convergence rate is readily
calculated. If you don’t have the true solution, it is nice if you can
do the following:
1. Solve on grids of size 4h, 2h and h.
2. Apply the error ansatz kuα − uβ k = C αp + C β p as follows:
C (4h)p + C (2h)p
4p + 2p
ku4h − u2h k
=
=
= 2p .
ku2h − uh k
C (2h)p + Chp
2p + 1
3. Take the logarithm of both sides to get p.
Note that this requires you to work with things defined on two
grids simultaneously to calculate the left side.