# 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!]

- Introduction
- List of tutorials
- The physical equations
- The heat transport equation - energy conservation equation
- The momentum conservation equations
- The mass conservation equations
- The equations in ASPECT manual
- the Boussinesq approximation: an Incompressible flow
- Stokes equation for elastic medium
- The strain rate tensor in all coordinate systems
- Boundary conditions
- Meaningful physical quantities

- The building blocks of the Finite Element Method
- Solving the flow equations with the FEM
- strong and weak forms
- Which velocity-pressure pair for Stokes?
- Families
- The bi/tri-linear velocity - constant pressure element (Q1P0)
- The bi/tri-quadratic velocity - discontinuous linear pressure element (Q2 P-1)
- The bi/tri-quadratic velocity - bi/tri-linear pressure element (Q2 Q1)
- The stabilised bi/tri-linear velocity - bi/tri-linear pressure element (Q1Q1-stab)
- The MINI triangular element (P1+P1)
- The quadratic velocity - linear pressure triangle (P2P1)
- The Crouzeix-Raviart triangle (P2+P-1)

- Other elements
- The penalty approach for viscous flow
- The mixed FEM for viscous flow
- Solving the elastic equations
- Solving the heat transport equation

- Additional techniques and features
- Picard and Newton
- The SUPG formulation for the energy equation
- Tracking materials and/or interfaces
- Dealing with a free surface
- Convergence criterion for nonlinear iterations
- Static condensation
- The method of manufactured solutions
- Assigning values to quadrature points
- Matrix (Sparse) storage
- Mesh generation
- Visco-Plasticity
- Pressure smoothing
- Pressure scaling
- Pressure normalisation
- The choice of solvers
- The GMRES approach
- The consistent boundary flux (CBF)
- The value of the timestep
- mappings
- Exporting data to vtk format
- Runge-Kutta methods
- Am I in or not?
- Error measurements and convergence rates
- The initial temperature field
- Single layer with imposed heat flux b.c.
- Single layer with imposed heat flux and temperature b.c.
- Kinematic boundary conditions

- fieldstone_01: simple analytical solution
- fieldstone_02: Stokes sphere
- fieldstone_03: Convection in a 2D box
- fieldstone_04: The lid driven cavity
- fieldstone_05: SolCx benchmark
- fieldstone_06: SolKz benchmark
- fieldstone_07: SolVi benchmark
- fieldstone_08: the indentor benchmark
- fieldstone_09: the annulus benchmark
- fieldstone_10: Stokes sphere (3D) - penalty
- fieldstone_11: stokes sphere (3D) - mixed formulation
- fieldstone_12: consistent pressure recovery
- fieldstone_13: the Particle in Cell technique (1) - the effect of averaging
- fieldstone_f14: solving the full saddle point problem
- fieldstone_f15: saddle point problem with Schur complement approach - benchmark
- fieldstone_f16: saddle point problem with Schur complement approach - Stokes sphere
- fieldstone_17: solving the full saddle point problem in 3D
- fieldstone_18: solving the full saddle point problem with Q2Q1 elements
- fieldstone_19: solving the full saddle point problem with Q3Q2 elements
- fieldstone_20: the Busse benchmark
- fieldstone_21: The non-conforming Q1 P0 element
- fieldstone_22: The stabilised Q1 Q1 element
- fieldstone_23: compressible flow (1) - analytical benchmark
- fieldstone_24: compressible flow (2) - convection box
- fieldstone_25: Rayleigh-Taylor instability (1) - instantaneous
- fieldstone_26: Slab detachment benchmark (1) - instantaneous
- fieldstone_27: Consistent Boundary Flux
- fieldstone_28: convection 2D box - Tosi et al, 2015
- fieldstone_29: open boundary conditions
- fieldstone_30: conservative velocity interpolation
- fieldstone_31: conservative velocity interpolation 3D
- fieldstone_32: 2D analytical sol. from stream function
- fieldstone_33: Convection in an annulus
- fieldstone_34: the Cartesian geometry elastic aquarium
- fieldstone_35: 2D analytical sol. in annulus from stream function
- fieldstone_36: the annulus geometry elastic aquarium
- fieldstone_37: marker advection and population control
- fieldstone_38: Critical Rayleigh number
- fieldstone: Gravity: buried sphere
- Problems, to do list and projects for students
- Three-dimensional applications
- Codes in geodynamics
- Matrix properties

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 Whatisaﬁeldstone?....................................... 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×P−1) . 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×P−1) ...................... 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 inﬁeldstone........................................ 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.1Couetteﬂow............................................ 180

