1 Introduction and Brief History
The IFISS software [
16] was developed by Elman, Ramage, and Silvester [
8]. It can be run in MATLAB (developed by the MathWorks
\({}^\copyright\) ) or Gnu Octave (free software). It is structured as a stand-alone package for studying discretisation algorithms for
partial differential equations (PDEs), and for exploring and developing algorithms in numerical linear and nonlinear algebra for solving the associated discrete systems. It can be used as a pedagogical tool for studying these issues or more elementary ones such as the properties of Krylov subspace iterative methods. Investigative numerical experiments in a teaching setting enable students to develop deduction and interpretation skills and are especially useful in helping students to remember critical ideas in the long term. IFISS is also an established starting point for developing code for specialised research applications (as evidenced by the variety of citations to it; see Reference [
17]) and is extensively used by researchers in numerical linear algebra as a source of reproducible test matrices of arbitrarily large dimension.
The development of the MATLAB functionality during the period 1990–2005 opened up the possibility of creating a problem-based-learning environment (notably, the IFISS package) that could be used together with standard teaching mechanisms to facilitate understanding of abstract theoretical concepts. The functionality of IFISS was significantly extended in the period between 2005 and 2015—culminating in the publication of the review article cited in Reference [
9], which coincided with the publication of the second edition of the monograph [
10].
A unique feature of IFISS is its comprehensive nature. For each problem it addresses, it enables the study of both discretisation and iterative solution algorithms, as well as the interaction between the two and the resulting effect on solution cost. However, it is restricted to the solution of PDEs on two-dimensional spatial domains. This limitation can be overcome by adding the new IFISS3D toolbox [
13] to the existing IFISS software. The three-dimensional finite element approximation and error estimation strategies included in the new software are specified in the next section. Section
3 describes three reference problems that provide a convenient starting point for studying rates of convergence of the approximations to the true solution. The structure of the IFISS3D package is discussed in Section
4. The directory structure is intended to simplify the task of extending the functionality to other PDE problems and higher-order finite element methods. Case studies of two important aspects of three-dimensional finite element approximation are presented in Section
5.
2 Discretisation and Error Estimation Specifics
The IFISS3D software generates approximations to the solution of PDEs modelling physical problems in three spatial dimensions. The starting point for the process is a finite element partitioning of a domain of interest
\(D\subset \mathbb {R}^3\) into
\(n_e\) hexahedral (brick) elements
\({\square }_e\subset D\) ,
\(e=1,2,3,\) \(\ldots ,n_e\) so
where the upper bar represents the closure of the union.
An arbitrary element
\({\square }_{e}\) is a hexahedron with six faces and with local vertex coordinates
\((x^e_i,y^e_i,z^e_i)\) ,
\(i=1,2, \ldots , 8\) ordered as shown in Figure
1.
The simplest choice of a conforming finite element space in
\(\mathbb {R}^3\) is the
\(\mathbb {Q}_1\) approximation space of piecewise
trilinear polynomials on each element
\({\square }_{e}\) . The continuity of the global approximation is ensured by defining a Lagrangian basis for
\(q^e_1\) at the eight vertices of the hexahedron, that is,
Additional geometric flexibility (
stretched grids) can be incorporated by constructing an
isoparametric transformation from the reference cube
\([-1,1]^3\) (denoted
\({\square }_{\star }\) ) to a general hexahedron
\({\square }_{e}\) . To this end, we define the following basis functions for the reference element:
where
\(\xi _i, \eta _i, \zeta _i\) are node values of
\({\square }_{\star }\) and
\(\xi ,\eta ,\zeta \in [-1,1]\) and map to an arbitrary hexahedral element
\({\square }_{e}\) with vertices
\((x^e_i,y^e_i,z^e_i)\) ,
\(i=1,2,\ldots , 8\) by the change of variables
This mapping is illustrated in Figure
2.
The IFISS3D software also provides the option of higher-order approximation using globally continuous piecewise
triquadratic polynomials (
\(\mathbb {Q}_2\) ). The continuity of the global
\(\mathbb {Q}_2\) approximation is ensured by defining a Lagrangian basis for
\(q_2^e\) at the 8 vertices of the hexahedron together with the 19 additional nodes shown in Figure
3.
The
\(\mathbb {Q}_2\) isoparametric transformation is given by
with reference basis functions
where
Finite element approximation of a linear elliptic PDE problem invariably results in a sparse system of algebraic equations (the so-called
Galerkin system). A typical Galerkin system matrix is assembled from element matrices that are associated with integrals of products of mapped derivatives defined on the reference cube (see Reference [
12]). These integrals are computed
exactly in IFISS3D using tensor-product Gauss rules of appropriate degree (see Reference [
10, p. 30]). A key feature of the IFISS code design is the inner-column indexing and data structure (two-dimensional element matrices are stored as a three-dimensional array). This indexing ensures that all element matrix calculations and subsequent assembly can be efficiently vectorised and multi-threaded using the BLAS functionality that is built into MATLAB.
Another fundamental feature of the IFISS software is the use of hierarchical error estimation. This strategy was developed for scalar elliptic PDEs by Bank & Smith [
2], but has been extended to more general PDE problems (including systems of PDEs) over the past two decades. Crucially, the hierarchical approach yields reliable estimates of the
error reduction that can be expected using an enhanced approximation. It also provides a rigorous setting for establishing the convergence of adaptive refinement strategies, such as those that are built into the T-IFISS package [
3], and the ML-SGFEM software [
6] associated with the work of Crowder et al. [
7] on multilevel stochastic Galerkin finite element methods for parametric PDEs.
A posteriori error estimation in practical finite element software (such as deal.II [
1], DUNE [
4], or FEniCS [
12]) is typically done using residual error estimation strategies. This requires the computation of norms of PDE residuals in the interior of each element and norms of flux jumps (edge residuals) on inter-element faces. The additional computational cost of hierachical error estimation is nontrivial but not overwhelming in practice. This is discussed further in Section
5. Having generated a solution using
\(\mathbb {Q}_1\) approximation,
1 computed interior and edge residuals are input as source data for
element PDE problems that are solved numerically using an enhanced approximation space. In IFISS3D, one can construct the enhanced space using triquadratic basis functions on the original element
\((\mathbb {Q}_2(h))\) or trilinear basis functions defined on a subdivision of the original element into eight smaller ones
\((\mathbb {Q}_1(h/2))\) . These basis or “bubble” functions are associated with the white nodes illustrated in Figure
4(a), leading to linear algebra systems of dimension 19. Alternatively,
reduced versions of these spaces of dimension 7, denoted
\(\mathbb {Q}^r_2(h)\) and
\(\mathbb {Q}^r_1(h/2)\) , can also be constructed by incorporating only the basis functions associated with the interior node and the central nodes on each face, as illustrated in Figure
4(b). In all four cases, a low-dimensional system must be solved for every element in the mesh. This calculation is efficiently vectorised in IFISS.
4 Structure of the Software Package
IFISS is designed for the MATLAB coding environment. This means that the source code is readable, portable, and easy to modify. All local calculations (quadrature in generating element matrices, application of essential boundary conditions, a posteriori error estimation) are vectorised over elements—thus the code runs efficiently on contemporary Intel processor architectures. IFISS3D has been developed for MATLAB (post 2018b) and tested with the current release (7.2) of Gnu Octave. The main directory is called diffusion3D, and this needs to be added as a subdirectory of the main IFISS directory. The subdirectories of diffusion3D are organised as follows:
\(\bullet\) /grids/
This directory contains all the functions associated with domain discretisation. Three types of domain are included in the first release. Introducing a new domain type is straightforward. A new function needs to be included that saves nodal information (arrays xyz, bound3D) and (triquadratic) element information (mv, mbound3D) in an appropriately named datafile. This file will be subsequently read by an appropriate driver function associated with the specific PDE being solved.
\(\bullet\) /graphs/
This directory contains the functions associated with the visualisation of the computed solution (nodal data) and the estimated errors (element data). The tensor-product subdivision structure simplifies the code structure substantially—plotting can be efficiently done using the built-in
slice functionality. Similarly, solution data defined on a one-dimensional incision into the domain of interest can be plotted using the function
xyzsectionplot. An illustration is shown in Figure
8.
\(\bullet\) /approximation/
This directory contains all the functions associated with setting up the discrete matrix system associated with the PDE of interest. The functions
femq1_diff3D and
femq2_diff3D set up the stiffness and mass matrices associated with the problems discussed in Section
3. Essential boundary conditions are imposed by a subsequent call to the function
nonzerobc3D. Extending the functionality by combining components of IFISS with IFISS3D to cover (a) nonisotropic diffusion and (b) Stokes flow problems (using
\(\mathbb {Q}_2\) –
\(\mathbb {Q}_1\) mixed approximation) is a straightforward exercise. Efficient approximation of the solution of the heat equation on a three-dimensional domain can also be done with ease: either using the adaptive time-stepping functionality built into IFISS or using one of the ODE integrators built into MATLAB. The functions associated with
a posteriori error estimation can be found in four separate subdirectories associated with the four options described in Section
2.
\(\bullet\) /solvers/
The MATLAB sparse direct solver (
\(\backslash\) ) has far from optimal complexity in a three-dimensional setting. This is explored in a case study in the next section.
Algebraic multigrid (AMG) functionality is included in this directory to enable exploration of an optimal solution strategy. If one does not have access to an efficient AMG setup routine, then the linear solver that is recommended when the dimension of the system exceeds
\(10^5\) is
MINRES (Minimum Residual) iteration preconditioned by an incomplete Cholesky factorisation of the system matrix with zero fill-in.
3 This strategy is encoded in the
it_solve3D function with a residual stopping tolerance of
\(10^{-10}\) . Solving the system in Example
3 using this strategy gives a solution in 66 iterations. The associated CPU time is half a second. This is over 10 times faster than the corresponding backslash solve!
\(\bullet\) /test_problems/
This directory contains all the high-level driver functions such as
diff3D_testproblem (the main driver). It also contains the functions associated with the problem data (right-hand side function and essential boundary specifications). The structure makes it straightforward to solve (
1) together with nonzero boundary data
\(u= g_D\) on
\(\partial D\) .
Help for IFISS is integrated into the MATLAB help facility. The command helpme_diff3D gives information on solving a Poisson problem in three dimensions. Starting from the main IFISS directory, typing help diffusion3D/ <subdirectory name> gives a complete list of the files in that subdirectory. Using MATLAB, the function names are “clickable” to give additional information. The initial release of IFISS3D comprises \(\sim\) 70 individual functions and script files. Simply type help <file-name> for further information on any of these.