# Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 242 [warning: Documents this large are best viewed by clicking the View PDF Link!]

The Finite Element Method in Geodynamics
C. Thieulot
April 19, 2019
Contents
1 Introduction 6
1.1 Philosophy ............................................ 6
1.2 Acknowledgments......................................... 6
1.3 Essentialliterature........................................ 6
1.4 Installation ............................................ 6
1.5 Whatisaeldstone?....................................... 6
1.6 Why the Finite Element method? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7 Notations ............................................. 7
1.8 Colour maps for visualisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 List of tutorials 8
3 The physical equations 10
3.1 The heat transport equation - energy conservation equation . . . . . . . . . . . . . . . . . 10
3.2 The momentum conservation equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3 The mass conservation equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.4 The equations in ASPECT manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.5 the Boussinesq approximation: an Incompressible ﬂow . . . . . . . . . . . . . . . . . . . . 13
3.6 Stokes equation for elastic medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.7 The strain rate tensor in all coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . 15
3.7.1 Cartesiancoordinates .................................. 15
3.7.2 Polarcoordinates..................................... 15
3.7.3 Cylindrical coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.7.4 Spericalcoordinates ................................... 15
3.8 Boundaryconditions....................................... 16
3.8.1 TheStokesequations................................... 16
3.8.2 The heat transport equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.9 Meaningful physical quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4 The building blocks of the Finite Element Method 19
4.1 Numericalintegration ...................................... 19
4.1.1 in1D-theory ...................................... 19
4.1.2 in1D-examples..................................... 21
4.1.3 in2D/3D-theory .................................... 22
4.2 Themesh ............................................. 22
4.3 AbitofFEterminology..................................... 22
4.4 Elements and basis functions in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.4.1 Linear basis functions (Q1) ............................... 23
4.4.2 Quadratic basis functions (Q2) ............................. 23
4.4.3 Cubic basis functions (Q3) ............................... 24
4.5 Elements and basis functions in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.5.1 Bilinear basis functions in 2D (Q1)........................... 27
4.5.2 Biquadratic basis functions in 2D (Q2)......................... 29
1
4.5.3 Eight node serendipity basis functions in 2D (Qs
2) .................. 30
4.5.4 Bicubic basis functions in 2D (Q3) ........................... 31
4.5.5 Linear basis functions for triangles in 2D (P1)..................... 32
4.5.6 Quadratic basis functions for triangles in 2D (P2)................... 32
4.5.7 Cubic basis functions for triangles (P3) ........................ 33
4.6 Elements and basis functions in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.6.1 Linear basis functions in tetrahedra (P1)........................ 34
4.6.2 Triquadratic basis functions in 3D (Q2) ........................ 35
5 Solving the ﬂow equations with the FEM 36
5.1 strongandweakforms...................................... 36
5.2 Which velocity-pressure pair for Stokes? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.2.1 The compatibility condition (or LBB condition) . . . . . . . . . . . . . . . . . . . . 36
5.3 Families .............................................. 36
5.3.1 The bi/tri-linear velocity - constant pressure element (Q1×P0)........... 36
5.3.2 The bi/tri-quadratic velocity - discontinuous linear pressure element (Q2×P1) . 36
5.3.3 The bi/tri-quadratic velocity - bi/tri-linear pressure element (Q2×Q1) ...... 36
5.3.4 The stabilised bi/tri-linear velocity - bi/tri-linear pressure element (Q1×Q1-stab) 37
5.3.5 The MINI triangular element (P+
1×P1)........................ 37
5.3.6 The quadratic velocity - linear pressure triangle (P2×P1).............. 37
5.3.7 The Crouzeix-Raviart triangle (P+
2×P1) ...................... 37
5.4 Otherelements .......................................... 37
5.5 The penalty approach for viscous ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.6 The mixed FEM for viscous ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.6.1 inthreedimensions.................................... 40
5.6.2 Goingfrom3Dto2D .................................. 46
5.7 Solving the elastic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.8 Solving the heat transport equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6 Additional techniques and features 50
6.1 PicardandNewton........................................ 50
6.2 The SUPG formulation for the energy equation . . . . . . . . . . . . . . . . . . . . . . . . 50
6.3 Tracking materials and/or interfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.4 Dealingwithafreesurface.................................... 50
6.5 Convergence criterion for nonlinear iterations . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.6 Staticcondensation........................................ 50
6.7 The method of manufactured solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.7.1 AnalyticalbenchmarkI ................................. 51
6.7.2 Analytical benchmark II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.7.3 Analytical benchmark III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.7.4 Analytical benchmark IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.8 Assigning values to quadrature points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.9 Matrix(Sparse)storage ..................................... 58
6.9.1 2D domain - One degree of freedom per node . . . . . . . . . . . . . . . . . . . . . 58
6.9.2 2D domain - Two degrees of freedom per node . . . . . . . . . . . . . . . . . . . . 59
6.9.3 ineldstone........................................ 60
6.10Meshgeneration ......................................... 61
6.11Visco-Plasticity.......................................... 64
6.11.1 Tensorinvariants..................................... 64
6.11.2 Scalarviscoplasticity................................... 65
6.11.3 about the yield stress value Y.............................. 65
6.12Pressuresmoothing........................................ 66
6.13Pressurescaling.......................................... 68
6.14Pressurenormalisation...................................... 69
6.15Thechoiceofsolvers....................................... 70
6.15.1 The Schur complement approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
2
6.16TheGMRESapproach...................................... 74
6.17 The consistent boundary ﬂux (CBF) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.17.1 applied to the Stokes equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.17.2 applied to the heat equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.17.3 implementation - Stokes equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.18Thevalueofthetimestep .................................... 79
6.19mappings ............................................. 80
6.20 Exporting data to vtk format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.21Runge-Kuttamethods ...................................... 83
6.22AmIinornot?.......................................... 84
6.22.1 Two-dimensionalspace.................................. 84
6.22.2 Three-dimensional space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.23 Error measurements and convergence rates . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.24 The initial temperature ﬁeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.24.1 Single layer with imposed temperature b.c. . . . . . . . . . . . . . . . . . . . . . . 88
6.25 Single layer with imposed heat ﬂux b.c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.26 Single layer with imposed heat ﬂux and temperature b.c. . . . . . . . . . . . . . . . . . . 89
6.26.1 Halfcoolingspace .................................... 89
6.26.2 Platemodel........................................ 89
6.26.3 McKenzieslab ...................................... 89
6.27 Kinematic boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.27.1 In-out ﬂux boundary conditions for lithospheric models . . . . . . . . . . . . . . . 90
7fieldstone 01: simple analytical solution 91
8fieldstone 02: Stokes sphere 93
9fieldstone 03: Convection in a 2D box 94
10 fieldstone 04: The lid driven cavity 97
10.1 the lid driven cavity problem (ldc=0) ............................. 97
10.2 the lid driven cavity problem - regularisation I (ldc=1).................... 97
10.3 the lid driven cavity problem - regularisation II (ldc=2) ................... 97
11 fieldstone 05: SolCx benchmark 99
12 fieldstone 06: SolKz benchmark 101
13 fieldstone 07: SolVi benchmark 102
14 fieldstone 08: the indentor benchmark 104
15 fieldstone 09: the annulus benchmark 106
16 fieldstone 10: Stokes sphere (3D) - penalty 108
17 fieldstone 11: stokes sphere (3D) - mixed formulation 109
18 fieldstone 12: consistent pressure recovery 110
19 fieldstone 13: the Particle in Cell technique (1) - the eﬀect of averaging 112
20 fieldstone f14: solving the full saddle point problem 116
21 fieldstone f15: saddle point problem with Schur complement approach - benchmark
118
3
22 fieldstone f16: saddle point problem with Schur complement approach - Stokes
sphere 121
23 fieldstone 17: solving the full saddle point problem in 3D 123
23.0.1 Constantviscosity .................................... 125
23.0.2 Variableviscosity..................................... 126
24 fieldstone 18: solving the full saddle point problem with Q2×Q1elements 128
25 fieldstone 19: solving the full saddle point problem with Q3×Q2elements 130
26 fieldstone 20: the Busse benchmark 132
27 fieldstone 21: The non-conforming Q1×P0element 134
28 fieldstone 22: The stabilised Q1×Q1element 135
28.1 The Donea & Huerta benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
28.2 The Dohrmann & Bochev benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
28.3 The falling block experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
29 fieldstone 23: compressible ﬂow (1) - analytical benchmark 138
30 fieldstone 24: compressible ﬂow (2) - convection box 141
30.1Thephysics ............................................ 141
30.2Thenumerics ........................................... 141
30.3Theexperimentalsetup ..................................... 143
30.4Scaling............................................... 143
30.5Conservationofenergy1..................................... 144
30.5.1 under BA and EBA approximations . . . . . . . . . . . . . . . . . . . . . . . . . . 144
30.5.2 under no approximation at all . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
30.6Conservationofenergy2..................................... 145
30.7 The problem of the onset of convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
30.8 results - BA - Ra = 104..................................... 148
30.9 results - BA - Ra = 105..................................... 150
30.10results - BA - Ra = 106..................................... 151
30.11results - EBA - Ra = 104.................................... 152
30.12results - EBA - Ra = 105.................................... 154
30.13Onsetofconvection........................................ 155
31 fieldstone 25: Rayleigh-Taylor instability (1) - instantaneous 157
32 fieldstone 26: Slab detachment benchmark (1) - instantaneous 159
33 fieldstone 27: Consistent Boundary Flux 161
34 fieldstone 28: convection 2D box - Tosi et al, 2015 165
34.0.1 Case 0: Newtonian case, a la Blankenbach et al., 1989 . . . . . . . . . . . . . . . . 166
34.0.2 Case1........................................... 167
34.0.3 Case2........................................... 169
34.0.4 Case3........................................... 171
34.0.5 Case4........................................... 173
34.0.6 Case5........................................... 175
35 fieldstone 29: open boundary conditions 177
4
36 fieldstone 30: conservative velocity interpolation 180
36.1Couetteow............................................ 180
36.2SolCx ............................................... 180
36.3Streamlineow.......................................... 180
37 fieldstone 31: conservative velocity interpolation 3D 181
38 fieldstone 32: 2D analytical sol. from stream function 182
38.1Backgroundtheory........................................ 182
38.2Asimpleapplication ....................................... 182
39 fieldstone 33: Convection in an annulus 186
40 fieldstone 34: the Cartesian geometry elastic aquarium 188
41 fieldstone 35: 2D analytical sol. in annulus from stream function 190
41.2 No slip boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
41.3 Free slip boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
42 fieldstone 36: the annulus geometry elastic aquarium 196
43 fieldstone 37: marker advection and population control 199
44 fieldstone 38: Critical Rayleigh number 200
45 fieldstone: Gravity: buried sphere 202
46 Problems, to do list and projects for students 204
A Three-dimensional applications 206
B Codes in geodynamics 207
C Matrix properties 209
C.1 Symmetricmatrices ....................................... 209
C.2 Schurcomplement ........................................ 209
5
WARNING: this is work in progress
1 Introduction
1.1 Philosophy
This document was writing with my students in mind, i.e. 3rd and 4th year Geology/Geophysics stu-
dents at Utrecht University. I have chosen to use jargon as little as possible unless it is a term that is
commonly found in the geodynamics literature (methods paper as well as application papers). There is
no mathematical proof of any theorem or statement I make. These are to be found in generic Numerical
Analysic, Finite Element and Linear Algebra books.
The codes I provide here are by no means optimised as I value code readability over code eﬃciency. I
have also chosen to avoid resorting to multiple code ﬁles or even functions to favour a sequential reading of
the codes. These codes are not designed to form the basis of a real life application: Existing open source
highly optimised codes shoud be preferred, such as ASPECT [297, 253], CITCOM, LAMEM, PTATIN,
PYLITH, ...
All kinds of feedback is welcome on the text (grammar, typos, ...) or on the code(s). You will have
my eternal gratitude if you wish to contribute an example, a benchmark, a cookbook.
All the python scripts and this document are freely available at
https://github.com/cedrict/fieldstone
1.2 Acknowledgments
I have beneﬁtted from many discussions, lectures, tutorials, coﬀee machine discussions, debugging ses-
sions, conference poster sessions, etc ... over the years. I wish to name these instrumental people in
particular and in alphabetic order: Wolfgang Bangerth, Jean Braun, Rens Elbertsen, Philippe Fullsack,
Menno Fraters, Anne Glerum, Timo Heister, Robert Myhill, John Naliboﬀ, E. Gerry Puckett, Lukas van
de Wiel, Arie van den Berg, Tom Weir, and the whole ASPECT family/team.
1.3 Essential literature
http://www-udc.ig.utexas.edu/external/becker/Geodynamics557.pdf
1.4 Installation
python3.6 -m pip install --user numpy scipy matplotlib
1.5 What is a ﬁeldstone?
Simply put, it is stone collected from the surface of ﬁelds where it occurs naturally. It also stands for the
bad acronym: nite element deformation of stones which echoes the primary application of these codes:
geodynamic modelling.
6
1.6 Why the Finite Element method?
The Finite Element Method (FEM) is by no means the only method to solve PDEs in geodynamics, nor
is it necessarily the best one. Other methods are employed very succesfully, such as the Finite Diﬀerence
Method (FDM), the Finite Volume Method (FVM), and to a lesser extent the Discrete Element Method
(DEM) [152, 153, 184], or the Element Free Galerkin Method (EFGM) .
1.7 Notations
Scalars such as temperature, density, pressure, etc ... are simply obtained in L
A
T
E
Xby using the math
mode: T,ρ,p. Although it is common to lump vectors and matrices/tensors together by using bold
fonts, I have decided in the interest of clarity to distinguish between those: vectors are denoted by an
arrow atop the quantity, e.g.
ν,g, while matrices and tensors are in bold M,σ, etc ...
Also I use the ·notation between two vectors to denote a dot product u ·v =uivior a matrix-vector
multiplication M·a =Mij aj. If there is no ·between vectors, it means that the result a
b=aibjis a
matrix. Case in point,
∇ ·
νis the velocity divergence while
1.8 Colour maps for visualisation
In an attempt to homogenise the ﬁgures obtained with ParaView, I have decided to use a ﬁxed colour
scale for each ﬁeld throughout this document. These colour scales were obtained from this link and are
Perceptually Uniform Colour Maps .
Field colour code
Velocity/displacement CET-D1A
Pressure CET-L17
Velocity divergence CET-L1
Density CET-D3
Strain rate CET-R2
Viscosity CET-R3
Temperature CET-D9
7
2 List of tutorials
tutorial number
element outer solver formulation physical problem
3D
temperature
time stepping
nonlinear
compressible
analytical benchmark
numerical benchmark
elastomechanics
1Q1×P0penalty analytical benchmark
2Q1×P0penalty Stokes sphere
3Q1×P0penalty Blankenbach et al., 1989 † †
4Q1×P0penalty Lid driven cavity
5Q1×P0penalty SolCx benchmark
6Q1×P0penalty SolKz benchmark
7Q1×P0penalty SolVi benchmark
8Q1×P0penalty Indentor
9Q1×P0penalty annulus benchmark
10 Q1×P0penalty Stokes sphere
11 Q1×P0full matrix mixed Stokes sphere
12 Q1×P0penalty analytical benchmark
+ consistent press recovery
13 Q1×P0penalty Stokes sphere
+ markers averaging
14 Q1×P0full matrix mixed analytical benchmark
15 Q1×P0Schur comp. CG mixed analytical benchmark
16 Q1×P0Schur comp. PCG mixed Stokes sphere
17 Q2×Q1full matrix mixed Burstedde benchmark
18 Q2×Q1full matrix mixed analytical benchmark
19 Q3×Q2full matrix mixed analytical benchmark
20 Q1×P0penalty Busse et al., 1993 † † †
21 Q1×P0R-T penalty analytical benchmark
22 Q1×Q1-
stab
full matrix mixed analytical benchmark
23 Q1×P0mixed analytical benchmark
24 Q1×P0mixed convection box † †
8
tutorial number
element outer solver formulation physical problem
3D
temperature
time stepping
nonlinear
compressible
analytical benchmark
numerical benchmark
elastomechanics
25 Q1×P0full matrix mixed Rayleigh-Taylor instability
26 Q1×P0full matrix mixed Slab detachment
27 Q1×P0full matrix mixed CBF benchmarks † †
28 Q1×P0full matrix mixed Tosi et al, 2015 † † †
29 Q1×P0full matrix mixed Open Boundary conditions
30 Q1, Q2X X Cons. Vel. Interp (cvi) † †
31 Q1, Q2X X Cons. Vel. Interp (cvi) † †
32 Q1×P0full matrix mixed analytical benchmark
33 Q1×P0penalty convection in annulus † † †
34 Q1elastic Cartesian aquarium † †
35
36 Q1elastic annulus aquarium † †
37 Q1, Q2X X population control, bmw test † †
38 Critical Rayleigh number
XX
Analytical benchmark means that an analytical solution exists while numerical benchmark means that a comparison with other code(s) has been carried
out.
9
3 The physical equations
Symbol meaning unit
tTime s
x, y, z Cartesian coordinates m
r, θ Polar coordinates m,-
r, θ, z Cylindrical coordinates m,-,m
r, θ, φ Spherical coordinates m,-,-
νvelocity vector m·s1
ρmass density kg/m3
ηdynamic viscosity Pa·s
λpenalty parameter Pa·s
Ttemperature K
∇· divergence operator m1
ppressure Pa
˙
ε(
ν) strain rate tensor s1
αthermal expansion coeﬃcient K1
kthermal conductivity W/(m ·K)
CpHeat capacity J/K
Hintrinsic speciﬁc heat production W/kg
βTisothermal compressibility Pa1
τdeviatoric stress tensor Pa
σfull stress tensor Pa
3.1 The heat transport equation - energy conservation equation
Let us start from the heat transport equation as shown in Schubert, Turcotte and Olson :
ρCp
DT
Dt αT Dp
Dt =
∇ · k
T+Φ+ρH
with D/Dt being the total derivatives so that
DT
Dt =T
t +
ν·
TDp
Dt =p
t +
ν·
p
Solving for temperature, this equation is often rewritten as follows:
ρCp
DT
Dt
∇ · k
T=αT Dp
Dt +Φ+ρH
A note on the shear heating term Φ: In many publications, Φ is given by Φ = τij jui=τ:
ν.
10
Φ = τij jui
= 2η˙εd
ij jui
= 2η1
2˙εd
ij jui+ ˙εd
jiiuj
= 2η1
2˙εd
ij jui+ ˙εd
ij iuj
= 2η˙εd
ij
1
2(jui+iuj)
= 2η˙εd
ij ˙εij
= 2η˙
εd:˙
ε
= 2η˙
εd:˙
εd+1
3(
∇ ·
ν)1
= 2η˙
εd:˙
εd+ 2η˙
εd:1(
∇ ·
ν)
= 2η˙
εd:˙
εd(1)
Finally
Φ = τ:
ν= 2η˙
εd:˙
εd= 2η( ˙εd
xx)2+ ( ˙εd
yy)2+ 2( ˙εd
xy)2
3.2 The momentum conservation equations
Because the Prandlt number is virtually zero in Earth science applications the Navier Stokes equations
reduce to the Stokes equation:
∇ · σ+ρg = 0
Since
σ=p1+τ
it also writes
p+
∇ · τ+ρg = 0
Using the relationship τ= 2η˙
εdwe arrive at
p+
∇ · (2η˙
εd) + ρg = 0
3.3 The mass conservation equations
The mass conservation equation is given by
Dρ
Dt +ρ
∇ ·
ν= 0
or,
ρ
t +
∇ · (ρ
ν)=0
In the case of an incompressible ﬂow, then ρ/∂t = 0 and
ρ= 0, i.e. Dρ/Dt = 0 and the remaining
equation is simply:
∇ ·
ν= 0
11
3.4 The equations in ASPECT manual
The following is lifted oﬀ the ASPECT manual. We focus on the system of equations in a d= 2- or
d= 3-dimensional domain Ω that describes the motion of a highly viscous ﬂuid driven by diﬀerences in
the gravitational force due to a density that depends on the temperature. In the following, we largely
follow the exposition of this material in Schubert, Turcotte and Olson .
Speciﬁcally, we consider the following set of equations for velocity u, pressure pand temperature T:
∇ · 2η˙
ε(
ν)1
3(
∇ ·
ν)1+
p=ρg in Ω,(2)
∇ · (ρv) = 0 in ,(3)
ρCpT
t +
ν·
T
∇ · k
T=ρH
+ 2η˙ε(v)1
3(
∇ ·
ν)1:˙ε(v)1
3(
∇ ·
ν)1(4)
+αT v·
pin Ω,
where ˙
ε(
ν) = 1
2(
ν+
νT) is the symmetric gradient of the velocity (often called the strain rate).
In this set of equations, (253) and (254) represent the compressible Stokes equations in which v=
v(x, t) is the velocity ﬁeld and p=p(x, t) the pressure ﬁeld. Both ﬁelds depend on space xand time
t. Fluid ﬂow is driven by the gravity force that acts on the ﬂuid and that is proportional to both the
density of the ﬂuid and the strength of the gravitational pull.
Coupled to this Stokes system is equation (255) for the temperature ﬁeld T=T(x, t) that contains
heat conduction terms as well as advection with the ﬂow velocity v. The right hand side terms of this
equation correspond to
internal heat production for example due to radioactive decay;
friction (shear) heating;
In order to arrive at the set of equations that ASPECT solves, we need to
neglect the p/∂t.WHY?
neglect the ρ/∂t .WHY?
from equations above.
—————————————-
Also, their deﬁnition of the shear heating term Φ is:
Φ = kB(·v)2+ 2η˙
εd:˙
εd
For many ﬂuids the bulk viscosity kBis very small and is often taken to be zero, an assumption known
as the Stokes assumption: kB=λ+ 2η/3 = 0. Note that ηis the dynamic viscosity and λthe second
viscosity. Also,
τ= 2η˙
ε+λ(·v)1
but since kB=λ+ 2η/3 = 0, then λ=2η/3 so
τ= 2η˙
ε2
3η(·v)1= 2η˙
εd
12
3.5 the Boussinesq approximation: an Incompressible ﬂow
[from aspect manual] The Boussinesq approximation assumes that the density can be considered constant
in all occurrences in the equations with the exception of the buoyancy term on the right hand side of
(253). The primary result of this assumption is that the continuity equation (254) will now read
·v= 0
This implies that the strain rate tensor is deviatoric. Under the Boussinesq approximation, the equations
are much simpliﬁed:
−∇ · [2η˙
ε(v)] + p=ρgin Ω,(5)
∇ · (ρv) = 0 in ,(6)
ρ0CpT
t +v· ∇T−∇·kT=ρH in Ω (7)
Note that all terms on the rhs of the temperature equations have disappeared, with the exception of the
source term.
13
3.6 Stokes equation for elastic medium
What follows is mostly borrowed from Becker & Kaus lecture notes.
The strong form of the PDE that governs force balance in a medium is given by
·σ+f=0
where σis the stress tensor and fis a body force.
The stress tensor is related to the strain tensor through the generalised Hooke’s law:
σij =X
kl
Cijklkl (8)
where Cis the fourth-order elastic tensor. In the case of an isotropic material, this relationship simpliﬁes
to
σij =λkkδij + 2µij or, σ=λ(·u)1+ 2µ(9)
where λis the Lam´e parameter and µis the shear modulus1. The term ·uis the isotropic dilation.
The strain tensor is related to the displacement as follows:
=1
2(u+uT)
The incompressibility (bulk modulus), K, is deﬁned as p=K·uwhere pis the pressure with
p=1
3T r(σ)
=1
3[λ(·u)T r+2µT r[]]
=1
3[λ(·u)3 + 2µ(·u)]
=[λ+2
3µ](·u) (10)
so that K=λ+2
3µ.
Remark : Eq. (8) and (9) are analogous to the ones that one has to solve in the context of viscous
ﬂow using the penalty method. In this case λis the penalty coeﬃcient, uis the velocity, and µis then
the dynamic viscosity.
The Lam´e parameter and the shear modulus are also linked to νthe poisson ratio, and E, Young’s
modulus:
λ=µ2ν
12ν=νE
(1 + ν)(1 2ν)with E= 2µ(1 + ν)
The shear modulus, expressed often in GPa, describes the material’s response to shear stress. The poisson
ratio describes the response in the direction orthogonal to uniaxial stress. The Young modulus, expressed
in GPa, describes the material’s strain response to uniaxial stress in the direction of this stress.
1It is also sometimes written G
14
3.7 The strain rate tensor in all coordinate systems
The strain rate tensor ˙
εis given by
˙
ε=1
2(
ν+
νT) (11)
3.7.1 Cartesian coordinates
˙εxx =u
x (12)
˙εyy =v
y (13)
˙εzz =w
z (14)
˙εyx = ˙εxy =1
2u
y +v
x (15)
˙εzx = ˙εxz =1
2u
z +w
x (16)
˙εzy = ˙εyz =1
2v
z +w
y (17)
3.7.2 Polar coordinates
˙εrr =vr
r (18)
˙εθθ =vr
r+1
r
vθ
θ (19)
˙εθr = ˙εrθ =1
2vθ
r vθ
r+1
r
vr
θ (20)
3.7.3 Cylindrical coordinates
http://eml.ou.edu/equation/FLUIDS/STRAIN/STRAIN.HTM
3.7.4 Sperical coordinates
˙εrr =vr
r (21)
˙εθθ =vr
r+1
r
vθ
θ (22)
˙εφφ =1
rsin θ
vφ
φ (23)
˙εθr = ˙εrθ =1
2r
r (vθ
r) + 1
r
vr
θ (24)
˙εφr = ˙ε=1
21
rsin θ
vr
φ +r
r (vφ
r)(25)
˙εφθ = ˙εθφ =1
2sin θ
r
θ (vφ
sin θ) + 1
rsin θ
vθ
φ (26)
15
3.8 Boundary conditions
In mathematics, the Dirichlet (or ﬁrst-type) boundary condition is a type of boundary condition, named
after Peter Gustav Lejeune Dirichlet. When imposed on an ODE or PDE, it speciﬁes the values that a
solution needs to take on along the boundary of the domain. Note that a Dirichlet boundary condition
may also be referred to as a ﬁxed boundary condition.
The Neumann (or second-type) boundary condition is a type of boundary condition, named after Carl
Neumann. When imposed on an ordinary or a partial diﬀerential equation, the condition speciﬁes the
values in which the derivative of a solution is applied within the boundary of the domain.
It is possible to describe the problem using other boundary conditions: a Dirichlet boundary condition
speciﬁes the values of the solution itself (as opposed to its derivative) on the boundary, whereas the Cauchy
boundary condition, mixed boundary condition and Robin boundary condition are all diﬀerent types of
combinations of the Neumann and Dirichlet boundary conditions.
3.8.1 The Stokes equations
You may ﬁnd the following terms in the computational geodynamics literature:
free surface: this means that no force is acting on the surface, i.e. σ·n =
0. It is usually used on
the top boundary of the domain and allows for topography evolution.
free slip:
ν·n = 0 and (σ·n)×n =
0. This condition ensures a frictionless ﬂow parallel to the
boundary where it is prescribed.
no slip: this means that the velocity (or displacement) is exactly zero on the boundary, i.e.
ν=
0.
prescribed velocity:
ν=
νbc
stress b.c.:
open .b.c.: see ﬁeldstone 29.
3.8.2 The heat transport equation
There are two types of boundary conditions for this equation: temperature boundary conditions (Dirichlet
boundary conditions) and heat ﬂux boundary conditions (Neumann boundary conditions).
16
3.9 Meaningful physical quantities
Velocity (m/s):
Root mean square velocity (m/s):
νrms =R|
ν|2d
Rd1/2
=1
VZ|
ν|2d1/2
(27)
In Cartesian coordinates, for a cuboid domain of size Lx ×Ly×Lz, the vrms is simply given by:
νrms = 1
LxLyLzZLx
0ZLy
0ZLz
0
(u2+v2+w2)dxdydz!1/2
(28)
In the case of an annulus domain, although calculations are carried out in Cartesian coordinates,
it makes sense to look at the radial velocity component vrand the tangential velocity component
vθ, and their respective root mean square averages:
vr|rms =1
VZ
v2
rd1/2
(29)
vθ|rms =1
VZ
v2
θd1/2
(30)
Pressure (Pa):
Stress (Pa):
Strain (X):
Strain rate (s1):
Rayleigh number (X):
Prandtl number (X):
Nusselt number (X): the Nusselt number (Nu) is the ratio of convective to conductive heat transfer
across (normal to) the boundary. The conductive component is measured under the same conditions
as the heat convection but with a (hypothetically) stagnant (or motionless) ﬂuid.
In practice the Nusselt number Nu of a layer (typically the mantle of a planet) is deﬁned as follows:
Nu = q
qc
(31)
where qis the heat transferred by convection while qc=kT/D is the amount of heat that would
be conducted through a layer of thickness Dwith a temperature diﬀerence ∆Tacross it with k
being the thermal conductivity.
For 2D Cartesian systems of size (Lx,Ly) the Nu is computed 
Nu =
1
LxRLx
0kT
y (x, y =Ly)dx
1
LxRLx
0kT (x, y = 0)/Lydx =LyRLx
0
T
y (x, y =Ly)dx
RLx
0T(x, y = 0)dx
i.e. it is the mean surface temperature gradient over the mean bottom temperature.
finish. not happy with definition. Look at literature
Note that in the case when no convection takes place then the measured heat ﬂux at the top is the
one obtained from a purely conductive proﬁle which yields Nu=1.
Note that a relationship Ra Nuαexists between the Rayleigh number Ra and the Nusselt number
Nu in convective systems, see  and references therein.
17
Turning now to cylindrical geometries with inner radius R1and outer radius R2, we deﬁne f=
R1/R2. A small value of fcorresponds to a high degree of curvature. We assume now that
R2R1= 1, so that R2= 1/(1 f) and R1=f/(1 f). Following , the Nusselt number at
the inner and outer boundaries are:
Nuinner =fln f
1f
1
2πZ2π
0T
r r=R1
(32)
Nuouter =ln f
1f
1
2πZ2π
0T
r r=R2
(33)
Note that a conductive geotherm in such an annulus between temperatures T1and T2is given by
Tc(r) = ln(r/R2)
ln(R1/R2)=ln(r(1 f))
ln f
so that Tc
r =1
r
1
ln f
We then ﬁnd:
Nuinner =fln f
1f
1
2πZ2π
0Tc
r r=R1
=fln f
1f
1
R1
1
ln f= 1 (34)
Nuouter =ln f
1f
1
2πZ2π
0Tc
r r=R2
=ln f
1f
1
R2
1
ln f= 1 (35)
As expected, the recovered Nusselt number at both boundaries is exactly 1 when the temperature
ﬁeld is given by a steady state conductive geotherm.
derive formula for Earth size R1 and R2
Temperature (K):
Viscosity (Pa.s):
Density (kg/m3):
Heat capacity cp(J.K1.kg1): It is the measure of the heat energy required to increase the
temperature of a unit quantity of a substance by unit degree.
Heat conductivity, or thermal conductivity k(W.m1.K1). It is the property of a material that
indicates its ability to conduct heat. It appears primarily in Fourier’s Law for heat conduction.
Heat diﬀusivity: κ=k/(ρcp) (m2.s1). Substances with high thermal diﬀusivity rapidly adjust
their temperature to that of their surroundings, because they conduct heat quickly in comparison
to their volumetric heat capacity or ’thermal bulk’.
thermal expansion α(K1): it is the tendency of a matter to change in volume in response to a
change in temperature.
check aspect manual The 2D cylindrical shell benchmarks by Davies et al. 5.4.12
18
4 The building blocks of the Finite Element Method
4.1 Numerical integration
As we will see later, using the Finite Element method to solve problems involves computing integrals
which are more often than not too complex to be computed analytically/exactly. We will then need to
compute them numerically.
[wiki] In essence, the basic problem in numerical integration is to compute an approximate solution
to a deﬁnite integral Zb
a
f(x)dx
to a given degree of accuracy. This problem has been widely studied and we know that if f(x) is a
smooth function, and the domain of integration is bounded, there are many methods for approximating
the integral to the desired precision.
There are several reasons for carrying out numerical integration.
The integrand f(x) may be known only at certain points, such as obtained by sampling. Some
embedded systems and other computer applications may need numerical integration for this reason.
A formula for the integrand may be known, but it may be diﬃcult or impossible to ﬁnd an an-
tiderivative that is an elementary function. An example of such an integrand is f(x) = exp(x2),
the antiderivative of which (the error function, times a constant) cannot be written in elementary
form.
It may be possible to ﬁnd an antiderivative symbolically, but it may be easier to compute a numerical
approximation than to compute the antiderivative. That may be the case if the antiderivative is
given as an inﬁnite series or product, or if its evaluation requires a special function that is not
available.
4.1.1 in 1D - theory
The simplest method of this type is to let the interpolating function be a constant function (a polynomial
of degree zero) that passes through the point ((a+b)/2, f((a+b)/2)).
This is called the midpoint rule or rectangle rule.
Zb
a
f(x)dx '(ba)f(a+b
2)
insert here figure
The interpolating function may be a straight line (an aﬃne function, i.e. a polynomial of degree 1)
passing through the points (a, f(a)) and (b, f(b)).
This is called the trapezoidal rule.
Zb
a
f(x)dx '(ba)f(a) + f(b)
2
insert here figure
For either one of these rules, we can make a more accurate approximation by breaking up the interval
[a, b] into some number n of subintervals, computing an approximation for each subinterval, then adding
up all the results. This is called a composite rule, extended rule, or iterated rule. For example, the
composite trapezoidal rule can be stated as
Zb
a
f(x)dx 'ba
n f(a)
2+
n1
X
k=1
f(a+kba
n) + f(b)
2!
where the subintervals have the form [kh, (k+ 1)h], with h= (ba)/n and k= 0,1,2, . . . , n 1.
19
a) b)
The interval [2,2] is broken into 16 sub-intervals. The blue lines correspond to the approximation of
the red curve by means of a) the midpoint rule, b) the trapezoidal rule.
There are several algorithms for numerical integration (also commonly called ’numerical quadrature’,
or simply ’quadrature’) . Interpolation with polynomials evaluated at equally spaced points in [a, b] yields
the NewtonCotes formulas, of which the rectangle rule and the trapezoidal rule are examples. If we allow
the intervals between interpolation points to vary, we ﬁnd another group of quadrature formulas, such
as the Gauss(ian) quadrature formulas. A Gaussian quadrature rule is typically more accurate than a
NewtonCotes rule, which requires the same number of function evaluations, if the integrand is smooth
(i.e., if it is suﬃciently diﬀerentiable).
An npoint Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule con-
structed to yield an exact result for polynomials of degree 2n1 or less by a suitable choice of the points
xiand weights wifor i= 1, . . . , n.
The domain of integration for such a rule is conventionally taken as [1,1], so the rule is stated as
Z+1
1
f(x)dx =
n
X
iq=1
wiqf(xiq)
In this formula the xiqcoordinate is the i-th root of the Legendre polynomial Pn(x).
It is important to note that a Gaussian quadrature will only produce good results if the function f(x)
is well approximated by a polynomial function within the range [1,1]. As a consequence, the method
is not, for example, suitable for functions with singularities.
Gauss-Legendre points and their weights.
As shown in the above table, it can be shown that the weight values must fulﬁll the following condition:
X
iq
wiq= 2 (36)
and it is worth noting that all quadrature point coordinates are symmetrical around the origin.
Since most quadrature formula are only valid on a speciﬁc interval, we now must address the problem
of their use outside of such intervals. The solution turns out to be quite simple: one must carry out a
change of variables from the interval [a, b] to [1,1].
We then consider the reduced coordinate r[1,1] such that
r=2
ba(xa)1
20
This relationship can be reversed such that when ris known, its equivalent coordinate x[a, b] can be
computed:
x=ba
2(1 + r) + a
From this it follows that
dx =ba
2dr
and then Zb
a
f(x)dx =ba
2Z+1
1
f(r)dr 'ba
2
n
X
iq=1
wiqf(riq)
4.1.2 in 1D - examples
example 1 Since we know how to carry out any required change of variables, we choose for simplicity
a=1, b= +1. Let us take for example f(x) = π. Then we can compute the integral of this function
over the interval [a, b] exactly:
I=Z+1
1
f(x)dx =πZ+1
1
dx = 2π
We can now use a Gauss-Legendre formula to compute this same integral:
Igq =Z+1
1
f(x)dx =
nq
X
iq=1
wiqf(xiq) =
nq
X
iq=1
wiqπ=π
nq
X
iq=1
wiq
| {z }
=2
= 2π
where we have used the property of the weight values of Eq.(36). Since the actual number of points was
never speciﬁed, this result is valid for all quadrature rules.
example 2 Let us now take f(x) = mx +pand repeat the same exercise:
I=Z+1
1
f(x)dx =Z+1
1
(mx +p)dx = [1
2mx2+px]+1
1= 2p
Igq =Z+1
1
f(x)dx=
nq
X
iq=1
wiqf(xiq)=
nq
X
iq=1
wiq(mxiq+p)= m
nq
X
iq=1
wiqxiq
| {z }
=0
+p
nq
X
iq=1
wiq
| {z }
=2
= 2p
since the quadrature points are symmetric w.r.t. to zero on the x-axis. Once again the quadrature is able
to compute the exact value of this integral: this makes sense since an n-point rule exactly integrates a
2n1 order polynomial such that a 1 point quadrature exactly integrates a ﬁrst order polynomial like
the one above.
example 3 Let us now take f(x) = x2. We have
I=Z+1
1
f(x)dx =Z+1
1
x2dx = [1
3x3]+1
1=2
3
and
Igq =Z+1
1
f(x)dx=
nq
X
iq=1
wiqf(xiq)=
nq
X
iq=1
wiqx2
iq
nq= 1: x(1)
iq = 0, wiq= 2. Igq = 0
nq= 2: x(1)
q=1/3, x(2)
q= 1/3, w(1)
q=w(2)
q= 1. Igq =2
3
It also works nq>2 !
21
4.1.3 in 2D/3D - theory
Let us now turn to a two-dimensional integral of the form
I=Z+1
1Z+1
1
f(x, y)dxdy
Igq '
nq
X
iq=1
nq
X
jq
f(xiq, yjq)wiqwjq
4.2 The mesh
4.3 A bit of FE terminology
We introduce here some terminology for eﬃcient element descriptions :
For triangles/tetrahedra, the designation Pm×Pnmeans that each component of the velocity is
approximated by continuous piecewise complete Polynomials of degree mand pressure by continuous
piecewise complete Polynomials of degree n. For example P2×P1means
ua1+a2x+a3y+a4xy +a5x2+a6y2
with similar approximations for v, and
pb1+b2x+b3y
Both velocity and pressure are continuous across element boundaries, and each triangular element
contains 6 velocity nodes and three pressure nodes.
For the same families, Pm×Pnis as above, except that pressure is approximated via piecewise
discontinuous polynomials of degree n. For instance, P2×P1is the same as P2P1except that
pressure is now an independent linear function in each element and therefore discontinuous at
element boundaries.
For quadrilaterals/hexahedra, the designation Qm×Qnmeans that each component of the velocity
is approximated by a continuous piecewise polynomial of degree min each direction on the quadri-
lateral and likewise for pressure, except that the polynomial is of degree n. For instance, Q2×Q1
means
ua1+a2x+a3y+a4xy +a5x2+a6y2+a7x2y+a8xy2+a9x2y2
and
pb1+b2x+b3y+b4xy
For these same families, Qm×Qnis as above, except that the pressure approximation is not
continuous at element boundaries.
Again for the same families, Qm×Pnindicates the same velocity approximation with a pressure
approximation that is a discontinuous complete piecewise polynomial of degree n(not of degree n
in each direction !)
The designation P+
mor Q+
mmeans that some sort of bubble function was added to the polynomial
approximation for the velocity. You may also ﬁnd the term ’enriched element’ in the literature.
Finally, for n= 0, we have piecewise-constant pressure, and we omit the minus sign for simplicity.
Another point which needs to be clariﬁed is the use of so-called ’conforming elements’ (or ’non-
conforming elements’). Following again , conforming velocity elements are those for which the
basis functions for a subset of H1for the continuous problem (the ﬁrst derivatives and their squares are
integrable in Ω). For instance, the rotated Q1×P0element of Rannacher and Turek (see section ??) is
such that the velocity is discontinous across element edges, so that the derivative does not exist there.
Another typical example of non-conforming element is the Crouzeix-Raviart element .
22
4.4 Elements and basis functions in 1D
4.4.1 Linear basis functions (Q1)
Let f(r) be a C1function on the interval [1 : 1] with f(1) = f1and f(1) = f2.
Let us assume that the function f(r) is to be approximated on [1,1] by the ﬁrst order polynomial
f(r) = a+br (37)
Then it must fulﬁll
f(r=1) = ab=f1
f(r= +1) = a+b=f2
a=1
2(f1+f2)b=1
2(f1+f2)
and then replacing a, b in Eq. (37) by the above values on gets
f(r) = 1
2(1 r)f1+1
2(1 + r)f2
or
f(r) =
2
X
i=1
Ni(r)f1
with
N1(r) = 1
2(1 r)
N2(r) = 1
2(1 + r) (38)
Let f(r) be a C1function on the interval [1 : 1] with f(1) = f1,f(0) = f2and f(1) = f3.
Let us assume that the function f(r) is to be approximated on [1,1] by the second order polynomial
f(r) = a+br +cr2(39)
Then it must fulﬁll
f(r=1) = ab+c=f1
f(r= 0) = a=f2
f(r= +1) = a+b+c=f3
23
a=f2b=1
2(f1 + f3) c=1
2(f1+f32f2)
and then replacing a, b, c in Eq. (39) by the above values on gets
f(r) = 1
2r(r1)f1+ (1 r2)f2+1
2r(r+ 1)f3
or,
f(r) =
3
X
i=1
Ni(r)fi
with
N1(r) = 1
2r(r1)
N2(r) = (1 r2)
N3(r) = 1
2r(r+ 1) (40)
4.4.3 Cubic basis functions (Q3)
The 1D basis polynomial is given by
f(r) = a+br +cr2+dr3
with the nodes at position -1,-1/3, +1/3 and +1.
f(1) = ab+cd=f1
f(1/3) = ab
3+c
9d
27 =f2
f(+1/3) = ab
3+c
9d
27 =f3
f(+1) = a+b+c+d=f4
Adding the ﬁrst and fourth equation and the second and third, one arrives at
f1+f4= 2a+ 2c f2+f3= 2a+2c
9
and ﬁnally:
a=1
16 (f1+ 9f2+ 9f3f4)
c=9
16 (f1f2f3+f4)
Combining the original 4 equations in a diﬀerent way yields
2b+ 2d=f4f1
2b
3+2d
27 =f3f2
so that
b=1
16 (f127f2+ 27f3f4)
d=9
16 (f1+ 3f23f3+f4)
Finally,
24
f(r) = a+b+cr2+dr3
=1
16(1 + r+ 9r29r3)f1
+1
16(9 27r9r2+ 27r3)f2
+1
16(9 + 27r9r227r3)f3
+1
16(1r+ 9r2+ 9r3)f4
=
4
X
i=1
Ni(r)fi
where
N1=1
16(1 + r+ 9r29r3)
N2=1
16(9 27r9r2+ 27r3)
N3=1
16(9 + 27r9r227r3)
N4=1
16(1r+ 9r2+ 9r3)
Veriﬁcation:
Let us assume f(r) = C, then
ˆ
f(r) = XNi(r)fi=X
i
NiC=CX
i
Ni=C
so that a constant function is exactly reproduced, as expected.
Let us assume f(r) = r, then f1=1, f2=1/3, f3= 1/3 and f4= +1. We then have
ˆ
f(r) = XNi(r)fi
=N1(r)1
3N2(r) + 1
3N3(r) + N4(r)
= [(1 + r+ 9r29r3)
1
3(9 27r9r227r3)
+1
3(9 + 27r9r2+ 27r3)
+(1r+ 9r2+ 9r3)]/16
= [r+ 9r+ 9rr]/16 + ...0...
=r(41)
The basis functions derivative are given by
25
N1
r =1
16(1 + 18r27r2)
N2
r =1
16(27 18r+ 81r2)
N3
r =1
16(+27 18r81r2)
N4
r =1
16(1 + 18r+ 27r2)
Veriﬁcation:
Let us assume f(r) = C, then
ˆ
f
r =X
i
Ni
r fi
=CX
i
Ni
r
=C
16[(1 + 18r27r2)
+(27 18r+ 81r2)
+(+27 18r81r2)
+(1 + 18r+ 27r2)]
= 0
Let us assume f(r) = r, then f1=1, f2=1/3, f3= 1/3 and f4= +1. We then have
ˆ
f
r =X
i
Ni
r fi
=1
16[(1 + 18r27r2)
1
3(27 18r+ 81r2)
+1
3(+27 18r81r2)
+(1 + 18r+ 27r2)]
=1
16[2 + 18 + 54r254r2]
= 1
4.5 Elements and basis functions in 2D
Let us for a moment consider a single quadrilateral element in the xy-plane, as shown on the following
ﬁgure:
26
Let us assume that we know the values of a given ﬁeld uat the vertices. For a given point Minside
the element in the plane, what is the value of the ﬁeld uat this point? It makes sense to postulate that
uM=u(xM, yM) will be given by
uM=φ(u1, u2, u3, u4, xM, yM)
where φis a function to be determined. Although φis not unique, we can decide to express the value
uMas a weighed sum of the values at the vertices ui. One option could be to assign all four vertices the
same weight, say 1/4 so that uM= (u1+u2+u3+u4)/4, i.e. uMis simply given by the arithmetic mean
of the vertices values. This approach suﬀers from a major drawback as it does not use the location of
point Minside the element. For instance, when (xM, yM)(x2, y2) we expect uMu2.
In light of this, we could now assume that the weights would depend on the position of Min a
continuous fashion:
u(xM, yM) =
4
X
i=1
Ni(xM, yM)ui
where the Niare continous (”well behaved”) functions which have the property:
Ni(xj, yj) = δij
or, in other words:
N3(x1, y1) = 0 (42)
N3(x2, y2) = 0 (43)
N3(x3, y3) = 1 (44)
N3(x4, y4) = 0 (45)
The functions Niare commonly called basis functions.
Omitting the Msubscripts for any point inside the element, the velocity components uand vare
given by:
ˆu(x, y) =
4
X
i=1
Ni(x, y)ui(46)
ˆv(x, y) =
4
X
i=1
Ni(x, y)vi(47)
Rather interestingly, one can now easily compute velocity gradients (and therefore the strain rate tensor)
since we have assumed the basis functions to be ”well behaved” (in this case diﬀerentiable):
˙xx(x, y) = u
x =
4
X
i=1
Ni
x ui(48)
˙yy(x, y) = v
y =
4
X
i=1
Ni
y vi(49)
˙xy (x, y) = 1
2
u
y +1
2
v
x =1
2
4
X
i=1
Ni
y ui+1
2
4
X
i=1
Ni
x vi(50)
How we actually obtain the exact form of the basis functions is explained in the coming section.
4.5.1 Bilinear basis functions in 2D (Q1)
In this section, we place ourselves in the most favorables case, i.e. the element is a square deﬁned by
1< r < 1, 1< s < 1 in the Cartesian coordinates system (r, s):
27
3===========2
| | (r_0,s_0)=(-1,-1)
| | (r_1,s_1)=(+1,-1)
| | (r_2,s_2)=(+1,+1)
| | (r_3,s_3)=(-1,+1)
| |
0===========1
This element is commonly called the reference element. How we go from the (x, y) coordinate system
to the (r, s) once and vice versa will be dealt later on. For now, the basis functions in the above reference
element and in the reduced coordinates system (r, s) are given by:
N1(r, s)=0.25(1 r)(1 s)
N2(r, s)=0.25(1 + r)(1 s)
N3(r, s)=0.25(1 + r)(1 + s)
N4(r, s)=0.25(1 r)(1 + s)
The partial derivatives of these functions with respect to rans sautomatically follow:
N1
r (r, s) = 0.25(1 s)N1
s (r, s) = 0.25(1 r)
N2
r (r, s) = +0.25(1 s)N2
s (r, s) = 0.25(1 + r)
N3
r (r, s) = +0.25(1 + s)N3
s (r, s) = +0.25(1 + r)
N4
r (r, s) = 0.25(1 + s)N4
s (r, s) = +0.25(1 r)
Let us go back to Eq.(47). And let us assume that the function v(r, s) = Cso that vi=Cfor
i= 1,2,3,4. It then follows that
ˆv(r, s) =
4
X
i=1
Ni(r, s)vi=C
4
X
i=1
Ni(r, s) = C[N1(r, s) + N2(r, s) + N3(r, s) + N4(r, s)] = C
This is a very important property: if the vfunction used to assign values at the vertices is constant, then
the value of ˆvanywhere in the element is exactly C. If we now turn to the derivatives of vwith respect
to rand s:
ˆv
r (r, s) =
4
X
i=1
Ni
r (r, s)vi=C
4
X
i=1
Ni
r (r, s) = C[0.25(1 s)+0.25(1 s)+0.25(1 + s)0.25(1 + s)] = 0
ˆv
s (r, s) =
4
X
i=1
Ni
s (r, s)vi=C
4
X
i=1
Ni
s (r, s) = C[0.25(1 r)0.25(1 + r)+0.25(1 + r)+0.25(1 r)] = 0
We reassuringly ﬁnd that the derivative of a constant ﬁeld anywhere in the element is exactly zero.
28
If we now choose v(r, s) = ar +bs with aand btwo constant scalars, we ﬁnd:
ˆv(r, s) =
4
X
i=1
Ni(r, s)vi(51)
=
4
X
i=1
Ni(r, s)(ari+bsi) (52)
=a
4
X
i=1
Ni(r, s)ri
| {z }
r
+b
4
X
i=1
Ni(r, s)si
| {z }
s
(53)
=a[0.25(1 r)(1 s)(1) + 0.25(1 + r)(1 s)(+1) + 0.25(1 + r)(1 + s)(+1) + 0.25(1 r)(1 + s)(1)]
+b[0.25(1 r)(1 s)(1) + 0.25(1 + r)(1 s)(1) + 0.25(1 + r)(1 + s)(+1) + 0.25(1 r)(1 + s)(+1)]
=a[0.25(1 r)(1 s)+0.25(1 + r)(1 s)+0.25(1 + r)(1 + s)0.25(1 r)(1 + s)]
+b[0.25(1 r)(1 s)0.25(1 + r)(1 s)+0.25(1 + r)(1 + s)+0.25(1 r)(1 + s)]
=ar +bs (54)
verify above eq. This set of bilinear shape functions is therefore capable of exactly representing a bilinear
ﬁeld. The derivatives are:
ˆv
r (r, s) =
4
X
i=1
Ni
r (r, s)vi(55)
=a
4
X
i=1
Ni
r (r, s)ri+b
4
X
i=1
Ni
r (r, s)si(56)
=a[0.25(1 s)(1) + 0.25(1 s)(+1) + 0.25(1 + s)(+1) 0.25(1 + s)(1)]
+b[0.25(1 s)(1) + 0.25(1 s)(1) + 0.25(1 + s)(+1) 0.25(1 + s)(+1)]
=a
4[(1 s) + (1 s) + (1 + s) + (1 + s)]
+b
4[(1 s)(1 s) + (1 + s)(1 + s)]
=a(57)
Here again, we ﬁnd that the derivative of the bilinear ﬁeld inside the element is exact: ˆv
r =v
r .
However, following the same methodology as above, one can easily prove that this is no more true
for polynomials of degree strivtly higher than 1. This fact has serious consequences: if the solution to
the problem at hand is for instance a parabola, the Q1shape functions cannot represent the solution
properly, but only by approximating the parabola in each element by a line. As we will see later, Q2
basis functions can remedy this problem by containing themselves quadratic terms.
4.5.2 Biquadratic basis functions in 2D (Q2)
This element is part of the so-called LAgrange family.
citation needed
Inside an element the local numbering of the nodes is as follows:
3=====6=====2
| | | (r_0,s_0)=(-1,-1) (r_4,s_4)=( 0,-1)
| | | (r_1,s_1)=(+1,-1) (r_5,s_5)=(+1, 0)
7=====8=====5 (r_2,s_2)=(+1,+1) (r_6,s_6)=( 0,+1)
| | | (r_3,s_3)=(-1,+1) (r_7,s_7)=(-1, 0)
| | | (r_8,s_8)=( 0, 0)
0=====4=====1
The basis polynomial is then
f(r, s) = a+br +cs +drs +er2+fs2+gr2s+hrs2+ir2s2
29
The velocity shape functions are given by:
N0(r, s) = 1
2r(r1)1
2s(s1)
N1(r, s) = 1
2r(r+ 1)1
2s(s1)
N2(r, s) = 1
2r(r+ 1)1
2s(s+ 1)
N3(r, s) = 1
2r(r1)1
2s(s+ 1)
N4(r, s) = (1 r2)1
2s(s1)
N5(r, s) = 1
2r(r+ 1)(1 s2)
N6(r, s) = (1 r2)1
2s(s+ 1)
N7(r, s) = 1
2r(r1)(1 s2)
N8(r, s) = (1 r2)(1 s2)
and their derivatives by:
N0
r =1
2(2r1)1
2s(s1) N0
s =1
2r(r1)1
2(2s1)
N1
r =1
2(2r+ 1)1
2s(s1) N1
s =1
2r(r+ 1)1
2(2s1)
N2
r =1
2(2r+ 1)1
2s(s+ 1) N2
s =1
2r(r+ 1)1
2(2s+ 1)
N3
r =1
2(2r1)1
2s(s+ 1) N3
s =1
2r(r1)1
2(2s+ 1)
N4
r = (2r)1
2s(s1) N4
s = (1 r2)1
2(2s1)
N5
r =1
2(2r+ 1)(1 s2)N5
s =1
2r(r+ 1)(2s)
N6
r = (2r)1
2s(s+ 1) N6
s = (1 r2)1
2(2s+ 1)
N7
r =1
2(2r1)(1 s2)N7
s =1
2r(r1)(2s)
N8
r = (2r)(1 s2)N8
s = (1 r2)(2s)
4.5.3 Eight node serendipity basis functions in 2D (Qs
2)
Inside an element the local numbering of the nodes is as follows:
3=====6=====2
| | | (r_0,s_0)=(-1,-1) (r_4,s_4)=( 0,-1)
| | | (r_1,s_1)=(+1,-1) (r_5,s_5)=(+1, 0)
7=====+=====5 (r_2,s_2)=(+1,+1) (r_6,s_6)=( 0,+1)
| | | (r_3,s_3)=(-1,+1) (r_7,s_7)=(-1, 0)
|||
0=====4=====1
The main diﬀerence with the Q2element resides in the fact that there is no node in the middle of the
element The basis polynomial is then
f(r, s) = a+br +cs +drs +er2+fs2+gr2s+hrs2
30
Note that absence of the r2s2term which was previously associated to the center node. We ﬁnd that
N0(r, s) = 1
4(1 r)(1 s)(rs1) (58)
N1(r, s) = 1
4(1 + r)(1 s)(rs1) (59)
N2(r, s) = 1
4(1 + r)(1 + s)(r+s1) (60)
N3(r, s) = 1
4(1 r)(1 + s)(r+s1) (61)
N4(r, s) = 1
2(1 r2)(1 s) (62)
N5(r, s) = 1
2(1 + r)(1 s2) (63)
N6(r, s) = 1
2(1 r2)(1 + s) (64)
N7(r, s) = 1
2(1 r)(1 s2) (65)
The shape functions at the mid side nodes are products of a second order polynomial parallel to side
and a linear function perpendicular to the side while shape functions for corner nodes are modiﬁcations
verify those
4.5.4 Bicubic basis functions in 2D (Q3)
Inside an element the local numbering of the nodes is as follows:
12===13===14===15 (r,s)_{00}=(-1,-1) (r,s)_{08}=(-1,+1/3)
|| || || || (r,s)_{01}=(-1/3,-1) (r,s)_{09}=(-1/3,+1/3)
08===09===10===11 (r,s)_{02}=(+1/3,-1) (r,s)_{10}=(+1/3,+1/3)
|| || || || (r,s)_{03}=(+1,-1) (r,s)_{11}=(+1,+1/3)
04===05===06===07 (r,s)_{04}=(-1,-1/3) (r,s)_{12}=(-1,+1)
|| || || || (r,s)_{05}=(-1/3,-1/3) (r,s)_{13}=(-1/3,+1)
00===01===02===03 (r,s)_{06}=(+1/3,-1/3) (r,s)_{14}=(+1/3,+1)
(r,s)_{07}=(+1,-1/3) (r,s)_{15}=(+1,+1)
The velocity shape functions are given by:
N1(r)=(1 + r+ 9r29r3)/16 N1(t)=(1 + t+ 9t29t3)/16
N2(r) = (+9 27r9r2+ 27r3)/16 N2(t) = (+9 27t9t2+ 27t3)/16
N3(r) = (+9 + 27r9r227r3)/16 N3(t) = (+9 + 27t9t227t3)/16
N4(r)=(1r+ 9r2+ 9r3)/16 N4(t)=(1t+ 9t2+ 9t3)/16
31
N01(r, s) = N1(r)N1(s)=(1 + r+ 9r29r3)/16 (1 + t+ 9s29s3)/16
N02(r, s) = N2(r)N1(s) = (+9 27r9r2+ 27r3)/16 (1 + t+ 9s29s3)/16
N03(r, s) = N3(r)N1(s) = (+9 + 27r9r227r3)/16 (1 + t+ 9s29s3)/16
N04(r, s) = N4(r)N1(s)=(1r+ 9r2+ 9r3)/16 (1 + t+ 9s29s3)/16
N05(r, s) = N1(r)N2(s) = (1 + r+ 9r29r3)/16 (9 27s9s2+ 27s3)/16
N06(r, s) = N2(r)N2(s) = (+9 27r9r2+ 27r3)/16 (9 27s9s2+ 27s3)/16
N07(r, s) = N3(r)N2(s) = (+9 + 27r9r227r3)/16 (9 27s9s2+ 27s3)/16
N08(r, s) = N4(r)N2(s)=(1r+ 9r2+ 9r3)/16 (9 27s9s2+ 27s3)/16
N09(r, s) = N1(r)N3(s) = (66)
N10(r, s) = N2(r)N3(s) = (67)
N11(r, s) = N3(r)N3(s) = (68)
N12(r, s) = N4(r)N3(s) = (69)
N13(r, s) = N1(r)N4(s) = (70)
N14(r, s) = N2(r)N4(s) = (71)
N15(r, s) = N3(r)N4(s) = (72)
N16(r, s) = N4(r)N4(s) = (73)
4.5.5 Linear basis functions for triangles in 2D (P1)
2
|\
| \ (r_0,s_0)=(0,0)
| \ (r_1,s_1)=(1,0)
| \ (r_2,s_2)=(0,2)
0=======1
The basis polynomial is then
f(r, s) = a+br +cs
and the shape functions:
N0(r, s)=1rs(74)
N1(r, s) = r(75)
N2(r, s) = s(76)
4.5.6 Quadratic basis functions for triangles in 2D (P2)
2
|\
| \ (r_0,s_0)=(0,0) (r_3,s_3)=(1/2,0)
5 4 (r_1,s_1)=(1,0) (r_4,s_4)=(1/2,1/2)
| \ (r_2,s_2)=(0,1) (r_5,s_5)=(0,1/2)
| \
0===3===1
The basis polynomial is then
f(r, s) = c1+c2r+c3s+c4r2+c5rs +c6s2
32
We have
f1=f(r1, s1) = c1
f2=f(r2, s2) = c1+c2+c4
f3=f(r3, s3) = c1+c3+c6
f4=f(r4, s4) = c1+c2/2 + c4/4
f5=f(r5, s5) = c1+c2/2 + c3/2
+c4/4 + c5/4 + c6/4
f6=f(r6, s6) = c1+c3/2 + c6/4
This can be cast as f=A·cwhere Ais a 6x6 matrix:
A=
100000
110100
101001
1 1/201/4 0 0
1 1/2 1/2 1/4 1/4 1/4
101/2 0 0 1/4
It is rather trivial to compute the inverse of this matrix:
A1=
1 0 0 0 0 0
31 0 4 0 0
3 0 1004
2204 0 0
4004 4 4
2 0 2 0 0 4
In the end, one obtains:
f(r, s) = f1+ (3f1f2+ 4f4)r+ (3f1f3+ 4f6)s
+(2f1+ 2f24f4)r2+ (4f14f4+ 4f54f6)rs
+(2f1+ 2f34f6)s2
=
6
X
i=1
Ni(r, s)fi(77)
with
N1(r, s)=13r3s+ 2r2+ 4rs + 2s2
N2(r, s) = r+ 2r2
N3(r, s) = s+ 2s2
N4(r, s)=4r4r24rs
N5(r, s)=4rs
N6(r, s)=4s4rs 4s2
4.5.7 Cubic basis functions for triangles (P3)
2
|\ (r_0,s_0)=(0,0) (r_5,s_5)=(2/3,1/3)
| \ (r_1,s_1)=(1,0) (r_6,s_6)=(1/3,2/3)
7 6 (r_2,s_2)=(0,1) (r_7,s_7)=(0,2/3)
| \ (r_3,s_3)=(1/3,0) (r_8,s_8)=(0,1/3)
8 9 5 (r_4,s_4)=(2/3,0) (r_9,s_9)=(1/3,1/3)
| \
0==3==4==1
33
The basis polynomial is then
f(r, s) = c1+c2r+c3s+c4r2+c5rs +c6s2+c7r3+c8r2s+c9rs2+c10s3
N0(r, s) = 9
2(1 rs)(1/3rs)(2/3rs) (78)
N1(r, s) = 9
2r(r1/3)(r2/3) (79)
N2(r, s) = 9
2s(s1/3)(s2/3) (80)
N3(r, s) = 27
2(1 rs)r(2/3rs) (81)
N4(r, s) = 27
2(1 rs)r(r1/3) (82)
N5(r, s) = 27
2rs(r1/3) (83)
N6(r, s) = 27
2rs(r2/3) (84)
N7(r, s) = 27
2(1 rs)s(s1/3) (85)
N8(r, s) = 27
2(1 rs)s(2/3rs) (86)
N9(r, s) = 27rs(1 rs) (87)
verify those
4.6 Elements and basis functions in 3D
4.6.1 Linear basis functions in tetrahedra (P1)
(r_0,s_0) = (0,0,0)
(r_1,s_1) = (1,0,0)
(r_2,s_2) = (0,2,0)
(r_3,s_3) = (0,0,1)
The basis polynomial is given by
f(r, s, t) = c0+c1r+c2s+c3t
f1=f(r1, s1, t1) = c0(88)
f2=f(r2, s2, t2) = c0+c1(89)
f3=f(r3, s3, t3) = c0+c2(90)
f4=f(r4, s4, t4) = c0+c3(91)
which yields:
c0=f1c1=f2f1c2=f3f1c3=f4f1
f(r, s, t) = c0+c1r+c2s+c3t
=f1+ (f2f1)r+ (f3f1)s+ (f4f1)t
=f1(1 rst) + f2r+f3s+f4t
=X
i
Ni(r, s, t)fi
Finally,
34
N1(r, s, t)=1rst
N2(r, s, t) = r
N3(r, s, t) = s
N4(r, s, t) = t
4.6.2 Triquadratic basis functions in 3D (Q2)
N1= 0.5r(r1) 0.5s(s1) 0.5t(t1)
N2= 0.5r(r+ 1) 0.5s(s1) 0.5t(t1)
N3= 0.5r(r+ 1) 0.5s(s+ 1) 0.5t(t1)
N4= 0.5r(r1) 0.5s(s+ 1) 0.5t(t1)
N5= 0.5r(r1) 0.5s(s1) 0.5t(t+ 1)
N6= 0.5r(r+ 1) 0.5s(s1) 0.5t(t+ 1)
N7= 0.5r(r+ 1) 0.5s(s+ 1) 0.5t(t+ 1)
N8= 0.5r(r1) 0.5s(s+ 1) 0.5t(t+ 1)
N9= (1.r2) 0.5s(s1) 0.5t(t1)
N10 = 0.5r(r+ 1) (1 s2) 0.5t(t1)
N11 = (1.r2) 0.5s(s+ 1) 0.5t(t1)
N12 = 0.5r(r1) (1 s2) 0.5t(t1)
N13 = (1.r2) 0.5s(s1) 0.5t(t+ 1)
N14 = 0.5r(r+ 1) (1 s2) 0.5t(t+ 1)
N15 = (1.r2) 0.5s(s+ 1) 0.5t(t+ 1)
N16 = 0.5r(r1) (1 s2) 0.5t(t+ 1)
N17 = 0.5r(r1) 0.5s(s1) (1 t2)
N18 = 0.5r(r+ 1) 0.5s(s1) (1 t2)
N19 = 0.5r(r+ 1) 0.5s(s+ 1) (1 t2)
N20 = 0.5r(r1) 0.5s(s+ 1) (1 t2)
N21 = (1 r2) (1 s2) 0.5t(t1)
N22 = (1 r2) 0.5s(s1) (1 t2)
N23 = 0.5r(r+ 1) (1 s2) (1 t2)
N24 = (1 r2) 0.5s(s+ 1) (1 t2)
N25 = 0.5r(r1) (1 s2) (1 t2)
N26 = (1 r2) (1 s2) 0.5t(t+ 1)
N27 = (1 r2) (1 s2) (1 t2)
35
5 Solving the ﬂow equations with the FEM
In the case of an incompressible ﬂow, we have seen that the continuity (mass conservation) equation takes
the simple form
·v = 0. In other word ﬂow takes place under the constraint that the divergence of its
velocity ﬁeld is exactly zero eveywhere (solenoidal constraint), i.e. it is divergence free.
We see that the pressure in the momentum equation is then a degree of freedom which is needed to
satisfy the incompressibilty constraint (and it is not related to any constitutive equation) . In other
words the pressure is acting as a Lagrange multiplier of the incompressibility constraint.
Various approaches have been proposed in the literature to deal with the incompressibility constraint
but we will only focus on the penalty method (section 5.5) and the so-called mixed ﬁnite element method
5.6.
5.1 strong and weak forms
The strong form consists of the governing equation and the boundary conditions, i.e. the mass, momentum
and energy conservation equations supplemented with Dirichlet and/or Neumann boundary conditions
on (parts of) the boundary.
To develop the ﬁnite element formulation, the partial diﬀerential equations must be restated in an
integral form called the weak form. In essence the PDEs are ﬁrst multiplied by an arbitrary function and
integrated over the domain.
5.2 Which velocity-pressure pair for Stokes?
The success of a mixed ﬁnite element formulation crucially depends on a proper choice of the local
interpolations of the velocity and the pressure.
5.2.1 The compatibility condition (or LBB condition)
5.3 Families
Taylor-Hood
USE  section 3.6.2
5.3.1 The bi/tri-linear velocity - constant pressure element (Q1×P0)
discussed in example 3.71 of 
5.3.2 The bi/tri-quadratic velocity - discontinuous linear pressure element (Q2×P1)
5.3.3 The bi/tri-quadratic velocity - bi/tri-linear pressure element (Q2×Q1)
This element, implemented in penalised form, is discussed in  and the follow-up paper . CHECK
Biquadratic velocities, bilinear pressure. See Hood and Taylor. The element satisﬁes the inf-sup
condition p215.
36
5.3.4 The stabilised bi/tri-linear velocity - bi/tri-linear pressure element (Q1×Q1-stab)
5.3.5 The MINI triangular element (P+
1×P1)
discussed in section 3.6.1 of 
5.3.6 The quadratic velocity - linear pressure triangle (P2×P1)
From . Taylor-Hood elements (Taylor and Hood 1973) are characterized by the fact that the
pressure is continuous in the region Ω. A typical example is the quadratic triangle (P2P1 element).
In this element the velocity is approximated by a quadratic polynomial and the pressure by a linear
polynomial. One can easily verify that both approximations are continuous over the element boundaries.
It can be shown, Segal (1979), that this element is admissible if at least 3 elements are used. The
quadrilateral counterpart of this triangle is the Q2×Q1element.
5.3.7 The Crouzeix-Raviart triangle (P+
2×P1)
From : seven-node Crouzeix-Raviart triangle with quadratic velocity shape functions enhanced
by a cubic bubble function and discontinuous linear interpolation for the pressure ﬁeld [e.g., Cuvelier et
al., 1986]. This element is stable and no additional stabilization techniques are required [Elman et al.,
2005].
From . These elements are characterized by a discontinuous pressure; discontinuous on element
boundaries. For output purposes (printing, plotting etc.) these discontinuous pressures are averaged in
vertices for all the adjoining elements.
The most simple Crouzeix-Raviart element is the non-conforming linear triangle with constant pressure
(P1×P0).
p248 of Elman book. satisﬁes LBB condition.
5.4 Other elements
P1P0 example 3.70 in 
37
P1P1
Q2Q2
P2P2
5.5 The penalty approach for viscous ﬂow
In order to impose the incompressibility constraint, two widely used procedures are available, namely the
Lagrange multiplier method and the penalty method [30, 259]. The latter is implemented in elefant,
which allows for the elimination of the pressure variable from the momentum equation (resulting in a
reduction of the matrix size).
Mathematical details on the origin and validity of the penalty approach applied to the Stokes problem
can for instance be found in ,  or .
The penalty formulation of the mass conservation equation is based on a relaxation of the incompress-
ibility constraint and writes
∇ ·
ν+p
λ= 0 (92)
where λis the penalty parameter, that can be interpreted (and has the same dimension) as a bulk
viscosity. It is equivalent to say that the material is weakly compressible. It can be shown that if one
chooses λto be a suﬃciently large number, the continuity equation
·
ν= 0 will be approximately
satisﬁed in the ﬁnite element solution. The value of λis often recommended to be 6 to 7 orders of
magnitude larger than the shear viscosity [144, 260].
Equation (92) can be used to eliminate the pressure in Eq. (??) so that the mass and momentum
conservation equations fuse to become :
∇ · (2η˙ε(
ν)) + λ
(
∇ ·
ν) = ρg= 0 (93)
 have established the equivalence for incompressible problems between the reduced integration of