36.2SolCx ............................................... 180

36.3Streamlineﬂow.......................................... 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.1Linkingwithourpaper...................................... 193

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) [247].

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

∇

νis the velocity gradient tensor.

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 [296].

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·s−1

ρmass density kg/m3

ηdynamic viscosity Pa·s

λpenalty parameter Pa·s

Ttemperature K

∇gradient operator m−1

∇· divergence operator m−1

ppressure Pa

˙

ε(

ν) strain rate tensor s−1

αthermal expansion coeﬃcient K−1

kthermal conductivity W/(m ·K)

CpHeat capacity J/K

Hintrinsic speciﬁc heat production W/kg

βTisothermal compressibility Pa−1

τ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 [416]:

ρ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

ji∂iuj

= 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 [416].

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)

ρCp∂T

∂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;

•adiabatic compression of material;

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)

ρ0Cp∂T

∂t +v· ∇T−∇·k∇T=ρ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

Cijklkl (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[1]+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ν

1−2ν=ν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

2∂u

∂y +∂v

∂x (15)

˙εzx = ˙εxz =1

2∂u

∂z +∂w

∂x (16)

˙εzy = ˙εyz =1

2∂v

∂z +∂w

∂y (17)

3.7.2 Polar coordinates

˙εrr =∂vr

∂r (18)

˙εθθ =vr

r+1

r

∂vθ

∂θ (19)

˙εθr = ˙εrθ =1

2∂vθ

∂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 = ˙ε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Ω

RΩdΩ1/2

=1

VΩZΩ|

ν|2dΩ1/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

VΩZΩ

v2

rdΩ1/2

(29)

vθ|rms =1

VΩZΩ

v2

θdΩ1/2

(30)

•Pressure (Pa):

•Stress (Pa):

•Strain (X):

•Strain rate (s−1):

•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=k∆T/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 [55]

Nu =

1

LxRLx

0k∂T

∂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 [505] 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

R2−R1= 1, so that R2= 1/(1 −f) and R1=f/(1 −f). Following [276], the Nusselt number at

the inner and outer boundaries are:

Nuinner =fln f

1−f

1

2πZ2π

0∂T

∂r r=R1

dθ (32)

Nuouter =ln f

1−f

1

2πZ2π

0∂T

∂r r=R2

dθ (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

1−f

1

2πZ2π

0∂Tc

∂r r=R1

dθ =fln f

1−f

1

R1

1

ln f= 1 (34)

Nuouter =ln f

1−f

1

2πZ2π

0∂Tc

∂r r=R2

dθ =ln f

1−f

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.K−1.kg−1): 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.m−1.K−1). 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.s−1). 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 α(K−1): 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 '(b−a)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 '(b−a)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 'b−a

n f(a)

2+

n−1

X

k=1

f(a+kb−a

n) + f(b)

2!

where the subintervals have the form [kh, (k+ 1)h], with h= (b−a)/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 n−point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule con-

structed to yield an exact result for polynomials of degree 2n−1 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

b−a(x−a)−1

20

This relationship can be reversed such that when ris known, its equivalent coordinate x∈[a, b] can be

computed:

x=b−a

2(1 + r) + a

From this it follows that

dx =b−a

2dr

and then Zb

a

f(x)dx =b−a

2Z+1

−1

f(r)dr 'b−a

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

2n−1 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

The equivalent Gaussian quadrature writes:

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 [239]:

•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

u∼a1+a2x+a3y+a4xy +a5x2+a6y2

with similar approximations for v, and

p∼b1+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×P−nis as above, except that pressure is approximated via piecewise

discontinuous polynomials of degree n. For instance, P2×P−1is 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

u∼a1+a2x+a3y+a4xy +a5x2+a6y2+a7x2y+a8xy2+a9x2y2

and

