Georgios Papanikos
David J. Silvester
Abstract
IFISS is an established MATLAB finite element software package for studying strategies for solving partial differential equations (PDEs). IFISS3D is a new add-on toolbox that extends IFISS capabilities for elliptic PDEs from two to three space dimensions. The open-source MATLAB framework provides a computational laboratory for experimentation and exploration of finite element approximation and error estimation, as well as iterative solvers. The package is designed to be useful as a teaching tool for instructors and students who want to learn about state-of-the-art finite element methodology. It will also be useful for researchers as a source of reproducible test matrices of arbitrarily large dimension.
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 \(\begin{equation*} \nonumber \nonumber {\left\lbrace \begin{array}{ll} {\square }_e \ \text{is open in} \ \mathbb {R}^3 \\ D = \overline{\cup _{e=1}^{n_e}{\square }_e} \\ {\square }_{i}\cap {\square }_{j} = \emptyset , \ \ i\ne j, \end{array}\right.} \end{equation*}\) 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.
Fig. 1. Hexahedral (brick) element (a) vertex numbering and (b) face numbering.
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, \(\begin{equation*} \phi ^e_i(\boldsymbol {x}) = {\left\lbrace \begin{array}{ll} 1 & \text{if} \ \ x = x^e_i, y = y^e_i, \ z= z^e_i, \\ 0 & \text{at the other vertices} \end{array}\right.}, \quad i =1,2,\ldots , 8 . \end{equation*}\)
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: \(\begin{eqnarray*} \ell ^i(\xi ,\eta ,\zeta) = 1/8 (1+\xi _i\xi)&(1+\eta _i\eta)(1+\zeta _i\zeta), \quad i=1,2,\ldots ,8, \nonumber \nonumber \end{eqnarray*}\) 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 \(\begin{equation*} x(\xi ,\eta ,\zeta) = \sum _{i=1}^8 x^e_i \ell ^i(\xi ,\eta ,\zeta), \ y(\xi ,\eta ,\zeta) = \sum _{i=1}^8 y^e_i \ell ^i(\xi ,\eta ,\zeta), \ z(\xi ,\eta ,\zeta) = \sum _{i=1}^8 z^e_i \ell ^i(\xi ,\eta ,\zeta). \end{equation*}\) This mapping is illustrated in Figure 2.
Fig. 2. Isoparametric mapping from a reference cube.
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.
Fig. 3. (a) \(\mathbb {Q}_1\) nodes and numbering and (b) \(\mathbb {Q}_2\) nodes and numbering.
The \(\mathbb {Q}_2\) isoparametric transformation is given by \(\begin{equation*} x(\xi ,\eta ,\zeta) = \sum _{i=1}^{27} x^e_i {\psi }^i(\xi ,\eta ,\zeta), \quad y(\xi ,\eta ,\zeta) = \sum _{i=1}^{27 }y^e_i {\psi }^i(\xi ,\eta ,\zeta), \quad z(\xi ,\eta ,\zeta) = \sum _{i=1}^{27} z^e_i {\psi }^i(\xi ,\eta ,\zeta), \end{equation*}\) with reference basis functions \(\begin{align*} {\psi }^i(\xi ,\eta ,\zeta) = {N}^k(\xi) \, {N}^l(\eta) \, {N}^m(\zeta),\ i=1,2,\ldots ,27; \ k,l,m =1,2,3, \end{align*}\) where \(\begin{equation*} {N}^1(\xi) = 1/2(\xi -1)\xi , \ {N}^2(\xi) = 1-\xi ^2, \ \ {N}^3(\xi) = 1/2(1+\xi)\xi . \end{equation*}\)
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.
Fig. 4. White nodes associated with bubble functions used to construct (a) full \(\mathbb {Q}_2(h)\) or \(\mathbb {Q}_1(h/2)\) and (b) reduced \(\mathbb {Q}^r_2(h)\) and \(\mathbb {Q}^r_1(h/2)\) error estimation spaces.
3 REFERENCE PROBLEMS
Three test problems are built into the IFISS3D toolbox. Illustrative results for these problems are discussed in this section. The reported timings were obtained on a 2.9 GHz six-Core Intel Core i9 MacBook using the
A \(\mathbb {Q}_1\) finite element solution to (1) and (2) defined on the domain \(D=[-1,1]^3\) is shown in Figure 5
Fig. 5. Solution and estimated energy error distribution for test problem 1.
Exploiting Galerkin orthogonality, the exact energy error can be estimated by comparing the energy of a reference solution2 with the energy of the computed finite element solution \(\Vert u_{\rm ref}-u_h \Vert _E^2 = \Vert u_{\rm ref}\Vert _E^2 - \Vert u_h \Vert _E^2.\) Estimated errors that are computed using this strategy are presented in Table 1.
Table 1. Estimated Energy Errors for Convex Domain Test Problem
We observe that the \(\mathbb {Q}_1\) energy errors are reducing by a factor of 2 with each grid refinement. This is consistent with the optimal \(O(h)\) rate of convergence predicted theoretically in the case of a \(H^2\)-regular problem. The \(\mathbb {Q}_2\) energy errors are reducing more rapidly with grid refinement. The observed rate is slightly less than \(O(h^2)\), which is the expected rate when solving a \(H^3\)-regular problem using triquadratic approximation.
A \(\mathbb {Q}_1\) finite element solution to the Poisson problem (1) and (2) defined on the domain \(D = {[-1,1]^3\setminus [-1,0)\times [-1,0) \times [-1,1] }\) is shown in Figure 6
Fig. 6. Solution and estimated energy error distribution for test problem 2.
Estimated errors for the second test problem are presented in Table 2. In this case, we observe that the \(\mathbb {Q}_1\) and \(\mathbb {Q}_2\) energy errors are both reducing by a factor of less than 2 with each grid refinement. This is exactly what one would expect—the edge singularity limits the rate of convergence that is possible using uniform grids. The second notable feature is that the \(\mathbb {Q}_2\) energy error is a factor of 2 smaller than the \(\mathbb {Q}_1\) error for the same number of degrees of freedom (the results on the same horizontal line). This behaviour is also consistent with expectations (see Schwab [15]).
Table 2. Estimated Energy Errors for the Staircase Domain Test Problem
A \(\mathbb {Q}_1\) finite element solution to the Poisson problem (1) and (2) defined on the cut domain \(D={[-1,1]^3\setminus (-\epsilon ,\epsilon)\times [0,1]\times (-\epsilon ,\epsilon)}\) with \(\epsilon =0.01\) is shown in Figure 7
Fig. 7. Solution and estimated energy error distribution for test problem 3.
In this computation, the borehole domain has been subdivided into a tensor-product grid of \(48\times 32 \times 48\) elements with geometric stretching in the x and z direction to capture the geometry without using an excessive number of elements. The grid spacing increases from \(h_{\min }=0.01=\epsilon\) within the hole to \(h_{\max }=0.0625\) next to the boundary, so the maximum element aspect ratio (adjacent to the hole) is 6.25. The dimension of the resulting linear algebra system is 85,833, and the MATLAB R2021b sparse direct solver solves the system in about 6 seconds. As anticipated, the error in the approximation is concentrated in a small region in the neighbourhood of the borehole, making this a very challenging problem to solve efficiently.
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
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
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
Fig. 8. Incisions through the borehole test problem solution visualised in Figure 7. The incisions show the solution in the x direction in the plane \(z=0\) for three different height values y.
This directory contains all the functions associated with setting up the discrete matrix system associated with the PDE of interest. The functions
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
This directory contains all the high-level driver functions such as
Help for IFISS is integrated into the MATLAB help facility. The command
5 CASE STUDIES
Two important aspects of three-dimensional finite element approximation that can be investigated easily in IFISS3D are discussed in this section.
5.1 Effectivity of a Posteriori Error Estimation Strategies
The effectiveness of hierarchical error estimation is well established in a two-dimensional setting (see, for example, Elman et al. [10, Table 1.4]). The IFISS3D software offers a choice of four such error estimation strategies in conjunction with \(\mathbb {Q}_1\) approximation. Computed error estimates obtained when solving the first test problem discussed in Section 3 are presented in Table 3. The four estimates are associated with the white nodes shown in Figure 4. The estimated energy errors should be compared with the reference energy errors listed in Table 1.
Table 3. Computed Error Estimates \(\eta _{\bullet }\) for the First Test Problem Using Four Different Hierarchical Strategies
Table 4 lists the associated effectivity indices. The indices get closer to 1 as the mesh is refined when the \(\mathbb {Q}_2(h)\) and \(\mathbb {Q}^r_2(h)\) strategies are employed. However, the effectivity indices for the \(\mathbb {Q}_1(h/2)\) and \({\mathbb {Q}^r_1(h/2)}\) strategies stagnate around 0.87 and 0.77, respectively. All four error estimates are correctly reducing by a factor of 2 with each grid refinement. In light of these results, the relatively cheap \({\mathbb {Q}^r_2(h)}\) is set to be the default option in IFISS3D. Extensive testing on other problems indicates that this estimator consistently underestimates the error by a small amount.
Table 4. Effectivity Indices \(\theta _\bullet := \eta _\bullet / \Vert {u_{\rm ref}} -u_h \Vert _E\) for Test Problem 1
All the errors reported in Table 3 were computed after making a boundary element correction. This is a postprocessing step wherein the local problems associated with elements that have one or more boundary faces are modified so the (zero) error on the boundary is enforced as an essential boundary condition. The motivation for making this correction is to recover the property of asymptotic exactness in special cases.4 The correction is, however, difficult to vectorise efficiently, raising the question as to whether it is worth including in a three-dimensional setting.
Computed effectivity indices for the special case of solving the Poisson problem
(3) \(\begin{align} -\nabla ^2 u &= f \quad \text{in } D =[-1,1]^3 \end{align}\)
(4) \(\begin{align} u &= 0 \quad \text{on } \partial D \end{align}\) with the right-hand side function f chosen so the exact solution is the triquadratic function
(5) \(\begin{align} u(\boldsymbol {x}) &= (1- x^2)(1-y^2)(1-z^2), \end{align}\) are presented in Table 5. The second and fourth columns are the results computed after making the boundary correction. The asymptotic exactness of the \(\mathbb {Q}_2(h)\) strategy can be clearly seen. The third and fifth columns list the results that are computed when the boundary correction is not made. Comparing results with the second and fourth columns, it is evident that the boundary correction reduces the estimated error and, more importantly, that the size of correction tends to zero in the limit \(h\rightarrow 0\). The \({\mathbb {Q}^r_2(h)}\) strategy is not asymptotically exact, so to speed up the computation the default setting in IFISS3D is to simply neglect the boundary correction. Thus, in the case of the finest grid in Table 5 (over 2 million elements) the \(\eta ^*_{\mathbb {Q}^r_2(h)}\) error estimate is computed in less than 9 seconds. This is significantly less than the time taken to compute the finite element solution itself (the
Table 5. Effectivity Indices \(\theta _\bullet := \eta _\bullet / { \Vert u-u_h \Vert _E}\) for the Poisson Problem (3)–(5)
5.2 Fast Linear Algebra
The solution of the (Galerkin) linear system is the computational bottleneck when solving a Poisson problem in three dimensions. To illustrate this point, a representative timing comparison of the distinct solution components when solving the first test problem using \(\mathbb {Q}_1\) approximation with default settings is presented in Table 6. The system assembly includes the grid generation. The overall time is the elapsed time from start to finish and includes the time taken to visualise the solution and the associated error. What is immediately apparent is the fact that the system assembly times and the error estimation times scale approximately linearly with the number of elements (or equivalently, the dimension n of the system matrix). The backslash solve times, in contrast, grow like the square of the number of elements on the most refined grids. The memory requirement for the sparse factors of the system matrix also increases at a much faster rate than \(O(n)\). Iteration counts and representative solution times for solving the linear systems in Table 6 using preconditioned MINRES are presented in Table 7.
Table 6. Representative Component Timings (in Seconds) when Solving Test Problem 1
Table 7. MINRES Iteration Counts and Timings (in Seconds) when Solving Test Problem 1 with Incomplete Cholesky Preconditioning and a Residual Stopping Tolerance of \(10^{-10}\)
The optimal \(O(n)\) complexity of the overall solution algorithm can be recovered by solving the linear system using a short-term Krylov subspace iteration such as MINRES in combination with an algebraic multigrid (AMG) preconditioning strategy (see Reference [10, sec. 2.5.3]). The set-up phase of AMG is a recursive procedure: Heuristics associated with algebraic relations (“strength of connections”) between the unknowns are used to generate a sequence of progressively coarser representations \(A_{\ell }\), \(\ell =1,2,\ldots , L\) of the Galerkin system matrix \(A := A_1\). The solution (preconditioning) phase approximates the action of the inverse of A on a vector by cycling through the associated grid sequence. At each level, a fixed-point iteration (typically point Gauss-Seidel) is applied to “smooth” the residual error that is generated by interpolation or restriction of the error vector generated at the previous level. If coarsening is sufficiently rapid, then the work associated with the preconditioning step will be proportional to the number of unknowns in the original system matrix.
The algorithmic complexity of any AMG coarsening strategy can be characterised by a few parameters. First, the grid complexity is defined as \(\begin{equation*} c_G := \frac{1}{n_{1}} \sum _{\ell =1}^L n_{\ell }, \end{equation*}\) where \(n_{\ell }\) is the dimension of the coarse grid matrix \(A_{\ell }\) at level \(\ell\). Starting from a uniform grid of \(\mathbb {Q}_{1}\) elements, if full coarsening (in each spatial direction) is done at each level, then \(n_{\ell }\) reduces by a factor of 8 at each level, in which case, we obtain \(c_{G} \approx {8}/{7}\). A value of \(c_{G}\) higher than this suggests that coarsening has not been done isotropically. The operator complexity is typically defined by \(\begin{equation*} c_A := \frac{1}{{\rm nnz}(A_{1})} \sum _{\ell =1}^{L} {\rm nnz}(A_{\ell }), \end{equation*}\) where \({\rm nnz}(A)\) is the number of nonzeros in the matrix. This parameter provides information about the associated storage requirements for the coarse grid matrices generated. If uniform coarsening is done and the coarse grid matrices correspond to the usual finite element discretisation on those grids, then we would expect \(c_{A} \approx c_{G}\). In practice, however, the coarse grid matrices become progressively denser, with larger stencil sizes, as the level number increases. If the matrices become too dense, then this may cause an issue with the computational cost of applying the chosen smoother. To quantify this, the average stencil size \(\begin{equation*} c_S: = \frac{1}{L} \sum _{\ell =1}^L \frac{{\rm nnz} (A_{\ell })}{n_{\ell }}, \end{equation*}\) should be compared with the average stencil size at the finest level, that is, \(\begin{equation*} c_1 := \frac{{\rm nnz}(A_{1})}{n_{1}}. \end{equation*}\)
An implementation of the coarsening strategy developed by Ruge and Stüben [14] is included in the original IFISS software. The corresponding IFISS3D function
AMG complexities and timings obtained when solving test problem 1 are presented in Table 8. The total times reported for the \(\mathbb {Q}_1\) approximation should be compared with the corresponding backslash solve timings recorded in Table 6 and Table 7. The setup times were recorded using the interface to the HSL_MI20 code and scale close to linearly with the problem dimension. This behaviour is consistent with the results that were reported for the same problem and discretisation by Boyle et al. [5, Ex. 4.5.1]. Looking at the AMG grid data in Table 8, we note that the grid complexity is under control and stays close to the optimal value using either of the two approximation strategies. The operator complexity results are also encouraging—increasing slowly as the dimension of the problem is increased. The growth in the average stencil size is not unexpected. The value of \(c_S\) is within a factor of 3 of \(c_1\) when using \(\mathbb {Q}_2\) approximation on the finest grid. The difference between the highlighted total time and the associated setup time is the time taken by preconditioned MINRES to reach the residual stopping tolerance of \(10^{-10}\) with the default smoothing parameters (that is, two pre- and post-smoothing sweeps using point Gauss-Seidel).
Table 8. AMG Grid Complexity Data and Representative Linear Solver Timings for Test Problem 1
6 SUMMARY AND FUTURE DEVELOPMENTS
The IFISS3D toolbox extends the capabilities of IFISS [16] to facilitate the numerical solution of elliptic PDEs on three-dimensional spatial domains that can be partitioned into hexahedra. In particular, it allows users to investigate the convergence properties of trilinear (\(\mathbb {Q}_{1}\)) and triquadratic (\(\mathbb {Q}_{2}\)) finite element approximation for test problems whose solutions have varying levels of spatial regularity and the performance of a range of iterative solution algorithms for the associated discrete systems, including an optimal AMG solver. For \(\mathbb {Q}_{1}\) approximation, the effectivity of four distinct state-of-the-art hierarchical error estimation schemes can also be explored. The IFISS3D software is structured in such a way that, when integrated into the existing IFISS software, users can easily solve a range of other PDE problems, including time-dependent ones, using \(\mathbb {Q}_{1}\) and \(\mathbb {Q}_{2}\) elements on three-dimensional spatial domains. IFISS together with IFISS3D is intended to be useful as a teaching tool and can be used to produce matrices of arbitrarily large dimension for testing linear algebra algorithms. Future developments of IFISS3D will be documented on GitHub.
Footnotes
1 Hierarchical error estimation for \(\mathbb {Q}_2\) approximation is not included in the current release of IFISS3D.
Footnote2 \(\Vert u_{\rm ref}\Vert ^2_E \, {:= \int _D \nabla u_{\rm ref}\cdot \nabla u_{\rm ref}} \, = 0.64539192\) computed using \(\mathbb {Q}_2\) approximation on a grid of 64\(^3\) elements.
Footnote3 The incomplete factorisation function
Footnoteichol provided in MATLAB R2021b is highly optimised.4 An estimator is said to be asymptotically exact if the effectivity of the estimator tends to 1 when \(h\rightarrow 0\).
Footnote5 The HSL_MI20 source code and associated MATLAB interface is freely available to staff and students of recognised educational institutions. The inclusion of the compiled
Footnotemex file is prohibited by the terms of the HSL academic licence.
- [1] . 2021. The deal.II finite element library: Design, features, and insights. Comput. Math. Appl. 81 (2021), 407–422.
DOI: Google Scholar
Cross Ref
- [2] . 1993. A posteriori error estimates based on hierarchical bases. SIAM J. Numer. Anal. 30, 4 (1993), 921–935. Retrieved from https://www.jstor.org/stable/2158183.Google ScholarDigital Library
- [3] . 2021. T-IFISS: A toolbox for adaptive FEM computation. Comput. Math. Appl. 81 (2021), 373–390.
DOI: Google ScholarCross Ref - [4] . 2016. The distributed and unified numerics environment, version 2.4. Arch. Num. Soft. 4, 100 (2016), 13–29.
DOI: Google ScholarCross Ref - [5] . 2010. HSL_MI20: An efficient AMG preconditioner for finite element problems in 3D. Int. J. Numer. Meth. Eng. 82 (2010), 64–98.Google ScholarCross Ref
- [6] . 2022. ML-SGFEM (Multilevel Stochastic Galerkin Finite Element Method) Software, Version 1.0. Retrieved from https://github.com/ceapowell/ML_SGFEM.Google Scholar
- [7] . 2019. Efficient adaptive multilevel stochastic Galerkin approximation using implicit a posteriori error estimation. SIAM J. Sci. Comput. 41, 3 (2019), A1681–A1705.Google ScholarDigital Library
- [8] . 2007. Algorithm 866: IFISS, a Matlab toolbox for modelling incompressible flow. ACM Trans. Math. Softw. 33 (2007), 2–14. Retrieved from Google ScholarDigital Library
- [9] . 2014. IFISS: A computational laboratory for investigating incompressible flow problems. SIAM Rev. 56, 2 (2014), 261–273. Retrieved from Google ScholarDigital Library
- [10] . 2014. Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics (2nd ed.). Oxford University Press, Oxford, UK.Google ScholarCross Ref
- [11] HSL 2015. HSL Mathematical Software Library. Retrieved from https://www.hsl.rl.ac.uk/catalogue/hsl_mi20.html.Google Scholar
- [12] . 2012. Automated Solution of Differential Equations by the Finite Element Method. Springer.
DOI: Google ScholarCross Ref - [13] . 2022. Incompressible Flow and Iterative Solver 3D Software (IFISS3D), version 1.0. Retrieved from https://github.com/mcbssds/IFISS3D.Google Scholar
- [14] . 1987. Algebraic multigrid. In Multigrid Methods, (Ed.)
, (Frontiers in Applied Mathematics, Vol. 3). SIAM, Philadelphia, PA, 73–130.Google ScholarCross Ref - [15] . 1998. p- and hp- Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics. Oxford University Press, Oxford, UK.Google Scholar
- [16] . 2019. Incompressible Flow and Iterative Solver Software (IFISS), version 3.6. Retrieved from http://www.manchester.ac.uk/ifiss/.Google Scholar
- [17] swMATH 2022. An information service for mathematical software. Retrieved from https://swmath.org/software/4398.Google Scholar
Index Terms
- IFISS3D: A Computational Laboratory for Investigating Finite Element Approximation in Three Dimensions
Recommendations
A Family of Finite Element Stokes Complexes in Three Dimensions
We construct finite element Stokes complexes on tetrahedral meshes. In the lowest-order case, the finite elements in the complex have 4, 18, 16, and 1 degrees of freedom on each tetrahedron, respectively. As a consequence, we obtain $\grad\curl$-...
Algebraic multigrid for higher-order finite elements
Two related approaches for solving linear systems that arise from a higher-order finite element discretization of elliptic partial differential equations are described. The first approach explores direct application of an algebraic-based multigrid ...
Finite Element Approximation of a Three Dimensional Phase Field Model for Void Electromigration
We consider a finite element approximation of a phase field model for the evolution of voids by surface diffusion in an electrically conducting solid. The phase field equations are given by the nonlinear degenerate parabolic system $$\gamma\frac{\...
Comments