A second-order face-centred finite volume method on general meshes with automatic mesh adaptation

A second-order face-centred finite volume strategy on general meshes is proposed. The method uses a mixed formulation in which a constant approximation of the unknown is computed on the faces of the mesh. Such information is then used to solve a set of problems, independent cell-by-cell, to retrieve the local values of the solution and its gradient. The main novelty of this approach is the introduction of a new basis function, utilised for the linear approximation of the primal variable in each cell. Contrary to the commonly used nodal basis, the proposed basis is suitable for computations on general meshes, including meshes with different cell types. The resulting approach provides second-order accuracy for the solution and first-order for its gradient, without the need of reconstruction procedures, is robust in the incompressible limit and insensitive to cell distortion and stretching. The second-order accuracy of the solution is exploited to devise an automatic mesh adaptivity strategy. An efficient error indicator is obtained from the computation of one extra local problem, independent cell-by-cell, and is used to drive mesh adaptivity. Numerical examples illustrating the approximation properties of the method and of the mesh adaptivity procedure are presented. The potential of the proposed method with automatic mesh adaptation is demonstrated in the context of microfluidics.


Introduction
Finite volume (FV) methods are one of the most popular computational methods for solving systems of conservation laws. [1][2][3][4] These methods are usually classified into two families, namely cell-centred FVs and vertex-centred FVs depending on the definition of the unknowns at the centroids or at the vertices of the cells respectively. The main attractive properties of FV methods are their numerical efficiency, local conservation and robustness which make them appealing solutions to treat flow problems of industrial interest. [5][6][7] However, one of the drawbacks of low-order cell-centred and vertex-centred FVs is the need for a reconstruction of the gradient. In this context, the quality of the reconstruction is directly linked to the quality of the mesh, leading to an important loss of accuracy, and even second-order convergence, in the presence of highly distorted or stretched cells. 8, 9

Problem statement and mixed formulation
An open bounded domain Ω ∈ R nsd is considered, where n sd denotes the number of spatial dimensions. The boundary of the domain is split into the non-overlapping Dirichlet boundary, Γ D , where the solution is known, and the Neumann boundary, Γ N , where the normal flux is known.
The computational domain is assumed to be partitioned in n e non-overlapping cells Ω e , for e = 1, . . . , n e . The boundary of each cell is expressed as the union of a set of faces (edges in two dimensions), Γ e,j , for j = 1, . . . , n e f , where n e f denotes the number of faces (edges in two dimensions) of the cell Ω e . The FCFV method considers the strong form of the Poisson equation written in mixed form, via the introduction of the variable q, and in a cell-by-cell fashion, namely in Ω e , and for e = 1, . . . , n e , ∇ · q = s in Ω e , and for e = 1, . . . , n e , u = u D on ∂Ω e ∩ Γ D , n · q = −t on ∂Ω e ∩ Γ N , un = 0 on Γ, n · q = 0 on Γ. (1) where s is a source term, n is the outward unit normal to the boundary, u D is the known value of the solution on the Dirichlet boundary, t is the known value of the flux on the Neumann boundary and Γ, defined by is the so-called internal mesh skeleton.
It is worth noting that the last two equations in (1) impose the continuity of the solution and the normal flux, respectively, across the internal faces of the mesh, the jump operator being defined as that is, the sum of the quantity inside the left and right cell, Ω l and Ω r respectively, sharing a face of the mesh skeleton. 25

Strong form of the local and global problems
Following the standard rationale of HDG 15,[26][27][28][29] and FCFV 10, 20, 21 methods, the strong mixed form (1) is split into the so-called local and global problems. The local problems are defined independently in each cell as          q e + ∇u e = 0 in Ω e , ∇ · q e = s in Ω e , u e = u D on ∂Ω e ∩ Γ D , u e =û on ∂Ω e \ Γ D , (4) for e = 1, . . . , n e . It is worth noting that each local problem contains only Dirichlet boundary conditions and it introduces a new independent variable,û, called the hybrid variable, that corresponds to the solution at the cell faces (edges in two dimensions).
The global problem is defined on the mesh skeleton and the Neumann boundary as      un = 0 on Γ, n · q = 0 on Γ, n · q = −t on Γ N . (5) The first equation in (5) can be henceforth omitted because the continuity of the solution is automatically satisfied due to the imposition of the Dirichlet boundary condition in the local problems and to the uniqueness of the hybrid variable on Γ.