p∼b1+b2x+b3y+b4xy

•For these same families, Qm×Q−nis as above, except that the pressure approximation is not

continuous at element boundaries.

•Again for the same families, Qm×P−nindicates 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 [239], 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 [125].

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) = a−b=f1

f(r= +1) = a+b=f2

This leads to

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)

4.4.2 Quadratic basis functions (Q2)

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) = a−b+c=f1

f(r= 0) = a=f2

f(r= +1) = a+b+c=f3

23

This leads to

a=f2b=1

2(−f1 + f3) c=1

2(f1+f3−2f2)

and then replacing a, b, c in Eq. (39) by the above values on gets

f(r) = 1

2r(r−1)f1+ (1 −r2)f2+1

2r(r+ 1)f3

or,

f(r) =

3

X

i=1

Ni(r)fi

with

N1(r) = 1

2r(r−1)

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) = a−b+c−d=f1

f(−1/3) = a−b

3+c

9−d

27 =f2

f(+1/3) = a−b

3+c

9−d

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+ 9f3−f4)

c=9

16 (f1−f2−f3+f4)

Combining the original 4 equations in a diﬀerent way yields

2b+ 2d=f4−f1

2b

3+2d

27 =f3−f2

so that

b=1

16 (f1−27f2+ 27f3−f4)

d=9

16 (−f1+ 3f2−3f3+f4)

Finally,

24

f(r) = a+b+cr2+dr3

=1

16(−1 + r+ 9r2−9r3)f1

+1

16(9 −27r−9r2+ 27r3)f2

+1

16(9 + 27r−9r2−27r3)f3

+1

16(−1−r+ 9r2+ 9r3)f4

=

4

X

i=1

Ni(r)fi

where

N1=1

16(−1 + r+ 9r2−9r3)

N2=1

16(9 −27r−9r2+ 27r3)

N3=1

16(9 + 27r−9r2−27r3)

N4=1

16(−1−r+ 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+ 9r2−9r3)

−1

3(9 −27r−9r2−27r3)

+1

3(9 + 27r−9r2+ 27r3)

+(−1−r+ 9r2+ 9r3)]/16

= [−r+ 9r+ 9r−r]/16 + ...0...

=r(41)

The basis functions derivative are given by

25

∂N1

∂r =1

16(1 + 18r−27r2)

∂N2

∂r =1

16(−27 −18r+ 81r2)

∂N3

∂r =1

16(+27 −18r−81r2)

∂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 + 18r−27r2)

+(−27 −18r+ 81r2)

+(+27 −18r−81r2)

+(−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 + 18r−27r2)

−1

3(−27 −18r+ 81r2)

+1

3(+27 −18r−81r2)

+(−1 + 18r+ 27r2)]

=1

16[−2 + 18 + 54r2−54r2]

= 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 uM→u2.

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(r−1)1

2s(s−1)

N1(r, s) = 1

2r(r+ 1)1

2s(s−1)

N2(r, s) = 1

2r(r+ 1)1

2s(s+ 1)

N3(r, s) = 1

2r(r−1)1

2s(s+ 1)

N4(r, s) = (1 −r2)1

2s(s−1)

N5(r, s) = 1

2r(r+ 1)(1 −s2)

N6(r, s) = (1 −r2)1

2s(s+ 1)

N7(r, s) = 1

2r(r−1)(1 −s2)

N8(r, s) = (1 −r2)(1 −s2)

and their derivatives by:

∂N0

∂r =1

2(2r−1)1

2s(s−1) ∂N0

∂s =1

2r(r−1)1

2(2s−1)

∂N1

∂r =1

2(2r+ 1)1

2s(s−1) ∂N1

∂s =1

2r(r+ 1)1

2(2s−1)

∂N2

∂r =1

2(2r+ 1)1

2s(s+ 1) ∂N2

∂s =1

2r(r+ 1)1

2(2s+ 1)

∂N3

∂r =1

2(2r−1)1

2s(s+ 1) ∂N3

∂s =1

2r(r−1)1

2(2s+ 1)

∂N4

∂r = (−2r)1

2s(s−1) ∂N4

∂s = (1 −r2)1

2(2s−1)