the penalty term and a mixed Finite Element approach if the pressure nodes coincide with the integration
points of the reduced rule.
In the end, the elimination of the pressure unknown in the Stokes equations replaces the original
saddle-point Stokes problem  by an elliptical problem, which leads to a symmetric positive deﬁnite
(SPD) FEM matrix. This is the major beneﬁt of the penalized approach over the full indeﬁnite solver
with the velocity-pressure variables. Indeed, the SPD character of the matrix lends itself to eﬃcient
solving stragegies and is less memory-demanding since it is suﬃcient to store only the upper half of the
matrix including the diagonal  .
list codes which use this approach
The stress tensor σis symmetric (i.e. σij =σji). For simplicity I will now focus on a Stokes ﬂow in
two dimensions.
Since the penalty formulation is only valid for incompressible ﬂows, then ˙
=˙
dso that the dsuper-
script is ommitted in what follows. The stress tensor can also be cast in vector format:
σxx
σyy
σxy
=
p
p
0
+ 2η
˙xx
˙yy
˙xy
=λ
˙xx + ˙yy
˙xx + ˙yy
0
+ 2η
˙xx
˙yy
˙xy
=
λ
110
110
000
| {z }
K
+η
200
020
001
| {z }
C
·
u
x
v
y
u
y +v
x
Remember that
u
x =
4
X
i=1
Ni
x ui
v
y =
4
X
i=1
Ni
y vi
38
u
y +v
x =
4
X
i=1
Ni
y ui+
4
X
i=1
Ni
x vi
so that
u
x
v
y
u
y +v
x
=
N1
x 0N2
x 0N3
x 0N4
x 0
0N1
y 0N2
y 0N3
y 0N4
y
N1
y
N1
x
N2
y
N2
x
N3
y
N3
x
N3
y
N4
x
| {z }
B
·
u1
v1
u2
v2
u3
v3
u4
v4
| {z }
V
Finally,
σ =
σxx
σyy
σxy
= (λK+ηC)·B·V
We will now establish the weak form of the momentum conservation equation. We start again from
∇ · σ+
b=
0
For the Ni’s ’regular enough’, we can write:
Ze
Ni
∇ · σdΩ + Ze
NibdΩ=0
We can integrate by parts and drop the surface term2:
Ze
Ni·σdΩ = Ze
Nibd
or,
Ze
Ni
x 0Ni
y
0Ni
y
Ni
x
·
σxx
σyy
σxy
dΩ = Ze
Nibd
Let i= 1,2,3,4 and stack the resulting four equations on top of one another.
Ze
N1
x 0N1
y
0N1
y
N1
x
·
σxx
σyy
σxy
dΩ = Ze
N1bx
bydΩ (94)
Ze
N2
x 0N2
y
0N2
y
N2
x
·
σxx
σyy
σxy
dΩ = Ze
Nibx
bydΩ (95)
Ze
N3
x 0N3
y
0N3
y
N3
x
·
σxx
σyy
σxy
dΩ = Ze
N3bx
bydΩ (96)
Ze
N4
x 0N4
y
0N4
y
N4
x
·
σxx
σyy
σxy
dΩ = Ze
N4bx
bydΩ (97)
We easily recognize BTinside the integrals! Let us deﬁne
NT
b= (N1bx, N1by, ...N4bx, N4by)
2We will come back to this at a later stage
39
then we can write
Ze
BT·
σxx
σyy
σxy
dΩ = Ze
Nbd
and ﬁnally: Ze
BT·[λK+ηC]·B·VdΩ = Ze
Nbd
Since Vcontains the velocities at the corners, it does not depend on the xor ycoordinates so it can be
taking outside of the integral:
Ze
BT·[λK+ηC]·Bd
| {z }
Ael(8×8)
·V
|{z}
(8x1)
=Ze
Nbd
| {z }
Bel(8×1)
or,
Ze
λBT·K·Bd
| {z }
Aλ
el(8×8)
+Ze
ηBT·C·Bd
| {z }
Aη
el(8×8)
·V
|{z}
(8x1)
=Ze
Nbd
| {z }
Bel(8×1)
INTEGRATION - MAPPING
reduced integration 
1. partition domain Ω into elements Ωe,e= 1, ...nel.
2. loop over elements and for each element compute Ael,Bel
3. a node belongs to several elements
need to assemble Ael and Bel in A,B
4. apply boundary conditions
5. solve system: x=A1·B
6. visualise/analyse x
5.6 The mixed FEM for viscous ﬂow
5.6.1 in three dimensions
In what follows the ﬂow is assumed to be incompressible, isoviscous and isothermal.
The methodology to derive the discretised equations of the mixed system is quite similar to the one
we have used in the case of the penalty formulation. The big diﬀerence comes from the fact that we are
now solving for both velocity and pressure at the same time, and that we therefore must solve the mass
and momentum conservation equations together. As before, velocity inside an element is given by
νh(r) =
mv
X
i=1
Nν
i(r)
νi(98)
40
where Nv
iare the polynomial basis functions for the velocity, and the summation runs over the mvnodes
composing the element. A similar expression is used for pressure:
ph(r) =
mp
X
i=1
Np
i(r)pi(99)
Note that the velocity is a vector of size while pressure (and temperature) is a scalar. There are then ndofv
velocity degrees of freedom per node and ndofppressure degrees of freedom. It is also very important
to remember that the numbers of velocity nodes and pressure nodes for a given element are more often
than not diﬀerent and that velocity and pressure nodes need not be colocated. Indeed, unless co-called
’stabilised elements’ are used, we have mv> mp, which means that the polynomial order of the velocity
ﬁeld is higher than the polynomial order of the pressure ﬁeld (usually by value 1).
insert here link(s) to manual and literature
Other notations are sometimes used for Eqs.(98) and (99):
uh(r) =
Nν·u vh(r) =
Nν·v wh(r) =
Nν·w ph(r) =
Np·p (100)
where
ν= (u, v, w) and
Nνis the vector containing all basis functions evaluated at location r:
Nv=Nν
1(r), Nν
2(r), Nν
3(r),...Nν
mv(r)(101)
Np=Np
1(r), Np
2(r), Np
3(r),...Np
mp(r)(102)
and with
u = (u1, u2, u3,...umv) (103)
v = (v1, v2, v3,...vmv) (104)
w = (w1, w2, w3,...wmv) (105)
p =p1, p2, p3,...pmp(106)
We will now establish the weak form of the momentum conservation equation. We start again from
∇ · σ+
b=
0 (107)
∇ ·v = 0 (108)
For the Nν
i’s and Np
i’regular enough’, we can write:
Ze
Nν
i
∇ · σdΩ + Ze
Nν
i
b d =
0 (109)
Ze
Np
i
∇ ·vdΩ = 0 (110)
We can integrate by parts and drop the surface term3:
Ze
Nν
i·σdΩ = Ze
Nν
i
bdΩ (111)
Ze
Np
i
∇ ·vdΩ = 0 (112)
or,
Ze
Nν
i
x 0 0 Nν
i
y
Nν
i
z 0
0Nν
i
y 0Nν
i
x 0Nν
i
z
0 0 Nν
i
z 0Nν
i
x
Nν
i
y
·
σxx
σyy
σzz
σxy
σxz
σyz
dΩ = Ze
Nν
i
b dΩ (113)
3We will come back to this at a later stage
41
As before (see section XXX) the above equation can ultimately be written:
Ze
BT·
σxx
σyy
σzz
σxy
σxz
σyz
dΩ = Ze
NbdΩ (114)
We have previously established that the strain rate vector
˙εis:
˙ε=
u
x
v
y
w
z
u
y +v
x
u
z +w
x
v
z +w
y
=
P
i
Nν
i
x ui
P
i
Nν
i
y vi
P
i
Nν
i
z wi
P
i
(Nν
i
y ui+Nν
i
x vi)
P
i
(Nν
i
z ui+Nν
i
x wi)
P
i
(Nν
i
z vi+Nν
i
y wi)
=
Nν
1
x 0 0 ··· Nν
mv
x 0 0
0Nν
1
y 0··· 0Nν
mv
y 0
0 0 Nν
1
z ··· 0 0 Nν
mv
z
Nν
1
y
Nν
1
x 0··· Nν
mv
x
Nν
mv
x 0
Nν
1
z 0Nν
1
x ··· Nν
mv
z 0Nν
mv
x
0Nν
1
z
Nν
1
y ··· 0Nν
mv
z
Nν
mv
y
| {z }
B
·
u1
v1
w1
u2
v2
w2
u3
v3
. . .
umv
vmv
wmv
| {z }
~
V
(115)
or,
˙ε=B·
Vwhere Bis the gradient matrix and
Vis the vector of all vector degrees of freedom for the
element. The matrix Bis then of size 3 ×mvand the vector
Vis mvlong. we have
σxx =p+ 2η˙εd
xx (116)
σyy =p+ 2η˙εd
yy (117)
σzz =p+ 2η˙εd
zz (118)
σxy = 2η˙εd
xy (119)
σxz = 2η˙εd
xz (120)
σyz = 2η˙εd
yz (121)
Since we here only consider incompressible ﬂow, we have ˙
εd=˙
εso
σ =
1
1
1
0
0
0
p+C·
˙ε=
1
1
1
0
0
0
Np·
P+C·B·
V(122)
with
C=η
200000
020000
002000
000100
000010
000001
˙ε=
˙εxx
˙εyy
˙εzz
2 ˙εxy
2 ˙εxz
2 ˙εyz
(123)
42
Let us deﬁne matrix Npof size 6 ×mp:
Np=
1
1
1
0
0
0
Np=
Np
Np
Np
0
0
0
(124)
so that
σ =Np·
P+C·B·
V(125)
ﬁnally Ze
BT·[Np·
P+C·B·
V]dΩ = Ze
NbdΩ (126)
or, Ze
BT·Npd
| {z }
G
·
P+Ze
BT·C·Bd
| {z }
K
·
V=Ze
Nbd
| {z }
~
f
(127)
where the matrix Kis of size (mvndofv×mvndofv), and matrix Gis of size (mvndofv×mpndofp).
43
Turning now to the mass conservation equation:
0 = Ze
Np
∇ ·v d
=Ze
Np
mv
X
i=1 Nν
i
x ui+N ν
i
y vi+N ν
i
z wid
=Ze
Np
1mv
P
i=1
Nν
i
x ui+
mv
P
i=1
Nν
i
y vi+
mv
P
i=1
Nν
i
z wi
Np
2mv
P
i=1
Nν
i
x ui+
mv
P
i=1
Nν
i
y vi+
mv
P
i=1
Nν
i
z wi
Np
3mv
P
i=1
Nν
i
x ui+
mv
P
i=1
Nν
i
y vi+
mv
P
i=1
Nν
i
z wi
. . .
Np
mpmv
P
i=1
Nν
i
x ui+
mv
P
i=1
Nν
i
y vi+
mv
P
i=1
Nν
i
z wi
d
=Ze
Np
1Np
1Np
1000
Np
2Np
2Np
2000
Np
3Np
3Np
3000
.
.
..
.
..
.
..
.
..
.
..
.
.
Np
mpNp
mpNp
mp000
·
P
i
Nν
i
x ui
P
i
Nν
i
y vi
P
i
Nν
i
z wi
P
i
(Nν
i
y ui+Nν
i
x vi)
P
i
(Nν
i
z ui+Nν
i
x wi)
P
i
(Nν
i
z vi+Nν
i
y wi)
d
=Ze
Np
1Np
1Np
1000
Np
2Np
2Np
2000
Np
3Np
3Np
3000
.
.
..
.
..
.
..
.
..
.
..
.
.
Np
mpNp
mpNp
mp000
|{z }
Np
·
˙ε d
=ZNp·Bd·
V
=GT
e·
V(128)
Ultimately we obtain the following system for each element:
KeGe
GT
e0·
V
P=
fe
0
Such a matrix is then generated for each element and then must me assembled into the global F.E. matrix.
Note that in this case the elemental Stokes matrix is antisymmetric. One can also deﬁne the following
symmetric modiﬁed Stokes matrix:
KeGe
GT
e0·
V
P=
fe
0
This matrix is symmetric, but indeﬁnite. It is non-singular if ker(GT) = 0, which is the case if t he
compatibility condition holds.
44
CHECK: Matrix Kis the viscosity matrix. Its size is (ndofvNv)×(ndofvNv) where ndofvis
the number of velocity degrees of freedom per node (typically 1,2 or 3) and Nvis the number of velocity
nodes. The size of matrix Gis (ndofvNv)×(ndofpNp) where ndofp(= 1) is the number of velocity
degrees of freedom per node and Npis the number of pressure nodes. Conversely, the size of matrix GT
is (ndofpNp)×(ndofvNv). The size of the global FE matrix is N=ndofvNv+ndofpNpNote
that matrix Kis analogous to a discrete Laplacian operator, matrix Gto a discrete gradient operator,
and matrix GTto a discrete divergence operator.
On the physical dimensions of the Stokes matrix blocks We start from the Stokes equations:
p+
∇ · (2η˙
ε) + ρg= 0 (129)
∇ ·
ν= 0 (130)
The dimensions of the terms in the ﬁrst equation are: ML2T2. The blocks Kand Gstem from
the weak form which obtained by multiplying the strong form equations by the (dimensionless) basis
funstions and integrating over the domain, so that it follows that
[K·
V]=[G·
P]=[
f] = ML2T2L3=MLT 2
We can then easily deduce:
[K] = MT 1[G] = L2
On elemental level mass balance. Note that in what is above no assumption has been made about
whether the pressure basis functions are continuous or discontinuous from one element to another.
Indeed, as mentioned in , since the weak formulation of the momentum equation involves inte-
gration by parts of
p, the resulting weak form contains no derivatives of pressure. This introduces the
possibility of approximating it by functions (piecewise polynomials, of course) that are not C0-continuous,
and indeed this has been done and is quite popular/useful.
It is then worth noting that only discontinuous pressure elements assure an element-level mass balance
: if for instance Np
iis piecewise-constant on element e(of value 1), the elemental weak form of the
mass conservervation equation is
Ze
Np
i
∇ ·
ν=Ze
∇ ·
ν=ZΓe
n ·
ν= 0
One potentially unwelcome consequence of using discontinuous pressure elements is that they do not
possess uniquely deﬁned pressure on the element boundaries; they are dual valued there, and often
multi-valued at certain velocity nodes.
On the Cmatrix The relationship between deviatoric stress and deviatoric strain rate tensor is
τ= 2η˙
εd(131)
= 2η˙
ε1
3(
∇ ·v)1(132)
= 2η
˙εxx ˙εxy ˙εxz
˙εyx ˙εyy ˙εyz
˙εzx ˙εzy ˙εzz
1
3( ˙εxx + ˙εyy + ˙εzz)
100
010
001
(133)
=2
3η
2 ˙εxx ˙εyy ˙εzz 3 ˙εxy 3 ˙εxz
3 ˙εyx ˙εyy + 2 ˙εyy ˙εyy 3 ˙εyz
3 ˙εzx 3 ˙εzy ˙εxx ˙εyy2 ˙εzz
(134)
45
so that
τ =2
3η
2 ˙εxx ˙εyy ˙εzz
˙εyy + 2 ˙εyy ˙εyy
˙εxx ˙εyy + 2 ˙εzz
3 ˙εxy
3 ˙εxz
3 ˙εyz
=η
3
422000
2 4 2000
22 4 0 0 0
0 0 0 3 0 0
0 0 0 0 3 0
0 0 0 0 0 3
| {z }
Cd
·
˙εxx
˙εyy
˙εzz
2 ˙εxy
2 ˙εxz
2 ˙εyz
=Cd·
˙ε(135)
In the case where we assume incompressible ﬂow from the beginning, i.e. ε=εd, then
τ =η
200000
020000
002000
000100
000010
000001
| {z }
C
·
˙εxx
˙εyy
˙εzz
2 ˙εxy
2 ˙εxz
2 ˙εyz
=C·
˙ε(136)
On the ’forgotten’ surface terms
5.6.2 Going from 3D to 2D
The world is three-dimensional. However, for many diﬀerent reasons one may wish to solve problems
which are two-dimensional.
Following ASPECT manual, we will think of two-dimensional models in the following way:
We assume that the domain we want to solve on is a two-dimensional cross section (in the xy
plane) that extends inﬁnitely far in both negative and positive zdirection.
We assume that the velocity is zero in the zdirection and that all variables have no variation in
the zdirection.
As a consequence, two-dimensional models are three-dimensional ones in which the zcomponent of the
velocity is zero and so are all zderivatives. This allows to reduce the momentum conservation equations
from 3 equations to 2 equations. However, contrarily to what is often seen, the 3D deﬁnition of the
deviatoric strain rate remains, i.e. in other words:
˙
εd=˙
ε1
3(
∇ ·v)1(137)
and not 1/2. In light of all this, the full strain rate tensor and the deviatoric strain rate tensor in 2D are
given by:
ε=
˙εxx ˙εxy ˙εxz
˙εyx ˙εyy ˙εyz
˙εzx ˙εzy ˙εzz
=
u
x
1
2u
y +v
x 0
1
2u
y +v
x v
y 0
0 0 0
(138)
˙
εd=1
3
2u
x v
y
1
2u
y +v
x 0
1
2u
y +v
x u
x + 2v
y 0
0 0 u
x v
y
(139)
Although the bottom right term may be surprising, it is of no consequence when this expression of the
deviatoric strain rate is used in the Stokes equation:
∇ · 2η˙
εd=
46
FINISH!
In two dimensions the velocity is then
ν= (u, v) and the FEM building blocks and matrices are
simply:
˙ε=
˙εxx
˙εyy
2 ˙εxy
=
u
x
v
y
u
y +v
x
=
Nν
1
x 0Nν
2
x 0Nν
3
x 0. . . Nν
mv
x 0
0Nν
1
y 0Nν
2
y 0Nν
3
y . . . 0Nν
mv
x
Nν
1
y
Nν
1
x
Nν
2
y
Nν
2
x
Nν
3
y
Nν
3
x . . . Nν
mv
y
Nν
mv
x
| {z }
B
·
u1
v1
u2
v2
u3
v3
. . .
umv
vmv
| {z }
~
V
(140)
we have
σxx =p+ 2η˙εxx (141)
σyy =p+ 2η˙εyy (142)
σxy = +2η˙εxy (143)
so
σ =
1
1
0
p+C·
˙ε=
1
1
0
Np·
P+C·B·
V(144)
with
C=η
200
020
001
or C=η
3
42 0
240
0 0 3
(145)
check the right C
Finally the matrix Npis of size 3 ×mp:
Np=
1
1
0
Np=
Np
Np
0
(146)
5.7 Solving the elastic equations
5.8 Solving the heat transport equation
We start from the ’bare-bones’ heat transport equation (source terms are omitted):
ρcpT
t +
ν·
T=
∇ · k
T(147)
In what follows we assume that the velocity vield
νis known so that temperature is the only unknwon.
Let Nθbe the temperature basis functions so that the temperature inside an element is given by4:
Th(r) =
mT
X
i=1
Nθ(r)Ti=
Nθ·
T(148)
where
Tis a vector of length mTThe weak form is then
Z
Nθ
iρcpT
t +
ν·
TdΩ = Z
Nθ
i
∇ · k
T dΩ (149)
4the θsuperscript has been chosen to denote temperature so as to avoid confusion with the transpose operator
47
Z
Nθ
iρcp
T
t d
| {z }
I
+Z
Nθ
iρcp
ν·
T d
| {z }
II
=Z
Nθ
i
∇ · k
T d
| {z }
III
i= 1, mT
Looking at the ﬁrst term:
Z
Nθ
iρcp
T
t dΩ = Z
Nθ
iρcp
Nθ·˙
T dΩ (150)
(151)
so that when we assemble all contributions for i= 1, mTwe get:
I=Z
Nθρcp
Nθ·˙
T dΩ = Z
ρcp
Nθ
Nθd·˙
T=MT·˙
T
where MTis the mass matrix of the system of size (mT×mT) with
MT
ij =Z
ρcpNθ
iNθ
jd
Turning now to the second term:
Z
Nθ
iρcp
ν·
T dΩ = Z
Nθ
iρcp(uT
x +vT
y )dΩ (152)
=Z
Nθ
iρcp(u
Nθ
x +v
Nθ
y )·
T dΩ (153)
(154)
so that when we assemble all contributions for i= 1, mTwe get:
II = Z
ρcp
Nθ(u
Nθ
x +v
Nθ
y )d!·
T=Ka·
T
where Kais the advection term matrix of size (mT×mT) with
(Ka)ij =Z
ρcpNθ
i uNθ
j
x +vN θ
j
y !d
Now looking at the third term, we carry out an integration by part and neglect the surface term for now,
so that
Z
Nθ
i
∇ · k
T dΩ = Z
k
Nθ
i·
T dΩ (155)
=Z
k
Nθ
i·
(
Nθ·
T)dΩ (156)
(157)
with
Nθ=
xNθ
1xNθ
2. . . xNθ
mT
yNθ
1yNθ
2. . . yNθ
mT
so that ﬁnally:
III =Z
k(
Nθ)T·
Nθd·
T=Kd·
T
where Kdis the diﬀusion term matrix:
Kd=Z
k(
Nθ)T·
Nθd
48
Ultimately terms I, II, III together yield:
Mθ·˙
T+ (Ka+Kd)·
T=
0
What now remains to be done is to address the time derivative on the temperature vector. The most
simple approach would be to use an explicit Euler one, i.e.:
T
t =
T(k)
T(k1)
δt
where
T(k)is the temperature ﬁeld at time step kand δt is the time interval between two consecutive
time steps. In this case the discretised heat transport equation is:
Mθ+ (Ka+Kd)δt·
T(k)=Mθ·
T(k1)
other time discr schemes!
49
6.1 Picard and Newton
6.2 The SUPG formulation for the energy equation
6.3 Tracking materials and/or interfaces
6.4 Dealing with a free surface
6.5 Convergence criterion for nonlinear iterations
6.6 Static condensation
50
6.7 The method of manufactured solutions
The method of manufactured solutions is a relatively simple way of carrying out code veriﬁcation. In
essence, one postulates a solution for the PDE at hand (as well as the proper boundary conditions),
inserts it in the PDE and computes the corresponding source term. The same source term and boundary
conditions will then be used in a numerical simulation so that the computed solution can be compared
with the (postulated) true analytical solution.
Examples of this approach are to be found in [142, 96, ?].
python codes/ﬁeldstone
python codes/ﬁeldstone burstedde
6.7.1 Analytical benchmark I
Taken from . We consider a two-dimensional problem in the square domain Ω = [0,1] ×[0,1],
which possesses a closed-form analytical solution. The problem consists of determining the velocity ﬁeld
ν= (u, v) and the pressure psuch that
η
ν+
p+
b=
0 in Ω
·v= 0 in Ω
v=0on Γ
where the ﬂuid viscosity is taken as η= 1. The components of the body force bare prescribed as
bx= (12 24y)x4+ (24 + 48y)x3+ (48y+ 72y248y3+ 12)x2
+(2 + 24y72y2+ 48y3)x+ 1 4y+ 12y28y3
by= (8 48y+ 48y2)x3+ (12 + 72y72y2)x2
+(4 24y+ 48y248y3+ 24y4)x12y2+ 24y312y4
With this prescribed body force, the exact solution is
u(x, y) = x2(1 x)2(2y6y2+ 4y3)
v(x, y) = y2(1 y)2(2x6x2+ 4x3)
p(x, y) = x(1 x)1/6
Note that the pressure obeys Rp d = 0 One can turn to the spatial derivatives of the ﬁelds:
u
x = (2x6x2+ 4x3)(2y6y2+ 4y3) (158)
v
y =(2x6x2+ 4x3)(2y6y2+ 4y3) (159)
with of course
∇ ·
ν= 0 and
p
x = 1 2x(160)
p
y = 0 (161)
The velocity and pressure ﬁelds look like:
51
http://ww2.lacan.upc.edu/huerta/exercises/Incompressible/Incompressible Ex1.htm
As shown in , If the LBB condition is not satisﬁed, spurious oscillations spoil the pressure
approximation. Figures below show results obtained with a mesh of 20x20 Q1P0 (left) and P1P1 (right)
elements:
]]
http://ww2.lacan.upc.edu/huerta/exercises/Incompressible/Incompressible Ex1.htm
Taking into account that the proposed problem has got analytical solution, it is easy to analyze
convergence of the diﬀerent pairs of elements:
http://ww2.lacan.upc.edu/huerta/exercises/Incompressible/Incompressible Ex1.htm
6.7.2 Analytical benchmark II
Taken from . It is for a unit square with ν=µ/ρ = 1 and the smooth exact solution is
u(x, y) = x+x22xy +x33xy2+x2y(162)
v(x, y) = y2xy +y23x2y+y3xy2(163)
p(x, y) = xy +x+y+x3y24/3 (164)
52
Note that the pressure obeys Rp dΩ=0
bx=(1 + y3x2y2) (165)
by=(1 3x2x3y) (166)
6.7.3 Analytical benchmark III
This benchmark begins by postulating a polynomial solution to the 3D Stokes equation :
v=
x+x2+xy +x3y
y+xy +y2+x2y2
2z3xz 3yz 5x2yz
(167)
and
p=xyz +x3y3z5/32 (168)
While it is then trivial to verify that this velocity ﬁeld is divergence-free, the corresponding body force
of the Stokes equation can be computed by inserting this solution into the momentum equation with a
given viscosity µ(constant or position/velocity/strain rate dependent). The domain is a unit cube and
velocity boundary conditions simply use Eq. (229). Note that the pressure fulﬁlls
Z
p(r)dΩ=0.
Constant viscosity In this case, the right hand side writes:
f=p+µ
2+6xy
2+2x2+ 2y2
10yz
=
yz + 3x2y3z
xz + 3x3y2z
xy +x3y3
+µ
2+6xy
2+2x2+ 2y2
10yz
We can compute the components of the strainrate tensor:
˙εxx = 1 + 2x+y+ 3x2y(169)
˙εyy = 1 + x+ 2y+ 2x2y(170)
˙εzz =23x3y5x2y(171)
˙εxy =1
2(x+y+ 2xy2+x3) (172)
˙εxz =1
2(3z10xyz) (173)
˙εyz =1
2(3z5x2z) (174)
Note that we of course have ˙εxx + ˙εyy + ˙εzz = 0.
Variable viscosity In this case, the right hand side is obtains through
f=p+µ
2+6xy
2+2x2+ 2y2
10yz
+
2 ˙εxx
2 ˙εxy
2 ˙εxz
µ
x +
2 ˙εxy
2 ˙εyy
2 ˙εyz
µ
y +
2 ˙εxz
2 ˙εyz
2 ˙εzz
µ
z (175)
53
The viscosity can be chosen to be a smooth varying function:
µ=exp(1 β(x(1 x) + y(1 y) + z(1 z))) (176)
Choosing β= 0 yields a constant velocity µ=e1(and greatly simpliﬁes the right-hand side). One can
easily show that the ratio of viscosities µ?in the system follows µ?= exp(3β/4) so that choosing β= 10
yields µ?'1808 and β= 20 yields µ?'3.269 ×106. In this case
µ
x =4β(1 2x)µ(x, y, z) (177)
µ
y =4β(1 2y)µ(x, y, z) (178)
µ
z =4β(1 2z)µ(x, y, z) (179)
 has carried out this benchmark for β= 4, i.e.:
µ(x, y, z) = exp(1 4(x(1 x) + y(1 y) + z(1 z)))
In a unit cube, this yields a variable viscosity such that 0.1353 < µ < 2.7182, i.e. a ratio of approx. 20
within the domain. We then have:
µ
x =4(1 2x)µ(x, y, z) (180)
µ
y =4(1 2y)µ(x, y, z) (181)
µ
z =4(1 2z)µ(x, y, z) (182)
sort out mess wrt Eq 26 of busa13
6.7.4 Analytical benchmark IV
From . The two-dimensional domain is a unit square. The body forces are:
fx= 128[x2(x1)212(2y1) + 2(y1)(2y1)y(12x212x+ 2)]
fy= 128[y2(y1)212(2x1) + 2(x1)(2x1)y(12y212y+ 2)]
(183)
The solution is
u=256x2(x1)2y(y1)(2y1)
v= 256x2(y1)2x(x1)(2x1)
p= 0 (184)
Another choice:
fx= 128[x2(x1)212(2y1) + 2(y1)(2y1)y(12x212x+ 2)] + y1/2
fy= 128[y2(y1)212(2x1) + 2(x1)(2x1)y(12y212y+ 2)] + x1/2
(185)
The solution is
u=256x2(x1)2y(y1)(2y1)
v= 256x2(y1)2x(x1)(2x1)
p= (x1/2)(y1/2) (186)
54
6.8 Assigning values to quadrature points
As we have seen in Section ??, the building of the elemental matrix and rhs requires (at least) to assign
a density and viscosity value to each quadrature point inside the element. Depending on the type of
modelling, this task can prove more complex than one might expect and have large consequences on the
solution accuracy.
Here are several options:
The simplest way (which is often used for benchmarks) consists in computing the ’real’ coordinates
(xq, yq, zq) of a given quadrature point based on its reduced coordinates (rq, sq, tq), and passing these
coordinates to a function which returns density and/or viscosity at this location. For instance, for
the Stokes sphere:
def rho(x,y):
if (x-.5)**2+(y-0.5)**2<0.123**2:
val=2.
else:
val=1.
return val
def mu(x,y):
if (x-.5)**2+(y-0.5)**2<0.123**2:
val=1.e2
else:
val=1.
return val
This is very simple, but it has been shown to potentially be problematic. In essence, it can introduce
of  and/or have a look at the section titled ”Averaging material properties” in the ASPECT
manual.
another similar approach consists in assigning a density and viscosity value to the nodes of the FE
mesh ﬁrst, and then using these nodal values to assign values to the quadrature points. Very often
,and quite logically, the shape functions are used to this eﬀect. Indeed we have seen before that for
any point (r, s, t) inside an element we have
fh(r, s, t) =
m
X
i
fiNi(r, s, t)
where the fiare the nodal values and the Nithe corresponding basis functions.
In the case of linear elements (Q1basis functions), this is straightforward. In fact, the basis functions
Nican be seen as moving weights: the closer the point is to a node, the higher the weight (basis
function value).
However, this is quite another story for quadratic elements (Q2basis functions). In order to
illustrate the problem, let us consider a 1D problem. The basis functions are
N1(r) = 1
2r(r1) N2(r)=1r2N3(r) = 1
2r(r+ 1)
Let us further assign: ρ1=ρ2= 0 and ρ3= 1. Then
ρh(r) =
m
X
i
ρiNi(r) = N3(r)
There lies the core of the problem: the N3(r) basis function is negative for r[1,0]. This means
that the quadrature point in this interval will be assigned a negative density, which is nonsensical
and numerically problematic!
use 2X Q1. write about it !
55
The above methods work ﬁne as long as the domain contains a single material. As soon as there are
multiple ﬂuids in the domain a special technique is needed to track either the ﬂuids themselves or their
interfaces. Let us start with markers. We are then confronted to the infernal trio (a menage a trois?)
which is present for each element, composed of its nodes, its markers and its quadrature points.
Each marker carries the material information (density and viscosity). This information must ul-
timately be projected onto the quadrature points. Two main options are possible: an algorithm is
designed and projects the marker-based ﬁelds onto the quadrature points directly or the marker ﬁelds
are ﬁrst projected onto the FE nodes and then onto the quadrature points using the techniques above.
————————–
At a given time, every element econtains nemarkers. During the FE matrix building process, viscosity
and density values are needed at the quadrature points. One therefore needs to project the values carried
by the markers at these locations. Several approaches are currently in use in the community and the
topic has been investigated by  and  for instance.
elefant adopts a simple approach: viscosity and density are considered to be elemental values, i.e.
all the markers within a given element contribute to assign a unique constant density and viscosity value
to the element by means of an averaging scheme.
While it is common in the literature to treat the so-called arithmetic, geometric and harmonic means
as separate averagings, I hereby wish to introduce the notion of generalised mean, which is a family of
functions for aggregating sets of numbers that include as special cases the arithmetic, geometric and
harmonic means.
If pis a non-zero real number, we can deﬁne the generalised mean (or power mean) with exponent p
of the positive real numbers a1, ... anas:
Mp(a1, ...an) = 1
n
n
X
i=1
ap
i!1/p
(187)
and it is trivial to verify that we then have the special cases:
M−∞ = lim
p→−∞ Mp= min(a1, ...an) (minimum) (188)
M1=n
1
a1+1
a2+··· +1
an
(harm.avrg.) (189)
M0= lim
p0Mp=n
Y
i=1
ai1/n
(geom.avrg.) (190)
M+1 =1
n
n
X
i=1
ai(arithm.avrg.) (191)
M+2 =v
u
u
t1
n
n
X
i=1
a2
i(root mean square) (192)
M+= lim
p+Mp= max(a1, ...an) (maximum) (193)
Note that the proofs of the limit convergence are given in .
An interesting property of the generalised mean is as follows: for two real values pand q, if p < q
then MpMq. This property has for instance been illustrated in Fig. 20 of .
One can then for instance look at the generalised mean of a randomly generated set of 1000 viscosity
values within 1018P a.s and 1023P a.s for 5p5. Results are shown in the ﬁgure hereunder and the
arithmetic, geometric and harmonic values are indicated too. The function Mpassumes an arctangent-
like shape: very low values of p will ultimately yield the minimum viscosity in the array while very high
values will yield its maximum. In between, the transition is smooth and occurs essentially for |p| ≤ 5.
56
1e+18
1e+19
1e+20
1e+21
1e+22
1e+23
-4 -2 0 2 4
M(p)
p
geom.
arithm.
harm.
python codes/ﬁeldstone markers avrg
57
6.9 Matrix (Sparse) storage
The FE matrix is the result of the assembly process of all elemental matrices. Its size can become quite
large when the resolution is being increased (from thousands of lines/columns to tens of millions).
One important property of the matrix is its sparsity. Typically less than 1% of the matrix terms is
not zero and this means that the matrix storage can and should be optimised. Clever storage formats
were designed early on since the amount of RAM memory in computers was the limiting factor 3 or 4
There are several standard formats:
compressed sparse row (CSR) format
compressed sparse column format (CSC)
the Coordinate Format (COO)
Skyline Storage Format
...
I focus on the CSR format in what follows.
6.9.1 2D domain - One degree of freedom per node
Let us consider again the 3 ×2 element grid which counts 12 nodes.
8=======9======10======11
||||
| (3) | (4) | (5) |
||||
4=======5=======6=======7
||||
| (0) | (1) | (2) |
||||
0=======1=======2=======3
In the case there is only a single degree of freedom per node, the assembled FEM matrix will look
like this:
X X X X
X X X X X X
X X X X X X
X X X X
X X X X X X
X X X X X X X X X
X X X X X X X X X
X X X X X X
X X X X
X X X X X X
X X X X X X
X X X X
where the Xstand for non-zero terms. This matrix structure stems from the fact that
node 0 sees nodes 0,1,4,5
node 1 sees nodes 0,1,2,4,5,6
node 2 sees nodes 1,2,3,5,6,7
...
58
node 5 sees nodes 0,1,2,4,5,6,8,9,10
...
node 10 sees nodes 5,6,7,9,10,11
node 11 sees nodes 6,7,10,11
In light thereof, we have
4 corner nodes which have 4 neighbours (counting themselves)
2(nnx-2) nodes which have 6 neighbours
2(nny-2) nodes which have 6 neighbours
(nnx-2)×(nny-2) nodes which have 9 neighbours
In total, the number of non-zero terms in the matrix is then:
NZ = 4 ×4+4×6+2×6+2×9 = 70
In general, we would then have:
NZ = 4 ×4 + [2(nnx 2) + 2(nny 2)] ×6+(nnx 2)(nny 2) ×9
Let us temporarily assume nnx =nny =n. Then the matrix size (total number of unknowns) is
N=n2and
NZ = 16 + 24(n2) + 9(n2)2
A full matrix array would contain N2=n4terms. The ratio of N Z (the actual number of reals to store)
to the full matrix size (the number of reals a full matrix contains) is then
R=16 + 24(n2) + 9(n2)2
n4
It is then obvious that when nis large enough R1/n2.
CSR stores the nonzeros of the matrix row by row, in a single indexed array A of double precision
numbers. Another array COLIND contains the column index of each corresponding entry in the A array.
A third integer array RWPTR contains pointers to the beginning of each row, which an additional pointer
to the ﬁrst index following the nonzeros of the matrix A. A and COLIND have length NZ and RWPTR
has length N+1.
In the case of the here-above matrix, the arrays COLIND and RWPTR will look like:
COLIND = (0,1,4,5,0,1,2,4,5,6,1,2,3,5,6,7, ..., 6,7,10,11)
RW P T R = (0,4,10,16, ...)
6.9.2 2D domain - Two degrees of freedom per node
When there are now two degrees of freedom per node, such as in the case of the Stokes equation in
two-dimensions, the size of the Kmatrix is given by
NfemV=nnpndofV
In the case of the small grid above, we have NfemV=24. Elemental matrices are now 8 ×8 in size.
We still have
4 corner nodes which have 4 ,neighbours,
2(nnx-2) nodes which have 6 neighbours
2(nny-2) nodes which have 6 neighbours
59
(nnx-2)x(nny-2) nodes which have 9 neighbours,
but now each degree of freedom from a node sees the other two degrees of freedom of another node too.
In that case, the number of nonzeros has been multiplied by four and the assembled FEM matrix looks
like:
X X X X X X X X
X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X
X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X X X X X X X
X X X X X X X X X X X X X X X X X X
X X X X X X X X X X X X X X X X X X
X X X X X X X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X
X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X X X X X
X X X X X X X X
X X X X X X X X
Note that the degrees of freedom are organised as follows:
(u0, v0, u1, v1, u2, v2, ...u11, v11)
In general, we would then have:
NZ = 4 [4 ×4 + [2(nnx 2) + 2(nny 2)] ×6+(nnx 2)(nny 2) ×9]
and in the case of the small grid, the number of non-zero terms in the matrix is then:
NZ = 4 [4 ×4+4×6+2×6+2×9] = 280
In the case of the here-above matrix, the arrays COLIND and RWPTR will look like:
COLIND = (0,1,2,3,8,9,10,11,0,1,2,3,8,9,10,11, ...)
RW P T R = (0,8,16,28, ...)
6.9.3 in ﬁeldstone
The majority of the codes have the FE matrix being a full array
a mat = np . z e r os ( ( Nfem , Nfem ) , dtype=np . f l o a t 6 4 )
and it is converted to CSR format on the ﬂy in the solve phase:
s o l = s ps . l i n a l g . s p s o l v e ( s ps . c s r m a t r i x ( a mat ) , r hs )
Note that linked list storages can be used (lil matrix). Substantial memory savings but much longer
compute times.
60
6.10 Mesh generation
Before basis functions can be deﬁned and PDEs can be discretised and solved we must ﬁrst tesselate the
domain with polygons, e.g. triangles and quadrilaterals in 2D, tetrahedra, prisms and hexahedra in 3D.
When the domain is itself simple (e.g. a rectangle, a sphere, ...) the mesh (or grid) can be (more or
less) easily produced and the connectivity array ﬁlled with straightforward algorithms . However,
real life applications can involve etremely complex geometries (e.g. a bridge, a human spine, a car chassis
and body, etc ...) and dedicated algorithms/softwares must be used (see [456, 180, 507]).
We usually distinguish between two broad classes of grids: structured grids (with a regular connec-
tivity) and unstructured grids (with an irregular connectivity).
On the following ﬁgure are shown various triangle- tetrahedron-based meshes. These are obviously
better suited for simulations of complex geometries:
add more examples coming from geo
Let us now focus on the case of a rectangular computational domain of size Lx ×Ly with a regular
mesh composed of nelx×nely=nel quadrilaterals. There are then nnx×nny=nnp grid points. The
elements are of size hx×hy with hx=Lx/nelx.
We have no reason to come up with an irregular/ilogical node numbering so we can number nodes
row by row or column by column as shown on the example hereunder of a 3×2 grid:
8=======9======10======11 2=======5=======8======11
||||||||
| (3) | (4) | (5) | | (1) | (3) | (5) |
||||||||
4=======5=======6=======7 1=======4=======7======10
||||||||
| (0) | (1) | (2) | | (0) | (2) | (4) |
||||||||
0=======1=======2=======3 0=======3=======6=======9
"row by row" "column by column"
The numbering of the elements themselves could be done in a somewhat chaotic way but we follow
the numbering of the nodes for simplicity. The row by row option is the adopted one in ﬁeldstoneand
the coordinates of the points are computed as follows:
x = np . empty (nnp , dtype=np . f l o a t 6 4 )
y = np . empty (nnp , dtype=np . f l o a t 6 4 )
61
co u nt er = 0
f o r jin range ( 0 , nny ) :
f o r ii n r a ng e ( 0 , nnx ) :
x [ c ou nt er ]= i hx
y [ c ou nt er ]= j hy
co u nt er += 1
The inner loop has iranging from 0to nnx-1 ﬁrst for j=0, 1, ... up to nny-1 which indeed corresponds
to the row by row numbering.
We now turn to the connectivity. As mentioned before, this is a structured mesh so that the so-called
connectivity array, named icon in our case, can be ﬁlled easily. For each element we need to store the
node identities of its vertices.