Introduction

The finite element (FE) method is a computational technique for obtaining approximate solutions to the partial differential equations that arise in scientific and engineering applications and is used widely in geodynamic modelling (see Christensen, 1984, 1992; Baumgardner, 1985; Naimark and Malevsky, 1988; King et al., 1990; Naimark and Ismail- Zadeh, 1995; Moresi and Solomatov, 1995; Moresi et al., 2003; Ismail-Zadeh et al., 1998, 2001a,b, 2004a,b, 2006, 2007). Introduced in the middle of the twentieth century (Hrennikoff, 1941; McHenry, 1943; Courant, 1943) the FE method has emerged as one of the most powerful numerical methods so far devised. Rather than approximating the partial differential equation directly as with finite difference methods (see Chapter 2), the FE method utilises a variational problem that involves an integral of the differential equation over the model domain. This domain is divided into a number of sub-domains called finite elements, and the solution of the partial differential equation is approximated by a simple polynomial function on each element. These polynomials have to be pieced together so that the approximate solution has an appropriate degree of smoothness over the entire domain. Once this has been done, the variational integral is evaluated as a sum of contributions from each finite element. The result is a set of algebraic equations for the approximate solution having a finite size rather than the original infinite-dimensional partial differential equation.