∂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(2r−1)(1 −s2)∂N7

∂s =1

2r(r−1)(−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)(−r−s−1) (58)

N1(r, s) = 1

4(1 + r)(1 −s)(r−s−1) (59)

N2(r, s) = 1

4(1 + r)(1 + s)(r+s−1) (60)

N3(r, s) = 1

4(1 −r)(1 + s)(−r+s−1) (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

of the bilinear quadrilateral element:

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+ 9r2−9r3)/16 N1(t)=(−1 + t+ 9t2−9t3)/16

N2(r) = (+9 −27r−9r2+ 27r3)/16 N2(t) = (+9 −27t−9t2+ 27t3)/16

N3(r) = (+9 + 27r−9r2−27r3)/16 N3(t) = (+9 + 27t−9t2−27t3)/16

N4(r)=(−1−r+ 9r2+ 9r3)/16 N4(t)=(−1−t+ 9t2+ 9t3)/16

31

N01(r, s) = N1(r)N1(s)=(−1 + r+ 9r2−9r3)/16 ∗(−1 + t+ 9s2−9s3)/16

N02(r, s) = N2(r)N1(s) = (+9 −27r−9r2+ 27r3)/16 ∗(−1 + t+ 9s2−9s3)/16

N03(r, s) = N3(r)N1(s) = (+9 + 27r−9r2−27r3)/16 ∗(−1 + t+ 9s2−9s3)/16

N04(r, s) = N4(r)N1(s)=(−1−r+ 9r2+ 9r3)/16 ∗(−1 + t+ 9s2−9s3)/16

N05(r, s) = N1(r)N2(s) = (−1 + r+ 9r2−9r3)/16 ∗(9 −27s−9s2+ 27s3)/16

N06(r, s) = N2(r)N2(s) = (+9 −27r−9r2+ 27r3)/16 ∗(9 −27s−9s2+ 27s3)/16

N07(r, s) = N3(r)N2(s) = (+9 + 27r−9r2−27r3)/16 ∗(9 −27s−9s2+ 27s3)/16

N08(r, s) = N4(r)N2(s)=(−1−r+ 9r2+ 9r3)/16 ∗(9 −27s−9s2+ 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)=1−r−s(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:

A−1=

1 0 0 0 0 0

−3−1 0 4 0 0

−3 0 −1004

220−4 0 0

400−4 4 −4

2 0 2 0 0 −4

In the end, one obtains:

f(r, s) = f1+ (−3f1−f2+ 4f4)r+ (−3f1−f3+ 4f6)s

+(2f1+ 2f2−4f4)r2+ (4f1−4f4+ 4f5−4f6)rs

+(2f1+ 2f3−4f6)s2

=

6

X

i=1

Ni(r, s)fi(77)

with

N1(r, s)=1−3r−3s+ 2r2+ 4rs + 2s2

N2(r, s) = −r+ 2r2

N3(r, s) = −s+ 2s2

N4(r, s)=4r−4r2−4rs

N5(r, s)=4rs

N6(r, s)=4s−4rs −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 −r−s)(1/3−r−s)(2/3−r−s) (78)

N1(r, s) = 9

2r(r−1/3)(r−2/3) (79)

N2(r, s) = 9

2s(s−1/3)(s−2/3) (80)

N3(r, s) = 27

2(1 −r−s)r(2/3−r−s) (81)

N4(r, s) = 27

2(1 −r−s)r(r−1/3) (82)

N5(r, s) = 27

2rs(r−1/3) (83)

N6(r, s) = 27

2rs(r−2/3) (84)

N7(r, s) = 27

2(1 −r−s)s(s−1/3) (85)

N8(r, s) = 27

2(1 −r−s)s(2/3−r−s) (86)