Weak form of the local and global problems
The recently proposed second-order FCFV 21  for all test functions v ∈ V 1 (Ω e ) and for e = 1, . . . , n e . Here, V 1 (Ω e ) is the space of at most linear functions in Ω e and V 0 (Ω e ) is the space of constant functions in Ω e . It is worth mentioning that the weak formulation (6a) is obtained by selecting an arbitrary test function in [V 0 (Ω e )] nsd .
The numerical flux, q h e , introduced in (6b) is defined as n e · q h e := n e · q h e + τ e (P 0 u h e − u D ) on ∂Ω e ∩ Γ D , n e · q h e + τ e (P 0 u h e −û h ) elsewhere, where τ e > 0 is the so-called stabilisation parameter 15,[26][27][28][29] and, similar to, 22, 23 the projection operator P 0 over the space of constant functions is introduced.
Remark 1. As discussed in, 21 the two ingredients required to obtain a second-order FCFV method are the use of piecewise linear functions to approximate the primal variable and the introduction of the projection operator in the definition of the numerical flux.
The weak formulation of the local problems is obtained after introducing the definition of the numerical flux into (6b), performing an integration by parts of the first term and exploiting that ∇ · q h e = 0, being q h e a piecewise constant function in each cell Ω e . Hence, it reads: for all v ∈ V 1 (Ω e ) and for e = 1, . . . , n e .
The weak formulation of the global problem is derived by following a similar procedure. It reads: find where an arbitrary test function has been selected fromV 0 (Γ ∪ Γ N ). By using the expression of the numerical flux introduced in (7), the discrete weak formulation becomes:

FCFV discretisation
Introduce the notation D e for the set of faces of cell Ω e on the Dirichlet boundary Γ D and B e for the faces of Ω e in Γ ∪ Γ N , that is the ones not on the Dirichlet boundary. Henceforth, the stabilisation parameter, Dirichlet and Neumann data are considered to assume constant values τ j , u D,j and t j respectively on each face/edge Γ e,j of the cell Ω e .
The discrete local problem obtained from the weak form (8) provides an explicit expression of the primal and mixed variables in each cell as a function of the hybrid variable on the cell faces/edges, namely where q e denotes the value of the mixed variable at the centroid of the cell and u e denotes the nodal values of the primal variable. The right hand side vectors depending on the problem data are defined as and the remaining matrices and vectors in the discrete equations are given by In the above expressions, n en and n fn are the numbers of nodes in each cell and face respectively and δ Ij is the Kronecker delta. The vector p e,j , defined as (p e,j ) I := 1 is used to compute the projection of the primal variable from the space of linear to the one of constant polynomial functions. In addition, F e,k denotes the set of nodes of the cell Ω e that belong to the face Γ e,k and the indicator function of a set is defined as The discrete global problem obtained by plugging the explicit expressions (11) in the weak form (10) results in a global system of equations where the unknown vector corresponds to the hybrid variable at the cell faces/edges. It can be written as where the global matrix K and vectorf are the result of assembling the contribution from each cell, given by for i, j ∈ B e and with δ ij denoting the Kronecker delta.

New basis functions for the second-order FCFV on general meshes
The second-order FCFV method proposed in, 21 summarised in the previous section, employs nodal shape functions to define the approximation of the primal variable, namely where {N e J } nen J=1 is the set of linear Lagrange polynomials in the cell Ω e and u e J , for J = 1, . . . , n en , are the corresponding nodal values of the unknown function.
This section shows that this approach is only applicable with simplicial (triangular and tetrahedral) cells and proposes a new approximation space for general polygons and polyhedrons in two and three dimensions respectively. The proposed formulation is first presented for the Poisson equation introduced in section 2, and is extended to Stokes equations in section 3.4.

The second-order FCFV with nodal basis functions
Consider the Poisson model problem (1) in two dimensions. The matrix m e used in the local problem to write the solution in the cell as a function of the solution on the faces can thus be computed analytically. For a triangular cell, and assuming a constant value of the stabilisation parameter τ e for all the faces, the matrix is given by This matrix is invertible, with determinant equal to (τ 3 e /16)|Γ e,1 ||Γ e,2 ||Γ e,3 |. In a similar fashion, assume a constant value of the stabilisation parameter τ e for all the faces of a quadrilateral cell. Using a nodal approximation with bilinear Lagrange polynomials, the matrix is given by Contrary to the matrix obtained in (20) for a triangular cell, the matrix m e for a quadrilateral cell is singular and the local problem in (11) cannot be solved. Hence, in order to devise a second-order FCFV method suitable to handle quadrilateral cells in 2D and hexahedral, prismatic and pyramidal cells in 3D, an alternative description of the primal variable needs to be considered.