N9(r, s) = 27rs(1 −r−s) (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=f2−f1c2=f3−f1c3=f4−f1

f(r, s, t) = c0+c1r+c2s+c3t

=f1+ (f2−f1)r+ (f3−f1)s+ (f4−f1)t

=f1(1 −r−s−t) + f2r+f3s+f4t

=X

i

Ni(r, s, t)fi

Finally,

34

N1(r, s, t)=1−r−s−t

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(r−1) 0.5s(s−1) 0.5t(t−1)

N2= 0.5r(r+ 1) 0.5s(s−1) 0.5t(t−1)

N3= 0.5r(r+ 1) 0.5s(s+ 1) 0.5t(t−1)

N4= 0.5r(r−1) 0.5s(s+ 1) 0.5t(t−1)

N5= 0.5r(r−1) 0.5s(s−1) 0.5t(t+ 1)

N6= 0.5r(r+ 1) 0.5s(s−1) 0.5t(t+ 1)

N7= 0.5r(r+ 1) 0.5s(s+ 1) 0.5t(t+ 1)

N8= 0.5r(r−1) 0.5s(s+ 1) 0.5t(t+ 1)

N9= (1.−r2) 0.5s(s−1) 0.5t(t−1)

N10 = 0.5r(r+ 1) (1 −s2) 0.5t(t−1)

N11 = (1.−r2) 0.5s(s+ 1) 0.5t(t−1)

N12 = 0.5r(r−1) (1 −s2) 0.5t(t−1)

N13 = (1.−r2) 0.5s(s−1) 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(r−1) (1 −s2) 0.5t(t+ 1)

N17 = 0.5r(r−1) 0.5s(s−1) (1 −t2)

N18 = 0.5r(r+ 1) 0.5s(s−1) (1 −t2)

N19 = 0.5r(r+ 1) 0.5s(s+ 1) (1 −t2)

N20 = 0.5r(r−1) 0.5s(s+ 1) (1 −t2)

N21 = (1 −r2) (1 −s2) 0.5t(t−1)

N22 = (1 −r2) 0.5s(s−1) (1 −t2)

N23 = 0.5r(r+ 1) (1 −s2) (1 −t2)

N24 = (1 −r2) 0.5s(s+ 1) (1 −t2)

N25 = 0.5r(r−1) (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) [142]. 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 [278] section 3.6.2

5.3.1 The bi/tri-linear velocity - constant pressure element (Q1×P0)

discussed in example 3.71 of [278]

5.3.2 The bi/tri-quadratic velocity - discontinuous linear pressure element (Q2×P−1)

5.3.3 The bi/tri-quadratic velocity - bi/tri-linear pressure element (Q2×Q1)

This element, implemented in penalised form, is discussed in [45] and the follow-up paper [46]. CHECK

Biquadratic velocities, bilinear pressure. See Hood and Taylor. The element satisﬁes the inf-sup

condition [259]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 [278]

5.3.6 The quadratic velocity - linear pressure triangle (P2×P1)

From [417]. 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×P−1)

From [130]: 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 [417]. 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 [278]

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 [129], [398] or [243].

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)

[330] 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 [43] 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 [225] .

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 0∂N2

∂x 0∂N3

∂x 0∂N4

∂x 0

0∂N1

∂y 0∂N2

∂y 0∂N3

∂y 0∂N4

∂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:

ZΩe

Ni

∇ · σdΩ + ZΩe

NibdΩ=0

We can integrate by parts and drop the surface term2:

ZΩe

∇Ni·σdΩ = ZΩe

NibdΩ

or,

ZΩe

∂Ni

∂x 0∂Ni

∂y

0∂Ni

∂y

∂Ni

∂x

·

σxx

σyy

σxy

dΩ = ZΩe

NibdΩ

Let i= 1,2,3,4 and stack the resulting four equations on top of one another.

ZΩe

∂N1

∂x 0∂N1

∂y

0∂N1

∂y

∂N1

∂x

·

σxx

σyy

σxy

dΩ = ZΩe

N1bx

bydΩ (94)

ZΩe

∂N2

∂x 0∂N2

∂y

0∂N2

∂y

∂N2

∂x

·

σxx

σyy

σxy

dΩ = ZΩe

Nibx

bydΩ (95)

ZΩe

∂N3

∂x 0∂N3

∂y

0∂N3

∂y

∂N3

∂x

·

σxx

σyy

σxy

dΩ = ZΩe

N3bx

bydΩ (96)

ZΩe

∂N4

∂x 0∂N4

∂y

0∂N4

∂y

∂N4

∂x

·

σxx

σyy

σxy

dΩ = ZΩe

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

ZΩe

BT·

σxx

σyy

σxy

dΩ = ZΩe

NbdΩ

and ﬁnally: ZΩe

BT·[λK+ηC]·B·VdΩ = ZΩe

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:

ZΩe

BT·[λK+ηC]·BdΩ

| {z }

Ael(8×8)

·V

|{z}

(8x1)

=ZΩe

NbdΩ

| {z }

Bel(8×1)

or,

ZΩe

λBT·K·BdΩ

| {z }

Aλ

el(8×8)

+ZΩe

ηBT·C·BdΩ

| {z }

Aη

el(8×8)

·V

|{z}

(8x1)

=ZΩe

NbdΩ

| {z }

Bel(8×1)

INTEGRATION - MAPPING

reduced integration [260]

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=A−1·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:

ZΩe

Nν

i

∇ · σdΩ + ZΩe

Nν

i

b dΩ =

0 (109)

ZΩe

Np

i

∇ ·vdΩ = 0 (110)

We can integrate by parts and drop the surface term3:

ZΩe

∇Nν

i·σdΩ = ZΩe

Nν

i

bdΩ (111)

ZΩe

Np

i

∇ ·vdΩ = 0 (112)

or,

ZΩe

∂Nν

i

∂x 0 0 ∂Nν

i

∂y

∂Nν

i

∂z 0

0∂Nν

i

∂y 0∂Nν

i

∂x 0∂Nν

i

∂z

0 0 ∂Nν

i

∂z 0∂Nν

i

∂x

∂Nν

i

∂y

·

σxx

σyy

σzz

σxy

σxz

σyz

dΩ = ZΩe

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:

ZΩe

BT·

σxx

σyy

σzz

σxy

σxz

σyz

dΩ = ZΩe

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

0∂Nν

1

∂y 0··· 0∂Nν

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 0∂Nν

1

∂x ··· ∂Nν

mv

∂z 0∂Nν

mv

∂x

0∂Nν

1

∂z

∂Nν

1

∂y ··· 0∂Nν

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 ZΩe

BT·[−Np·

P+C·B·

V]dΩ = ZΩe

NbdΩ (126)

or, −ZΩe

BT·NpdΩ

| {z }

G

·

P+ZΩe

BT·C·BdΩ

| {z }

K

·

V=ZΩe

NbdΩ

| {z }

~

f

(127)

where the matrix Kis of size (mv∗ndofv×mv∗ndofv), and matrix Gis of size (mv∗ndofv×mp∗ndofp).

43

Turning now to the mass conservation equation:

0 = ZΩe

Np

∇ ·v dΩ

=ZΩe

Np

mv

X

i=1 ∂Nν

i

∂x ui+∂N ν

i

∂y vi+∂N ν

i

∂z widΩ

=ZΩe

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Ω

=ZΩe

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Ω

=ZΩe

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)

say something about minus sign?

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 (ndofv∗Nv)×(ndofv∗Nv) 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 (ndofv∗Nv)×(ndofp∗Np) 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 (ndofp∗Np)×(ndofv∗Nv). The size of the global FE matrix is N=ndofv∗Nv+ndofp∗NpNote

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: ML−2T−2. 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] = ML−2T−2L3=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 [239], 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

[239]: if for instance Np

iis piecewise-constant on element e(of value 1), the elemental weak form of the

mass conservervation equation is

ZΩe

Np

i

∇ ·

ν=ZΩe

∇ ·

ν=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

4−2−2000

−2 4 −2000

−2−2 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 x−y

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

2∂u

∂y +∂v

∂x 0

1

2∂u

∂y +∂v

∂x ∂v

∂y 0

0 0 0

(138)

˙

εd=1

3

2∂u

∂x −∂v

∂y

1

2∂u

∂y +∂v

∂x 0

1

2∂u

∂y +∂v

∂x −∂u

∂x + 2∂v

∂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 0∂Nν

2

∂x 0∂Nν

3

∂x 0. . . ∂Nν

mv

∂x 0

0∂Nν

1

∂y 0∂Nν

2

∂y 0∂Nν

3

∂y . . . 0∂Nν

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

4−2 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):

ρcp∂T

∂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ρcp∂T

∂t +