New linear basis functions for the second-order FCFV method
This work proposes the use of a linear approximation of the primal variable in each cell, irrespective of its shape. More precisely, the approximation of the primal variable is defined as where the number of terms of the expansion is selected as M = n sd + 1, { N e J } M J=1 is a set of basis functions that span the space of polynomials of, at most, degree one and c e J , for J = 1, . . . , M are coefficients appropriately defined to describe the unknown function.
The following basis functions are proposed in this work

Second-order FCFV discretisation of the Poisson equation
Considering the approximation of the primal variable proposed in the previous section and a constant approximation of the mixed and hybrid variables, the discrete local problem obtained from the weak form (8) is It is worth noting that equation (24a), expressing the mixed variable in terms of the hybrid variable, is identical to the discrete equation (11a) of the original second-order FCFV, whereas equation (24b) requires the following definitionsb The vectorp e,j , given by is introduced to compute the projection of the primal linear variable on the space of the constant functions over a face Γ e,j , wherex e,j = (x e,j 1 , . . . ,x e,j nsd ) denotes the centroid of the face Γ e,j . Analogously, the discrete global problem obtained from the weak form (10) using the expressions of primal and mixed variable in (24), results in a global system of equations where the unknown vector corresponds to the hybrid variable at the cell faces/edges. It can be written as where the global matrix K and vectorf are the result of assembling the contribution from each cell, given by for i, j ∈ B e .

Second-order FCFV discretisation of the Stokes equations
In this section, the procedure described above for the derivation of the second-order FCFV formulation of the Poisson problem is extended to Stokes equations. By introducing the mixed variable L, the strong form of the mixed problem is written cell-by-cell as in Ω e and for e = 1, . . . , n e , in Ω e and for e = 1, . . . , n e , ∇ · u = 0 in Ω e and for e = 1, . . . , n e , u = u D on ∂Ω e ∩ Γ D , with the solvability constraint for the uniqueness of pressure given by where ρ e is the mean value of the pressure on the boundary of the cell Ω e . 30,31 As for the scalar problem, the last two equations in (31) impose the continuity of the velocity and the normal flux, respectively, across the internal faces of the mesh.
Following the rationale in, 10, 21 an additional unknown, the hybrid velocity u, is introduced on the cell faces and the condition u = u is enforced on the boundary ∂Ω e \ Γ D for e = 1, . . . , n e . The weak formulation of the Stokes problem in each cell is: for all test functions w ∈ [V 1 (Ω e )] nsd , where the definition of the numerical normal flux featuring the projection operator P 0 is utilised In equation (33), [V 1 (Ω e )] nsd is the space of n sd dimensional vectors whose components are at most linear functions in Ω e , whereas [V 0 (Ω e )] nsd×nsd and V 0 (Ω e ) are the spaces of constant n sd × n sd tensorial and scalar functions in the cell Ω e , respectively. It is worth mentioning that the weak formulations in equations (33a) and (33c) are obtained by selecting an arbitrary constant test function in the spaces [V 0 (Ω e )] nsd×nsd and V 0 (Ω e ), respectively.
The weak formulation of the global problem is derived in a similar way by imposing Neumann and transmission conditions on Γ N and Γ, respectively, and a compatibility condition for the incompressiblity constraint in each cell Ω e , for e = 1, . . . , n e . It reads: where an arbitrary test function in the space [V 0 (Γ ∪ Γ N )] nsd has been selected in equation (35a).
Remark 3. Equation (33c) of the local problem coincides with the compatibility condition (35b) enforcing the divergence-free nature of the velocity field cell-by-cell. Following the strategy discussed for the first and second-order FCFV method, 10, 21 this equation is omitted from the local problems and imposed solely in the global one since it involves only the global variable u h .
The discrete FCFV local problems are obtained starting from the linear discretisation described in section 3.2 for the velocity u h e , a constant cell-by-cell approximation of the pressure p h e , the mixed variable L h e and the mean pressure ρ h e and a constant face-by-face approximation of the hybrid velocity u h . It follows that where It is worth noticing that equations (36a) and (36c) are identical to the original second-order FCFV 21 and only equation (36b) is affected by the change of basis discussed above. More precisely, The discrete FCFV global system is obtained by plugging the expressions (36) of u e , p e and L e into equation (35), leading to It is straightforward to observe that the matrix of problem (40) features a saddle-point structure with a symmetric block Kûû in the top-left position, as classical in the approximation of incompressible Stokes equations. Both the left and right hand sides of the linear system above are assembled by computing for i, j ∈ B e the contributions of each cell as The operator utilised to project the linear velocity in each cell on the space of constant functions on the face Γ e,i is defined via the matrix P e,i P e,i : in 2D and in 3D,p e,i being the vector introduced in (28) for the projection in the scalar case.