ν·

∇TdΩ = 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(u∂T

∂x +v∂T

∂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 u∂Nθ

j

∂x +v∂N θ

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θ

1∂xNθ

2. . . ∂xNθ

mT

∂yNθ

1∂yNθ

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(k−1)

δ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(k−1)

other time discr schemes!

49

6 Additional techniques and features

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 saddlepoint

python codes/ﬁeldstone saddlepoint q2q1

python codes/ﬁeldstone saddlepoint q3q2

python codes/ﬁeldstone burstedde

6.7.1 Analytical benchmark I

Taken from [144]. 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+ 72y2−48y3+ 12)x2

+(−2 + 24y−72y2+ 48y3)x+ 1 −4y+ 12y2−8y3

by= (8 −48y+ 48y2)x3+ (−12 + 72y−72y2)x2

+(4 −24y+ 48y2−48y3+ 24y4)x−12y2+ 24y3−12y4

With this prescribed body force, the exact solution is

u(x, y) = x2(1 −x)2(2y−6y2+ 4y3)

v(x, y) = −y2(1 −y)2(2x−6x2+ 4x3)

p(x, y) = x(1 −x)−1/6

Note that the pressure obeys RΩp dΩ = 0 One can turn to the spatial derivatives of the ﬁelds:

∂u

∂x = (2x−6x2+ 4x3)(2y−6y2+ 4y3) (158)

∂v

∂y =−(2x−6x2+ 4x3)(2y−6y2+ 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 [144], 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 [141]. It is for a unit square with ν=µ/ρ = 1 and the smooth exact solution is

u(x, y) = x+x2−2xy +x3−3xy2+x2y(162)

v(x, y) = −y−2xy +y2−3x2y+y3−xy2(163)

p(x, y) = xy +x+y+x3y2−4/3 (164)

52

Note that the pressure obeys RΩp dΩ=0

bx=−(1 + y−3x2y2) (165)

by=−(1 −3x−2x3y) (166)

6.7.3 Analytical benchmark III

This benchmark begins by postulating a polynomial solution to the 3D Stokes equation [141]:

v=

x+x2+xy +x3y

y+xy +y2+x2y2

−2z−3xz −3yz −5x2yz

(167)

and

p=xyz +x3y3z−5/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 =−2−3x−3y−5x2y(171)

˙εxy =1

2(x+y+ 2xy2+x3) (172)

˙εxz =1

2(−3z−10xyz) (173)

˙εyz =1

2(−3z−5x2z) (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)

[96] 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 [45]. The two-dimensional domain is a unit square. The body forces are:

fx= 128[x2(x−1)212(2y−1) + 2(y−1)(2y−1)y(12x2−12x+ 2)]

fy= 128[y2(y−1)212(2x−1) + 2(x−1)(2x−1)y(12y2−12y+ 2)]

(183)

The solution is

u=−256x2(x−1)2y(y−1)(2y−1)

v= 256x2(y−1)2x(x−1)(2x−1)

p= 0 (184)

Another choice:

fx= 128[x2(x−1)212(2y−1) + 2(y−1)(2y−1)y(12x2−12x+ 2)] + y−1/2

fy= 128[y2(y−1)212(2x−1) + 2(x−1)(2x−1)y(12y2−12y+ 2)] + x−1/2

(185)

The solution is

u=−256x2(x−1)2y(y−1)(2y−1)

v= 256x2(y−1)2x(x−1)(2x−1)

p= (x−1/2)(y−1/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

very large contrasts inside a single element and perturb the quadrature. Please read section 3.3

of [253] 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(r−1) N2(r)=1−r2N3(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 [139] and [149] 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)

M−1=n

1

a1+1

a2+··· +1

an

(harm.avrg.) (189)

M0= lim

p→0Mp=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 [83].

An interesting property of the generalised mean is as follows: for two real values pand q, if p < q

then Mp≤Mq. This property has for instance been illustrated in Fig. 20 of [414].

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 −5≤p≤5. 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

decades ago. [407]

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(n−2) + 9(n−2)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(n−2) + 9(n−2)2

n4

It is then obvious that when nis large enough R∼1/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=nnp∗ndofV

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 [452]. 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.