Computational aspects
The implementation of a FCFV method features three steps. First, a preprocess routine to compute all the elemental quantities required by the method. Second, the computation, assembly and solution of a global problem on the mesh faces with the hybrid variable as unknown. Finally, the solution of n e local problems to retrieve the primal and mixed variable in each cell.
The current implementation assumes that the stabilisation parameter is constant in the whole domain. Therefore, for the Poisson problem, the entries of the matricesm e in equation (26) and the vectorsd j and r j in equation (27) can be precomputed. Analogously, for the Stokes problem, the entries of the matrices M e and R j in equation (38) and (39) respectively and the vector D j in equation (39) can be precomputed. These terms only require the calculation of the integral of the new basis functions over a generic face of a cell. More precisely, for a generic face in two dimensions and for a triangular face in three dimensions, the integrals are given by For a generic quadrilateral face in three dimensions, analytical integration is only feasible for the first basis function, N 1 , whereas for the remaining basis functions, N k for k = 2, . . . M , a numerical quadrature is employed.
Similarly, the entries of the vectorsf e and F e in equations (26) and (38), respectively for the Poisson and Stokes problems, can be precomputed via an integral of the new basis functions over a generic cell. For a triangular or tetrahedral cell, such integral is given by whereas for other cell types, the integral is computed numerically. Hence, the computational cost of the preprocess routine to determine the terms in (26)- (27) and (38)- (39) is limited.
The proposed second-order FCFV method requires the solution of a global system of equations with exactly the same number of unknowns and non-zero entries of the global matrix of the first-order FCFV proposed in. 10,20 This is because the second-order FCFV uses the same constant approximation space for the hybrid variable, which is the only variable featuring in the global systems of equations (30) and (41), respectively for the Poisson and Stokes problems.
The extra cost of the second-order FCFV compared to the first-order FCFV is due to the extra operations required to assemble the global system. Tables 1 and 2 detail the number of operations required to compute one elemental matrix and one elemental right hand side vector for the Poisson and Stokes problems respectively, for the proposed second-order FCFV method and the first-order FCFV method introduced in. 10 More precisely, table 1 presents the cost of computing one elemental contribution of equation (30) and one of equation 22 in, 10 whereas table 2 compares the cost of building one elemental block of equation (41) with one of equation 39 in. 10 Although the second-order method requires more operations for a given spatial discretisation, the extra accuracy provided results in a more efficient method, when the computational cost required to achieve a given accuracy is considered. To quantify this gain, Figure 1 shows the relative error of the velocity field measured in the L 2 (Ω) norm as a function of the CPU time for the first and second-order methods using different cell types. The test case considered involves the solution of the Stokes equations, in two and three dimensions, for a problem described in the next section, where the analytical solution is known. The results show that to achieve a 1% error, the second-order method is almost one order of magnitude Method   Triangle Quadrilateral Tetrahedron Hexahedron Prism Pyramid   First-order  126  212  252  534  380  380  Second-order  354  592  932  1,962 1,400 1,400    faster in two dimensions whereas in three dimensions the second-order method is more than four order of magnitude faster.
Finally, the solution of equation (24) and (36) for the local Poisson and Stokes problems respectively, relies on precomputed quantities and can be easily performed in parallel being such computation independent cell-by-cell.

An automatic mesh adaptivity strategy for the second-order FCFV
In, 21 the authors devised an error indicator using the higher convergence rate of the second-order FCFV method with respect to the original first-order FCFV approach. 10,20 The main drawback of the proposed error indicator is the high computational cost that induces the solution of an extra global problem to estimate the error of a numerical solution.
This work proposes a new error indicator that does not require the solution of an extra global problem. Instead, only a local computation cell-by-cell is required to devise an accurate and efficient indicator to drive mesh adaptivity. The error indicator proposed here is significantly cheaper because it exploits the solution of the global system already computed for the second-order FCFV method to solve an extra local problem cell-by-cell. Hence, its cost is negligible when compared to the cost of assembling and solving an extra global problem as discussed in. 21 Hereafter, the proposed strategy is described for the Poisson problem but it is also applicable to the Stokes equations. After the global problem given by equation (29) is solved, two approximations of the primal variable are computed using the same hybrid variableû. On the one hand, a first approximation of the primal variable, u e , is computed by solving the local, cell-by-cell, problem of equation (24). On the other hand, a second approximation u e is obtained by solving an extra local, cell-by-cell, problem corresponding to the first-order FCFV, first presented in, 10 namely where α e := j∈Ae |Γ e,j |τ j , β e := |Ω e |s e + j∈De |Γ e,j |τ j u D,j , and A e denotes the set of all faces of cell Ω e .
The local error indicator for the cell Ω e is thus defined as and it is employed to devise an automatic mesh adaptivity process. First, the following a priori local error estimate for elliptic problems is recalled 32-35 where u ex and u h are the exact solution and its constant approximation in the cell Ω e respectively, h e is the characteristic cell size and C is an unknown constant. Then, using a classical Richardson extrapolation, the desired cell size is calculated as where ε is the user defined target error in each cell. It is worth noting that the corresponding formula in, 21 namely equation (47), presents a typo in the exponent of ε/E e .

Numerical studies
This section presents an extensive set of experiments to numerically validate the optimal convergence properties of the proposed method and to show its robustness against the choice of the stabilisation parameter and the mesh properties such as cell distortion and stretching. The examples presented involve meshes of different cell types, namely triangular and quadrilateral cells in two dimensions, and tetrahedral, hexahedral, prismatic and pyramidal cells in three dimensions. In addition, results with hybrid meshes are presented for the first time in the context of the FCFV method.
The model problem considered for the Poisson equation involves the numerical solution of (1) in Ω = [0, 1] nsd . In two dimensions, the source term and boundary data are selected such that the analytical solution is known and given by u ex (x 1 , x 2 ) = exp α sin(ax 1 + cx 2 ) + β cos(bx 1 + dx 2 ) ,  are set on the rest of the boundary. For the three dimensional Poisson problem, the source term and boundary data are selected such that the analytical solution is The domain Ω = [0, 1] nsd is also utilised for the Stokes equations (31) with viscosity ν = 1. For the two dimensional case, the source term and boundary conditions are devised in order for the analytical velocity and pressure fields to be On Γ N = {(x 1 , x 2 ) ∈ R 2 | x 2 = 0}, a Neumann condition representing a pseudo-traction is imposed, whereas the analytical velocity enforcing Dirichlet conditions is set on the rest of the boundary. Similarly, in three dimensions, Neumann boundary conditions are imposed on and Dirichlet conditions on the rest of the boundary, to match the analytical expressions of velocity and pressure given by To perform the mesh convergence study in two dimensions, a set of eight uniform triangular and quadrilateral meshes are generated. The meshes corresponding to the third level of refinement are displayed in Figures 2 (a) and (b). To demonstrate the flexibility of the proposed method, hybrid meshes made of quadrilateral and triangular cells are also considered. The third hybrid mesh is displayed in Figure 2   In three dimensions, a set of six uniform tetrahedral, hexahedral, prismatic and pyramidal meshes are generated. In addition, hybrid meshes containing a mixture of different cell types are also considered. Figure 3 shows the fourth level of mesh refinement for the different types of meshes considered and table 5 shows the statistics for the hybrid meshes.  Table 5: Details of the six hybrid meshes for the convergence study in three dimensions.

Optimal convergence of the FCFV method
The first experiment involves a mesh convergence study for the Poisson and Stokes problems both in two and three dimensions and using meshes with different cell types. In all cases the stabilisation parameter is selected to be τ = 10 4 in two dimensions and τ = 10 2 for three dimensional problems. A detailed study of the influence of the stabilisation parameter on the accuracy of the proposed FCFV method is provided in section 5.2.
For the two dimensional Poisson problem, Figure 4 shows the relative error, measured in the L 2 (Ω) norm, of the primal and mixed variables as a function of the characteristic cell size. The results show optimal, quadratic, convergence of the primal variable for triangular, quadrilateral and hybrid meshes with almost identical accuracy in the three cases. For the mixed variable an optimal, linear, convergence is also observed, again, with almost identical accuracy in the three cases.
The same mesh convergence study is performed for the three dimensional Poisson problem. Figure 5 shows the relative error of the primal and mixed variables, measured in the L 2 (Ω) norm, as a function of the characteristic cell size. The results display a similar qualitative behaviour when compared to the two dimensional case. An optimal, quadratic, convergence is observed for the primal variable and an optimal, linear, rate of convergence is observed for the mixed variable, for all cell types and hybrid meshes.   Tetrahedral and hybrid meshes provide slightly more accurate results when compared to hexahedral, prismatic and pyramidal meshes.
Next, the mesh convergence study is performed for the Stokes problem in two and three dimensions. Figure 6 displays the relative error of velocity, pressure and gradient of velocity, measured in the L 2 (Ω) norm, as a function of the characteristic cell size. The results show again an optimal quadratic convergence of the error of the velocity for all the different types of meshes. For the pressure and the gradient of the velocity, optimal, linear, convergence of the error is also observed for all types of meshes. The same conclusions are observed from the results of the Stokes problem in three dimensions, displayed in Figure 7.

Influence of the stabilisation parameter
The influence of the stabilisation parameter τ is studied numerically. The results in this section only consider the Stokes problem as further numerical examples, not reported here for brevity, have shown that identical conclusions are obtained for the Poisson problem.   In all cases, the results show that a low value of the stabilisation parameter leads to a high error for the primal variable, whereas a large value, namely τ = 10 4 in two dimensions and τ = 10 2 in three dimensions, provides the maximum accuracy for the velocity. When triangular or tetrahedral cells are considered, the error of the gradient of the velocity and the pressure is found to be independent on the value of the stabilisation parameter used. In contrast, for quadrilateral and hexahedral cells, the error of the gradient of the velocity and the pressure decreases as the value of τ increases. For these cell types, the maximum accuracy is reached for a value of τ = 10 4 in two dimensions and τ = 10 2 in three dimensions. Finally, for prismatic and pyramidal cells, the accuracy of the gradient of the velocity and the pressure is less sensitive to the choice of τ and the qualitative behaviour is extremely similar to the one observed for quadrilateral and hexahedral cells.
Henceforth, the stabilisation parameter is selected as τ = 10 4 in two dimensions and τ = 10 2 in three dimensions, for both Poisson and Stokes problems and for any cell type.

Influence of cell distortion and stretching
Previous experiments, employed to test the optimal approximation properties of the proposed method, involved regular meshes. In this section, the effect of cell distortion and stretching on the accuracy of the proposed method is studied. This is of major importance for the method to be applicable to more complicated problems involving complex geometries and for its extension to computational fluid dynamics applications involving boundary layers.
Cell distortion is introduced by perturbing the internal nodes of the mesh according to a random variation of maximum magnitude h min /4, where h min is the minimum edge of the regular mesh. For cells with quadrilateral faces, the motion is constrained to ensure that all faces on the distorted mesh are planar. 37 Figure 9    The mesh convergence results for the Stokes problem in two and three dimensions using meshes with distorted cells are shown in Figures 10 and 11, respectively. In two dimensions, the optimal convergence properties are observed for velocity, pressure and gradient of velocity both using triangular and quadrilateral meshes. Furthermore, it can be observed that, for the same level of mesh refinement, quadrilateral cells provide more accurate results when compared to meshes with triangular cells. Similar conclusions are obtained in three dimensions, where optimal rate of convergence is achieved in all cases for velocity, pressure and gradient of velocity and comparable accuracy is provided by all cell types.
The convergence study on stretched meshes with stretching factor s = 10 and s = 100, is performed for the Poisson problem. Figure 12 shows the relative error, measured in the L 2 (Ω) norm, of the primal and mixed variables as a function of the characteristic cell size. The results reveal that the accuracy of the proposed method is not dependent upon the stretching factor. The conclusions also hold for the Poisson problem in three dimensions, as illustrated by the results in figure 13. Optimal convergence is observed for all the variables and all cell types. The accuracy is again found to be almost insensitive to the stretching factor.

Computational cost
The last numerical experiment involves a study of the computational cost of the proposed method for meshes with different cell types. The computational efficiency is compared by directly measuring the CPU time (in seconds) required to assemble and solve the global system of equations, as this is the dominant cost of the proposed methodology. Figure 14 displays the evolution of the relative error of velocity, pressure and gradient of velocity in the L 2 (Ω) norm, as a function of the CPU time. The results reveal that quadrilateral and hybrid meshes provide the same accuracy as triangular meshes with slightly less computational effort. The better performance of quadrilateral cells is clearly observed when measuring the error of the velocity, whereas for pressure and gradient of velocity, all types of cells provide the same accuracy with a similar computational effort. It is worth noting that the advantages of using quadrilateral cells are not only observed when high accuracy is required. Even for an accuracy of 1% quadrilateral cells require nearly one order of magnitude less CPU time than triangles.
In three dimensions the conclusions are similar, as illustrated in Figure 15. Hexahedral and hybrid meshes are able to provide the solution with a given accuracy with slightly less computational effort when compared to tetrahedral, prismatic and pyramidal meshes. The most important differences are appreciated when the error of the velocity is considered.
It is worth noting that quadrilateral and hexahedral meshes seem to provide the maximum performance, in terms of achieving the desired accuracy with the minimum computational effort. However, it is well known that the mesh generation of complex objects using unstructured hexahedral meshes is still today an open problem. 38 In this scenario, the ability of the proposed method to handle hybrid meshes will be of use. As demonstrated by the results of Figures 14 and 15, the use of hybrid meshes is still beneficial when compared to pure triangular or tetrahedral meshes.

Two dimensional heat transfer problem with localised source
The model problem (1) in Ω = [0, 1] 2 is considered, where the source term and boundary data are selected such that the analytical solution is known and given by After computing the error indicator, as detailed in section 4, a desired size is determined for each cell of the coarse mesh. With this information, new meshes are generated, the first and second-order solutions are recomputed and the mesh adaptivity procedure is repeated. Figures 16 (b) and (e) display the meshes after two adaptivity iterations. The triangular mesh has 1,271 cells, whereas the quadrilateral mesh has 871 cells. As it can be observed, the approximations computed at the second iteration of the mesh adaptive process already capture the main feature of the solution, which is a Gaussian profile centred at (0.7,0.7). The adaptive process converges in seven iterations for triangular meshes and six iterations for quadrilateral meshes. show that the error indicator devised in section 4 produces a very accurate estimate of the error of the approximation u , computed by solving an inexpensive extra local problem given by equation (46), for both types of meshes. A slight difference is observed between the error indicator and the exact error of u for triangular meshes for the first four iterations of the mesh adaptivity process, whereas for quadrilateral cells a perfect agreement is observed.
To quantify this difference, the so-called indicator efficiency, computed as the ratio between the exact error of u and the error indicator (48), is reported in figure 18 (c). On the first iteration, the efficiency for triangular meshes is 0.76 and slowly improves during the adaptive process, being 0.82 in the third iteration, 0.98 in the fifth iteration and 1.01 in the final iteration. For quadrilateral cells, the indicator efficiency is already 1.01 in the first iteration and takes a value of 1.00 from the second to the sixth iteration. This clearly indicates that in this example, the use of quadrilateral meshes is beneficial compared to triangular meshes. Not only less iterations of the mesh adaptivity process are required but, in addition, the final error of the approximate solution is lower when using quadrilateral cells.
It is worth noting that despite the mesh adaptivity process is driven by computing an error indicator for the approximate solution u , the error of the more accurate approximation u also decreases monotonically during the adaptive process, as shown in figures 18 (a) and (b). This is expected due to the higher accuracy of the approximation computed with the proposed second-order FCFV, when compared to the accuracy of the first-order solution obtained by solving the extra local problem of equation (46). The difference in accuracy can be observed in figure 1, where the error of the first and second-order FCFV is compared for a Stokes problem in two and three dimensions.

Three dimensional Stokes flow around complex microswimmers
The last example considers the three dimensional Stokes flow around microswimmers. The geometry of the microswimmers, taken from, 40 is given in parametric form as where the curve C is a parametrisation of the centreline of the swimmer, namely C(λ) = (β cos(κλ), β sin(κλ), αλ) .
In (56), n 1 and n 2 denote the unit normal vectors to the centreline tangent and serve as the short and long axis respectively of the propeller cross-section. More precisely, they are defined as in terms of the Serret-Frenet normal N and bi-normal B. The radii of the long and short axis of the propeller cross-sections of the swimmer, denoted by R b and R n respectively, are defined as Figure 19: Geometry of three microswimmers for the parameters given in table 6, showing two perspectives of the same geometry.
where the function F s is defined, for s = {0, 1}, as All the parameters appearing in the previous expressions are given in table 6.  Table 6: Parameters used to define the geometry of the microswimmers.
Three microswimmers, obtained by varying the parameter γ in equation (58), are considered to demonstrate the potential of the automatic mesh adaptivity framework proposed in section 4. The geometries of the three cases considered are displayed in figure 19. 1 − x 2 1 )(L 2 2 − x 2 2 )/L 2 1 L 2 2 )µm/s is imposed on the inlet, at x 3 = L 3 , and a free-traction condition is enforced on the outlet, at x 3 = −L 3 . On the remaining lateral walls of B and on the surface of the microswimmer, a no-slip boundary condition, corresponding to material walls, is enforced. The kinematic viscosity is taken as ν = 2.65 mm 2 /s, which is the value corresponding to blood at 37 • C. The initial meshes displayed in figure 20 are generated, using the technique described in, 41 to start the automatic adaptivity process. To ensure a good geometric representation, the initial meshes are generated by imposing a desired element size of 0.04 µm along the curves S(λ, 0) and S(λ, π), for λ ∈ [−L, L]. The desired element size in the rest of the domain is 0.5 µm. The resulting meshes have 80,024, 80,294 and 77,533 elements for γ = 0, γ = π/4 and γ = π/2 respectively. It is worth noting that, due to the complexity of the geometry, only tetrahedral meshes are considered in this example.
The automatic mesh adaptive process is launched with a desired relative error of ε = 5 × 10 −2 . Convergence of the adaptivity procedure is achieved in four iterations for the three geometries considered. For the geometry corresponding to γ = π/4, figure 21 shows the second, third and fourth mesh obtained during the automatic adaptive process. The meshes in figure 21 have 193,159, 463,342, and 1,101,623 tetrahedrons respectively. As it can be observed, the refinement introduced by the automatic mesh adaptivity process concentrates the cells in the regions where the flow is more complex. The size of the corresponding global system to be solved to compute the velocity on the cell faces and the mean pressure on each cell is 1,338,781, 3,216,565 and 7,662,725, respectively.
The pressure distribution on the last mesh is displayed in figure 22. To offer a visual comparison between the three cases, the colour scale is adjusted in the three simulations to be in between −26.60 mPa and 33.24 mPa. The results clearly show the significant variation in the pressure distribution as the geometric configuration, controlled by the parameter γ, is changed. As the value of γ increases the reference area of the body of the swimmer increases, generating higher pressure values in the body for the geometry with γ = π/2 when compared to the geometry with γ = 0. In all cases, the maximum pressure is observed on the head of the swimmer and the magnitude of the maximum pressure is similar in all three configurations.
Finally, figure 23 shows the streamlines coloured with the magnitude of the velocity, with a minimum value of 0 µm/s and a maximum value of 4.92 µm/s. The complexity of the flow around the microswimmers  can be appreciated as well as the influence of the geometric parameter γ on the flow features.

Concluding remarks
This paper proposed a formulation of the second-order FCFV method for elliptic problems suitable for application in general meshes of triangular and quadrilateral cells in two dimensions and tetrahedral, hexahedral, prismatic and pyramidal cells in three dimensions. The computational cost of the resulting problem is comparable to the one of the original first-order FCFV method in terms of number of operations, since in both cases the global unknown is approximated with constant functions on the mesh faces. As in the original FCFV method, optimal first-order convergence of the gradient of the solution is achieved without the need to perform a reconstruction procedure. In addition, when CPU time is compared, the proposed method guarantees an improved approximation of the primal variable, which is now second-order accurate, of almost two orders of magnitude in three dimensional problems.
The proposed approach also inherits the robustness of the original first-order FCFV method in the incompressible limit. In addition, the proposed method is insensitive to the choice of the type of cells utilised in the discretisation, to their distortion and stretching. More precisely, successful simulations with hybrid meshes were also presented, paving the path towards the application of the discussed methodology to more complex problems requiring boundary layer meshes.
Finally, an automatic mesh adaptivity strategy is devised by means of a local error indicator obtained via the solution of one extra local problem per cell. The resulting cost of this strategy is thus limited, whereas its advantages are shown in presence of complex geometries and localised phenomenons as, starting from a coarse discretisation, the method is able to automatically construct a set of meshes to achieve a user defined target accuracy.
Extensive numerical simulations in two and three dimensions are presented to validate the proposed FCFV methodology and the mesh adaptivity procedure. Moreover, a three dimensional incompressible Stokes problem featuring geometries of interest in microfluidics applications is presented, showing the potential of the method, enhanced by the automatic mesh adaptivity strategy, to treat complex large scale flow problems.