Interior boundary-aligned unstructured grid generation and cell-centered versus vertex-centered CVD-MPFA performance

Grid generation for reservoir simulation must honor classical key constraints and be boundary aligned such that control-volume boundaries are aligned with geological features such as layers, shale barriers, fractures, faults, pinch-outs, and multilateral wells. An unstructured grid generation procedure is proposed that automates control-volume and/or control point boundary alignment and yields a PEBI-mesh both with respect to primal and dual (essentially PEBI) cells. In order to honor geological features in the primal configuration, we introduce the idea of protection circles, and to generate a dual-cell feature based grid, we construct halos around key geological features. The grids generated are employed to study comparative performance of cell-centred versus cell-vertex control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulations using equivalent degrees of freedom. The formulation of CVD-MPFA schemes in cell-centred and cell-vertex modes is analogous and requires switching control volume from primal to dual or vice versa together with appropriate data structures and boundary conditions. The relative benefits of both types of approximation, i.e., cell-centred versus vertex-centred, are made clear in terms of flow resolution and degrees of freedom required.


Introduction
Subsurface reservoirs are often comprised of complex geometric and geologic objects and features. In addition to robust numerical methods for solving the flow equations, grid generation methods are required which can honor geometric complexity and permit local grid cell density control. Grid generation for large-scale porous media poses the challenge of complex geometries and random distribution of spatial heterogeneities in the domain, e.g., [42,69].
Standard reservoir simulators were originally based on simple grid blocks, i.e., squares and cubes primarily using structured grids [6]. Although it is relatively easy to implement simulation techniques on such simple grids, they inherently lack the ability to adapt to general geological features and complex geometries [31]. Unstructured grid generation offers the desired flexibility by employing triangles and tetrahedra as grid elements. Unstructured grids allow grid cells to adapt to various flow and geometric constraints and also local refinement with smooth transition [31,58]. However, unstructured grids require special data structures and are computationally more involved.
Despite unstructured grid generation methods having been successfully employed in modelling complex giant reservoirs, in field applications, there is still increased inclination toward the use of structured grids. Fung et al. [31] have reported that this might be the result of novelty of these methods in the field compared with structured grids for which well-established simulation tools exist, and consequently, more research work is required for simulation on unstructured grids.
The grids generated must be compatible with the flux approximation schemes employed. The two-point flux approximation (TPFA) involving a two-point pressure difference to approximate the flux across the face of each control-volume, is still regarded as the standard reservoir simulation flux approximation [7]. However, TPFA has consistency limitations which are discussed below. The schemes employed here have been designed to overcome the TPFA limitations and involve families of control-volume distributed multi-point flux approximations (CVD-MPFA), with flow variables and rock properties being controlvolume distributed and assigned to share the same controlvolume locations. Both cell-vertex and cell-centred CVD-MPFA formulations [19][20][21][22][23][24][25][26][27][28] are compared in this work on unstructured grids. Related cell-centred MPFA methods are presented in [1-3, 44, 48, 80]. We note that alternative methods have been proposed to improve over standard TPFA including [63] which involves a hybrid approximation and [50] which imposes a maximum principle via a non-linear formulation. The unstructured grids used in reservoir modelling, commonly employ Delaunay-voronoi (PEBI) grids for spatial discretization of domain. The grids generated are comprised of simplices. The dual voronoi mesh is obtained by joining circumcentres of primal grid simplices attached to a common grid vertex. Consequently, a natural choice is for primal grid cells to act as control-volumes, then grid generation can be performed with primal grid cell boundaries being aligned with key interior constraint boundaries. This naturally leads to cell-centred approximation, where flow variables and rock properties are associated with grid cell centres.
The alternative is to employ cell-vertex approximation which uses far fewer approximation points on a given unstructured grid. In this case control-volumes are constructed around primal grid vertices. The grid generation process is less straight forward since control-volumes must be constrained to satisfy interior boundary alignment. A novel grid generation procedure is proposed that automates control-volume boundary alignment and yields an essentially voronoi mesh. The actual grid is then generated such that dual-cell boundaries are aligned with key internal constraint boundaries and the cell-vertex approximation becomes the natural choice. In this case, flow variables and rock properties are associated with grid cell vertices and their associated control-volumes resulting from the dual mesh.
In this work, development of boundary aligned unstructured grid generation technique is presented, which can be employed to generate quality meshes for reservoir geometries. The unstructured grid generation techniques presented offer the advantage of being equally capable of generating boundary aligned grids(BAGs) and well aligned grids(WAGs) both with respect to primal and dual(voronoi) cells. Grids thus generated are employed to simulate pressure fields, and a comparative study of cell-centred versus cell-vertex control volume distributed(CVD) flux approximation schemes is presented. We begin with a brief review of methods for generating feature based unstructured grids and Delaunay triangulations. In the next section, description of the grid generation techniques together with generation of BAGs and WAGs is presented. This is followed by a brief description of the control-volume distributed multipoint flux approximation (CVD-MPFA) schemes, including the earlier piecewise-linear based formulations [19,20,23,24], and more robust piecewise-bilinear based formulations [22,25,27,28]. A comparison between the cell-centred and cell-vertex CVD-MPFA results and corresponding grids is presented, followed by conclusions.

Methods for geological feature based grids
In order to minimize the effect of grid orientation and discretization errors in simulation of hydrocarbon flows, unstructured grids should conform as closely as possible to geological features such as faults, fractures, layers, pinch-outs and wells. A suitable grid generation method for reservoir simulation (e.g., [58]) should honor the following geological features: -Upscaling consistency: Geological layers are relatively homogeneous in nature with respect to permeability and porosity. However, across the layers, these physical properties may jump by several orders of magnitudes. Therefore, control-volumes should completely reside in layers, and not contain interfaces between layers or features. -Faults: Faults are either barriers or drains and often modelled with transmissibility multipliers [58]. Thus, in reservoir modelling, grids generated should be fault aligned. -Fractures: Fractured reservoirs are characterized by the presence of two distinct types of porous media: matrix and fracture [5]. The matrix and fractured media have very different fluid storage and conductivity characteristics; consequently, fracture aligned grids are required, to minimize discretization error. -Pinch-outs recovery: In reservoir simulation, pinch-outs are frequent, so any grid generation method must ensure that control-volumes conform to the respective interior geometric boundaries. -Well trajectories: Wells act as either point or line sources and often drive the flow field to be simulated. In order to minimize discretization error, grids are gene-rated such that control points of control volumes honour the resulting sequentially defined well path trajectories.
Unstructured grids used in reservoir modelling, commonly employ Delaunay-voronoi (PEBI) grids for spatial discretization of domain [60]. Voronoi grids can be made to conform to geological features by special treatment such that their cells become aligned to these features. In contrast, flow-based grids concentrate local grid refinement in high flow regions and apply to heterogeneous media, e.g., [11,20,79], vorticity based grids have a similar objective [54]. Here, we focus on constructing feature based voronoi meshes, some of the most commonly employed techniques are outlined below.
Equidistant seeds In this method, the features to be recovered, can be thought of as those aligned with the edges of voronoi cells. Thus, to recover a feature, two equidistant generating points(seeds) are introduced on either side of the feature, in a direction normal to the feature surface. Construction of voronoi cells corresponding to these seeds, ensure that the desired local feature is retained as a boundary edge of the voronoi cells [76]. This method works well for simple features; however, for complex geometries and pinch-outs, limitations are noted [58].

Geological features as Planar Straight Line Graph (PSLG)
The idea behind this method is to construct a protection area, usually created by advancing from PSLG thereby building up a layer of voronoi cells on each side of PSLG [13]. Once such a protection layer is built, the remaining domain is meshed by constrained (non-conformal) Delaunay triangulation, and consequently the perpendicular bisectional(PEBI) property is not preserved. In [9] a method has been proposed which also relies on the construction of protection region to honour geological features.
Optimized coordinates of voronoi seeds This method is described in [58], according to which, grids involving voronoi cells are first generated, which are then conformed to geological features by recursive optimization of voronoi seeds based on the concept of centroidal voronoi tessellation(CVT) [17]. Since this method relies on an iterative technique, it can only recover features with a close proximity.

Region-based optimum triangulation
In [31], a method has been proposed which allows cells to be optimally placed, within the solution domain in the region of interest only. This method employs a structured background mesh, which is then converted into an unstructured mesh in regions of interest such as wells, where polygons in two and polyhedra in three dimensions are extracted. These near well regions are resolved by employing local grid refinement. This method allows the grid to retain the properties of a structured grid, with unstructured grid dominant in the region of interest only. A growing region technique which also involves a region-based optimization has been presented in [42].
Proposed grid generation methods In this work, we present methods to generate both primal and dual-cell feature preserved Delaunay admissible triangulations [33,55,66,72,73,78]. In order to construct primal-cell BAG, we introduce the idea of protection circles build around key geological features, this ensures integrity of features in the final Delaunay triangulation. The idea of protection-circles is unique and does not impair the PEBI property associated with Delaunay meshes. For dual-cell feature-based grids, it is a prerequisite to build a (protection) layer of quadrilaterals, with a geological feature defining the median line inside the quadrilateral layer, we call this additional protection layer a halo. A halo protected primal mesh is essentially Delaunay and honours geological features with respect to dual-cells. For comparison, a list of some opensource meshing tools are given in Appendix B together with a brief discussion.

Delaunay triangulation and voronoi (PEBI) meshes
Grids based on voronoi diagrams remain predominant in reservoir simulation [9,31,58,60,61,64,76]. Voronoibased grids are locally orthogonal and permit two-point flux approximation if the permeability field is isotropic, or if the grid is generated to be K-orthogonal [64]. Voronoi grids are also called Dirichlet or Thiessen tessellations. Mathematically, each voronoi region v i , associated with a site s i is defined by: where d(p, s i ) is the Euclidean distance between a point p belonging to v i and associated site s i . The boundary between two sites s 1 and s 2 , also called a voronoi edge, and is a line segment given by points(p) which satisfy d(p, s 1 ) = d(p, s 2 ). Delaunay triangulations lead to the construction of PEBI(PErpendicular BIsectional) grids, e.g., see Fig. 1c.

Duality of Delaunay triangulations and voronoi grids
Delaunay triangulation(DT) is the dual of the voronoi diagram, given a DT or a voronoi grid, the dual can be constructed e.g. [66,72,73,78]. From the grid generation view point, it is relatively easy to work with Delaunay triangulation directly, and voronoi grids if required are obtained as the circumcentre dual of the associated Delaunay  triangulation. The Delaunay triangulation is equivalently defined by Boris Delaunay through the empty circle property [46]. A triangulation is Delaunay if no simplex circumcircle contains any site in its interior e.g. [33,34,38,78]. Figure 1 shows a simple example of Delaunay triangulation and associated voronoi grid, constructed by joining circumcentres of cells having an edge adjacency.

Construction of Delaunay triangulation
There are several algorithms for construction of Delaunay triangulation [29,35,47], among others the incremental insertion [8,36,74,77] is the most widely used technique. It is a simple and flexible technique in that its extension to higher space dimensions is relatively straight forward [75]. The incremental insertion algorithm, involves triangulating each point one at a time, in an existing triangulation. In order to start the triangulation three/four phantom points are created thereby leading to the formation of a convex bounding box, which if not simplex, is first divided into simplices. Then, in an existing background mesh, for every newly introduced point incremental insertion requires: locating the containing element/elements(base); wherein the point at hand is triangulated by locally breaking the Delaunay chain; followed by local reconstruction carried out to restore the Delaunay property. This process is repeated until all points are exhausted. The two variants of incremental insertion algorithm, namely Watson's [77] and Green-Sibson [36], are the most commonly employed algorithms for construction of Delaunay triangulation. The Green-Sibson algorithm is more general in a sense that it can be used with any user defined connection optimization criterion to construct a data dependent triangulation [8]. In Appendix C, the procedure Incremental-Insertion outlines key steps involved in the Green-Sibson algorithm employed to construct Delaunay triangulation, these steps are illustrated in Fig. 2. In the Green-Sibson algorithm, the connection optimization phase is carried out through local reconstruction thereby swapping edges(2D)/faces(3D) subject to violation of a quality criterion [41,45]. If a Delaunay triangulation is desired then connections are improved subject to Delaunay measure(α).
(a) Highlighted is the base of candidate   Delaunay measure A triangulation is Delaunay if no simplex circumcircle contains any site in its interior. In order to ensure that the circumcircle of every simplex is empty, Delaunay triangulation relies on the Delaunay measure(α) [32]. Figure 3a represents a triangle v 1 v 2 v 3 with r being the radius of the circle circumscribing v 1 v 2 v 3 , and there exists a query point q. The Delaunay measure(α) is defined as: where d(q, O) is the Euclidean distance between query point(q) and circumcentre(O) of v 1 v 2 v 3 , and r v 1 v 2 v 3 represents circumradius of v 1 v 2 v 3 . Note that α < 1 implies that query point(q) falls inside the circumcircle. A triangulation comprised of n simplices and m vertices is represented as: where t corresponds to triangle(2D)/tetrahedron(3D) and v represents vertices. In D dimensions each t i is defined by a unique and distinct set of D + 1 vertices called connectivity of simplex t i . A triangulation T is regarded as a Delaunay triangulation provided where α t i is the Delaunay measure for simplex t i comprised of (v a , v b , v c ) vertices in 2D. For a Delaunay triangulation α t i > 1 holds true computed with any grid point(v p ) other than those constituting simplex(t i ) under consideration, e.g. Fig 3b.

Boundary integrity with conformal and non-conformal Delaunay triangulation
Delaunay triangulation is the triangulation of a convex hull of a predefined data set. In a Delaunay triangulation, connections are improved such that none of the simplices contains any site in its interior, but it cannot be guaranteed that connections between the given point set are present in a prescribed manner [78]. Thus, to ensure integrity of interior boundaries and those of outer boundaries, it is mandatory to couple the Delaunay triangulation algorithm with a boundary/feature recovery technique. This limitation of DT is well known and several methods have been proposed which can be employed to retrieve missing boundary connections. In what follows we briefly outline two commonly employed methods including local refinement [78] also called the stitching method and edge swapping [71] which are used to constrain the mesh to honour a feature. These methods are equally applicable to recover boundary and/or desired field(interior) edges. Consider the simple Delaunay triangulation shown in Fig. 4, where a field edge ab has been honoured The local refinement technique requires introduction of new points midway between each missing connection, recursively, until the desired edge is recovered. On the other hand, edge swapping involves swapping edges intersected by segment ab recursively until the feature is recovered.  are pictorial representations of field edge(ab) constrained triangulations, achieved through local refinement and edge swapping respectively. It can be noticed that the former of these boundary recovery approaches results in the conformal Delaunay triangulation, i.e., PEBI-Grid, while the later yields a non-conformal Delaunay triangulation. In the meshing technique (described subsequently), we employ the former, i.e., use local refinement to recover missing boundaries, and/or features in the empty mesh, which are then preserved throughout the field mesh generation process.

Delaunay triangulation and field point placement
In order to start with any triangulation, a prerequisite is to construct an empty mesh, comprised of a prescribed set of points defining domain boundaries and/or geological features. We use the Green-Sibson algorithm in conjunction with the incircle criterion for connection optimization, which yields the Delaunay empty mesh of predefined set of points, e.g., see Fig. 5a-b. Delaunay is a criterion to connect a predefined set of points, which of course is optimal in many aspects [45,70]. The empty mesh constructed from prescribed boundary discretization comprises of lowquality elements, e.g., Fig. 5b. In order to generate a quality well-resolved grid, the empty mesh is locally refined by introducing new field points. Several options exist relating to introducing field (interior) points in an empty mesh, so as to design an automatic Delaunay triangulation algorithm. In general point placement strategies involve introducing field points at either the centroid [38], circumcentre [67], or by edge subdivision [33] of the elements in the background mesh. This process of point placement is used repetitively, until the desired edge length or some other criterion measuring element size is satisfied. In Delaunay triangulation whenever point placement strategies are discussed, one must include another class of unstructured mesh generation techniques called the advancing front method [52], which provides an optimal point placement strategy. The advancing front method starts with a valid boundary discretization comprised of edges(2D) regarded as initial fronts. The advancing front method operates on each front one at a time locate an existing point or introduces a new field point in a direction normal to the front so as to construct a quality simplicial mesh [52,59]. The advancing front method constructs a mesh by generating each element one at a time, where in order to validate each newly created element a check for intersection with existing elements is required. Due to these reasons the advancing front method is not only inefficient but also suffers from robustness issues [56,57]. It is well established that DT has a sound mathematical basis, while its counter part, the advancing front method provides optimal point placement. The idea to combine these two methods into one technique was introduced in the nineties [39,56,57,59]. In such a hybrid method, field points are introduced in a manner similar to the advancing front method, while their connections are improved by enforcing the Delaunay criterion. The advancing front method used in conjunction with Delaunay criterion both provides optimal point placement, and makes it relatively easy to obtain boundary aligned grids, where integrity of boundaries is a further advantage.
Uniform isotropic mesh of a unit square domain with different field point placements A test case involving an empty mesh formed over a unit square domain is selected to illustrate different field point creation strategies, e.g., see Fig. 5b. In order to start with field mesh generation in the empty mesh a metric [32] is assigned to each point, and boundary point spacing is used for each point to construct its associated metric. For the case under consideration, metric specification turns out to be isotropic and uniform, this is because all boundary points are equally spaced(h = 0.10 units). The domain is meshed either by introducing candidate(field) points at the centroid, or circumcentre, or through edge subdivisions of the elements of the background(empty) mesh, or via the advancing front point placement technique. The grids generated are shown in Fig. 5c-f, respectively. For comparison purposes in all cases, field points are introduced iteratively until a unit metric length mesh [34] is generated. The candidate points introduced during each level are filtered to remove conflict with existing and/or previously accepted points. We use the same filter in all of the cases, and accepted points are connected by using Green-Sibson algorithm subject to the incircle quality criterion test. Note that no mesh cosmetics are applied to the final triangulations. From Fig. 5f, it is observed that, with the exception of regions where fronts collide with each other, the advancing front point creation technique provides optimal point placement and thus yields a quality mesh compared to other counterpart field point creation methods. Furthermore, in the advancing front point placement method, fronts can be initiated from interior boundaries and is therefore well suited for meshing reservoir geometries. In the following, we use advancing front point placement in conjunction with incircle criterion to generate Delaunay triangulation.

Geological feature-based grid generation
In general, reservoir geometries are comprised of various features such as faults, fractures, pinch-outs and layered media, with a wide range of variations in porosity and permeability across different layers e.g. [69], in addition reservoirs can have a complex spatial distribution of wells in place. In order to minimize discretization error, it is common practice to generate meshes which are aligned with such features, thereby leading to feature based triangulations. The grids generated, provide a means of spatial discretization required by flux approximation schemes. The flux approximation schemes used in reservoir simulation are control volume distributed, where a piecewise constant representation of flow variables and rock properties is assigned to control volumes, so that field variables are located at their centres, and/or circumcentres also called control points [64].
Classification of geological feature-based grids For the purpose of grid generation, the geological features can be classified into two groups; the first of which deals with domains involving layers, fractures and/or faults, and the second treats well distribution (Fig. 6). In the former, control-volume aligned grids also known as boundary aligned grids(BAGs) are used, while in the latter, controlpoint aligned grids will be used, which we will call well aligned grids(WAGs). Furthermore, the geological features can be honored either with respect to elements constituting a simplicial(primal) mesh, or with respect to dual(voronoi) cells constructed from an underlying simplicial mesh. The grids generated are termed primal-cell and dual-cell feature aligned grids respectively. Figure 6 displays the classification of feature based grids. Referring to Fig. 6 it is observed that the dual of a primal-cell BAG corresponds to dualcell WAG, and primal-cell WAG becomes dual-cell BAG viewed in dual configuration. In the following we describe methods for generating boundary and well aligned grids  Fig. 6 Classification of geological feature-based grids. In BAG control volume faces are aligned with the geological features, whereas in WAG control volumes are generated such that their control points when joined sequentially respect well trajectories with respect to primal-cells, a dual constructed from an underlying primal-cell feature based grid is used to honour features in the corresponding dual configuration. To honor geological features in primal and dual frameworks, we introduce the ideas of protection circles and quadrilateral halo construction, respectively.
Empty mesh generation and ensuring integrity of boundaries and/or embedded features Regardless of the type of feature aligned Delaunay grid, i.e., BAG or WAG, a prerequisite is to generate a feature preserved Delaunay empty mesh. Curves characterizing domain boundaries and geological features are embedded in discrete form, e.g., see Fig. 7a. The Delaunay triangulation of a point set defining domain boundaries and features leads to an empty mesh, where integrity of features may not be preserved. Figure 7b displays the Delaunay empty mesh of a test case, where a desired connection constituting features to be honoured between highlighted points is missing. In an empty mesh, the integrity of embedded curves is ensured by employing the local refinement method(described above), this yields a feature preserved Delaunay empty mesh, e.g., see Fig. 7c.

Primal-cell BAG (dual-Cell WAG) generation
Boundary-aligned grids (BAGs) are those grids in which control volumes are generated such that their boundaries honour geological features, control-volume faces crossing the features are not allowed. In order to honor features with respect to primal-cell control volumes faces, i.e., triangles and/or quadrilaterals cells, a prerequisite is to generate a feature preserved Delaunay empty mesh (Fig. 7). The feature-aligned Delaunay empty mesh comprises of low quality elements. In order to generate a well-resolved mesh comprised of quality simplices, the empty mesh is refined by introducing new(field) points. The Delaunay triangulation of new field points encroaching simplices honouring geological features can lead to reconfiguration of these connections, destroying boundary integrity. In order to avoid swapping and to preserve integrity of features honoured in the empty mesh, we introduce the idea of protection circles that pass through the simplices constituting geological features. In the field mesh generation integrity of features is maintained provided new points encroaching the protection circles are not accepted.

Protecting interior boundaries in the empty mesh by protection circles
In order to generate primal-cell boundary aligned Delaunay grids, starting from a feature honoured empty mesh, field mesh generation is carried out such that edges defining geological features are not reconfigured. To this end, we protect edges constituting features by passing empty-circles through the respective edge vertices. The empty-circles are called protection circles, and they are either diametric circles of the feature edges, or as explained below, they are simplex circles, where a feature edge forms one side of the simplex.

Delaunay admissible simplexes and idea of protectioncircle
In two dimensions, a simplicial mesh is comprised of points (0-D simplexes), edges (1-D simplexes), and triangles (2-D simplexes). In Delaunay triangulation, connections between simplexes are established such that circles circumscribing triangles are empty, i.e., circumcircle of a triangle does not contain any site in its interior. A simplex (edge/triangle) whose smallest circle is empty is Delaunay admissible and exists in Delaunay triangulation [33]. In two dimensions (D = 2), the D−1 entities(simplexes) are edges. Consider an arbitrary 2D triangulation shown in Fig. 8a. In Fig. 8a, the edge ab is not Delaunay admissible because the diametric circle of the edge contains points of simplexes sharing edge ab; edge cd is local Delaunay, but not global, since the diametric circle contains point p; and edge de preserves the Delaunay property both in a local and global (a) Domain boundaries and feature discretization.
(b) Empty Delaunay mesh with feature to be honoured has been triangulated, also shown is the close up view. Note that a desired connection between highlighted points is missing.
(c) Local refinement is employed to retrieve missing connection. It requires one new point introduced as mid point of the missing connection to recover missing edge, as displayed in the close up view  perspective. In an arbitrary triangulation of a data set, if connections between D − 1 simplexes are (global) Delaunay admissible then these connections exist in the Delaunay triangulation [33] of the data set. The Delaunay triangulation of the point set defining arbitrary triangulation displayed in Fig. 8a, generated by using the Green-Sibson's algorithm subject to the incircle test as a quality criterion is shown in Fig. 8b. We note that among three representative selected edges, i.e. ab, cd, and de ( Fig. 8a), the edge de is Delaunay admissible and thus exists in the Delaunay triangulation, e.g., see Fig. 8b. In two dimensions in an arbitrary triangulation of a given point set, edges whose diametric circles are initially empty, always exist in the Delaunay triangulation of the same data set [14]. Nevertheless, in the Delaunay triangulation of a data set, the diametric circle of an edge can have one or more points contained in it, e.g. see Fig. 8c the diametric circle of edge pq contains point c in its interior. For a triangulation where the diametric circle of an edge has points in its interior, then for it to be Delaunay, there must exist a triangle formed by joining the edge at hand with a point contained in the diametric circle such that its circumcircle is empty. In Fig. 8c, the circumcircle of triangle pqc is empty. The Delaunay admissibility of the simplexes constituting a Delaunay mesh is summarised in following theorem:

Theorem 1 In two dimensions, for a D − 1 simplex(edge) to be part of Delaunay triangulation either its smallest circle (diametric-circle) is empty or there exists a 2-D-simplex with empty smallest circle (circumcircle), defined by connecting the edge vertices to the nearest point inside the diametric circle.
Notion of protection circle To start with meshing, a prerequisite is to generate a feature preserved empty Delaunay mesh. Delaunay triangulation of a point set defining domain boundaries and features (i.e., data set defining input surface and/or curve mesh) leads to an empty mesh, where integrity of features may not be preserved. To generate PEBI grids, in an empty mesh integrity of embedded surfaces and/or curves is ensured by employing the local refinement (conformal boundary recovery) method(described above), this yields a feature preserved Delaunay empty mesh. The feature recovered Delaunay empty mesh is comprised of low quality elements. In order to generate a well-resolved mesh comprised of quality simplexes, the empty mesh is refined by introducing new(field) points. The Delaunay triangulation of new field points encroaching simplexes honoring geological features can lead to reconfiguration of these connections, destroying boundary integrity. In order to avoid swapping and to preserve integrity of features honoured in the empty mesh, we introduce the idea of a protection circle that passes through the simplexes constituting geological features. In field mesh generation, integrity of features is maintained provided new points encroaching the protection circle are not accepted. The protection circle used is either diametric or a simplex-circle.
Protection circle is either a diametric circle or circumcircle In carrying out field mesh generation, the diametric circle can only ensure integrity of a feature (edge) provided it is empty in the corresponding empty-mesh, e.g. see Fig. 9. For edges to be protected the protection circles are defined by the respective (edge) diametric circles if they are empty, otherwise they are defined by the circumcircles of the triangles resulting from connecting the (feature) edge vertices with the corresponding point contained in the respective diametric circles. If edges whose diametric circle contains point(s) in their interior constitute Delaunay (b) Triangulation of a field point p, direct insertion followed by optimization subject to incircle test. Note that the edge ab has been reconfigured.
triangulation then existence of triangles formed by joining these edges with a point contained in their diametric circle is guaranteed [14,18,33]. This is because the Delaunay triangulation is the nearest neighbour graph [14]. Note that the circumcircle is always empty in a Delaunay triangulation.

Primal-cell feature honoured mesh generation
Starting with an empty mesh, a field triangulation is carried out, such that the Delaunay property of edges defining geological features is preserved. To this end edges constituting features are marked as Delaunay admissible edges and are protected by constructing protection circles around them. Figure 10 displays a test case feature honoured empty mesh, with diametric circles drawn around the edges constituting the feature. We note (Fig. 10a) that there exists two diametric circles (highlighted) containing a point p in their interior. In order to maintain integrity of features the diametric circle of edges containing a point in their interior are replaced with associated circumcircles, e.g., see Fig. 10b. The empty mesh delineated in Fig. 11a with disks drawn around simplices constituting features to be honoured, defines a region where no field point should be accepted. The field mesh is generated by employing advancing front point placement used in conjunction with incircle criterion. The empty mesh is refined by introducing new points in an iterative manner. The mesh after first level of point placement is shown in Fig. 11b. The final primal-cell feature honoured mesh is displayed in Fig. 11c. The above grid generation processes provide a framework for primal-cell boundary aligned grid (BAG) generation. As described above, boundary-aligned grids(BAGs) are those grids in which control-volumes are generated such that their boundaries honour geological features, controlvolume faces crossing the features are not allowed. The proposed triangulation method can be used to obtain meshes which are aligned with respect to geological layers, faults and fractures. In this work, we first generate an empty mesh and then honour any missing connections corresponding to boundaries and/or geological features. The features (a) Delaunay empty mesh with diametric circle as protection circle around the features. Note that diametric circles(highlighted) contain a point p in their interior.
(b) Delaunay empty mesh with protection circles. Note that diametric circle of edges containing points in their interior are replaced with associated circumcircles.  honoured in the empty mesh, are protected during field mesh generation. To this end, the honoured edges are marked as Delaunay admissible simplices and are protected by enclosing them with protection circles. During field mesh generation, new field points encroaching any protection circles are not accepted.

Construct dual-cell WAG from underlying primal-cell BAG
In constructing a dual-cell well aligned grid(WAG), the first step requires embedding curves representing well trajectories in the reservoir domain, in discretized form. Then, in a manner similar to primal cell boundary-aligned grid generation, meshing of the reservoir domain is carried out. In this case, the target mesh is dual, as opposed to primal. Figure 12 shows a primal-cell BAG honouring highlighted feature together with its dual mesh. It can be seen ( Fig. 12) that control points of the highlighted dual cells define the embedded feature as a well trajectory. Figure 13 highlights key steps involved in generating a dual-cell well aligned grid from an underlying primal-cell BAG.

Primal-Cell WAG (Dual-Cell BAG) generation:
A procedure for generating primal-cell well aligned meshes can be constructed by employing the proposed meshing technique is outlined. While generating dual-cell WAGs, it has been demonstrated that well-aligned grids(WAGs) require control point alignment to a well pattern (Fig. 13).
Regardless of the selection of control-volume, i.e., primal/dual, the control point is always defined as the centre of the control-volume. For a well-aligned mesh, we are interested in positioning control points such that when joined sequentially, they follow a predefined well trajectory. The curves characterizing well trajectories require special treatment. This is because well alignment is ensured via control points, and not with control volume faces. This naturally leads to protecting well trajectories by building a protection layer such that control points of the elements constituting the protection layer define well paths. To this end, we propose the idea of a quadrilateral protection layer, constructed such that it encloses a well trajectory as its medial line, the additional protection is called a halo construction.
Halo construction In order to honor well paths (or generate dual-cell BAG), curves characterizing geological features require special treatment. Beginning with an empty mesh where features are embedded in discretized form, before field mesh generation takes place, halos comprised of quadrilateral cells are constructed such that honoured features can be retrieved as medial lines of halo elements.
To elucidate construction of a halo, consider a feature discretely defined by points p 1 p 2 p 3 p 4 as shown in Fig. 14a. For simplicity halo construction is first demonstrated independently of an empty mesh. The halo construction is performed by pushing the underlying feature point in directions normal (upward and downward) to the feature edges in an average sense, a distance assumed to be a fraction (γ ≈ 0.10) of average length of the feature edges attached (c) Dual-cell(voronoi) well aligned grid .
to the point to be split. The unit normal used to split a feature point say p i is computed from associated edge normals by: where n e represents edge normal, e.g. see Fig. 14a. While splitting we update the position of underling point (from p i to p i ) and simultaneously generate a new point(q i ), given by: where the factor β is defined as: and l e represents edge length (see Fig. 14a). This simple halo construction procedure provides a framework for generating primal-cell WAG (or dual-cell BAG).
To start with halo construction, a spear comprised of two triangles is first built as shown in Fig. 14b, which is then propagated until the last point (p 4 ) of the feature is encountered, e.g., see Fig. 14c. The halo is propagated via the use of spear triangles which enables quadrilateral cell construction (one by one). Construction of the halo is demonstrated in Fig. 14b-c, where points defining features are split by the spear and enclosed by the resulting quadrilaterals. Finally, Since the feature honored can be retrieved as the medial line of the halo, and is thus honoured in the dual configuration. The dual of a primal-cell WAG becomes a dual-cell BAG viewed in dual framework, e.g., see Fig. 16. The halo construction is reminiscent of the advancing front method, and its construction in two dimensions (i.e., D = 2) can be explained as being constructed by sweeping the 1-D simplex (edge) along the geological feature.
Halo construction in the empty-mesh As described earlier, to start with triangulation of a reservoir domain, we construct an empty mesh comprising of a prescribed data set. From the mesh construction view point, it is relatively simple to construct halos around the desired features in the empty mesh, i.e., halo construction is dealt with before invoking field mesh generation. Figure 15 displays key steps involved in constructing halo in an empty mesh. We start with the spear, which is propagated constructing a quad cell (one by one), until the feature is fully enclosed by opening spear ends, these steps are delineated in Fig. 15a-d, respectively. After features to be protected are enclosed with halos, field mesh generation is performed without impairing connectivity of halo quads. To this end, the edges constituting halo quads are constrained, and during field mesh generation any new(field) point falling inside halo quads is not accepted. Furthermore, for field meshing, we use advancing front point placement which initiates fronts from these quad edges, providing protection around halo quad edges. The resulting halo protected feature based grid is shown in Fig. 15e. The dual of the halo protected grid corresponds to dual-cell BAG, e.g., see Fig. 15f.
The halo used is comprised of quadrilateral elements, in order to yield a PEBI mesh each quad element should be allowed to have only four co-circular points [34] (Fig. 14d). For simple features, e.g., Fig. 16 it is possible to generate a PEBI mesh comprised of triangles and quads. However, in complex configurations and in higher space dimensions, it is almost impossible to yield a mixed element PEBI mesh [9,12]. In order to design a robust halo construction method, we relax the PEBI property locally for halo quads only. Therefore, halo protected meshes in general lead to essentially PEBI grids.  Halo construction around intersecting features The halo construction around intersecting features can be easily dealt with by using spear triangles. The intersection point where all of the features meet is also called a junction point. Figure 17a depicts a test case mesh (empty) where features intersect at junction point. We propose to start the halo construction from the junction point, for the case at hand the halo method leads to the mesh shown in Fig. 17b. Once the halo construction has taken place, the junction point is merged thereby introducing a polygon with number of edges equal to number of curves emanating from the junction   Fig. 17b, there are four branches intersecting at the junction point, thus a quadrilateral is used to merge the junction point, e.g., see Fig. 17c. Finally, spear ends are opened to fully enclose the honoured feature (Fig. 17d). The halo protected field mesh generated by employing advancing front field point placement is shown in Fig. 17e. The dual of a halo protected grid corresponds to a dual-cell BAG, e.g., see Fig. 17f.
Halo construction as an integral part of boundary/curve meshing A mesh generation process involves three key steps, i.e., (i) curve or boundary meshing, (ii)empty mesh generation coupled with boundary recovery, and (iii) field mesh generation. Halo construction can be performed at any stage of the mesh generation process. Halo construction if performed after field meshing is computationally more involved, since it requires to interact with full field mesh. In the grid generation methods described above, halo construction is performed in the empty mesh and does not require boundary recovery as an additional step. Figure 18 displays a case where halo construction has been performed as an integral part of boundary or curve meshing process.
In the empty mesh boundary recovery of the halo-quads is performed, followed by field meshing, e.g., see Fig. 18. Figure 19 summarizes the proposed cell-centred and/or vertex-centred grid generation methods. Geological featurebased grid generation involves three key steps, described below.

Key components of the proposed gridding methods
-Prerequisites: Data set defining domain boundaries and/or geological features is assumed as input. In D dimensions D − 1 (boundary) mesh generation is performed defining domain boundaries and features of interest, achieving desired uniform and/or non-uniform refinement. To start with a triangulation process, first the point set defining D −1 boundaries and features are triangulated, this is followed by boundary and/or feature recovery. -Novel components: We propose the idea of protection circles and halo construction performed around key geological features. Grid generation is performed by employing primal-cells (triangle and/or quad) as grid elements. Halo construction is required to honour geological features with respect to primal-cell control points (centroids), whereas protection circles are used to protect primal-cell faces. Halo construction and protection circles are required both in cell-centred and vertex-centred grid generation methods, e.g., see Fig. 19. -Field meshing: After geological features are enclosed with halos and/or protection circles, field meshing is performed generating a quality mesh with smooth

Examples of BAGs and WAGs generated using the new meshing techniques
The algorithm designed to mesh reservoir geometries is summarized in the procedure Feature Based PEBI-Grid Generator, listed in Appendix C. In order to illustrate the capabilities of the proposed meshing technique, the following examples are considered.

Both boundary and well-aligned grids
As illustrated above, from the construction view point, boundary-aligned and well-aligned grids are quite similar, they only differ in the setting that they are viewed, i.e., primal/dual. Nevertheless, when they appear in a group, their differences and requirements are evident. This case is selected so as to model a situation which involves both interior boundaries and predefined well paths. Figure 20, is the  Fig. 21a, b, respectively. In the grids generated, we note that control points follow prescribed well trajectories  embedded at the centre of domain, whereas along the layers control-volume alignment has been ensured.

Layered system honored boundary-aligned grids
Many oil and gas reservoirs are stratified, caused by sedimentation spanned over a long period of time. In order to check the capabilities of the triangulation technique, a synthetic layered reservoir is designed, motivated by [40]. A cross sectional view of the model, mainly driven by seismic attributes, contains a layered system with embedded features. We generate grids respecting the highlighted geological features both with respect to primal and dual cells.
Curves characterizing highlighted features with close proximity, are embedded in discretized form, and the resulting primal and dual-cell BAGs generated by employing the proposed gridding methods are shown in Fig. 22a, b, respectively.

Complex multilateral-well aligned grids
A representative complex multilateral well-path has been selected similar to [31], and the resulting feature honoured primal and dual-cell grids generated with the above techniques are shown in Fig. 23. We note that grids generated conform to the prescribed well trajectory both with respect to primal and dual cells. Note that for the boundary aligned grid, where the complex well pattern can be thought similar to a complex fracture, the same procedure holds, however, settings will now be switched, i.e., a primal cell wellaligned grid, becomes a dual-cell boundary aligned grid when viewed in a dual configuration.

Fracture network honoured grids
Naturally fractured reservoirs are characterized by the presence of two distinct types of porous media: matrix (a) Layered system honoured primal-cell BAG.
(b) Layered system honoured dual-cell BAG, derived from halo protected primal mesh. (b) Dual-cell BAG, derived from halo protected primal mesh (when viewed in primal setting will become primal-cell-WAG).

Fig. 23
A complex feature aligned grid; highlighted are the features honoured by the meshing technique and fracture [5]. As a result of the different fluid storage and conductivity characteristics of the matrix and fractures, boundary-aligned grids are required to honor geological boundaries within the finite-volume approximation. In the following, to demonstrate the capabilities of the grid generation method when employed to grid a fractured reservoir, a synthetic reservoir model embedded with a fractured network is meshed. Figure 24a, b display grids honoring a fracture network with respect to primal and dual-cells, respectively.

Extension of proposed grid generation to field applications
Reservoir geometries are often comprised of a layered structure and extend horizontally in the order of kilometers, whereas vertically they can be only a few feet. By incorporating geological information of a field reservoir, a discrete model corner-point-grid (CPG) [65] comprised of hexahedra is constructed, which honours geological features involving faults, fractures, pinch-outs and layers. By  default, a CPG grid does not honor well-paths [31]. Data sets defining well-paths are comprised of well-perforations and well-trajectories. 2.5D grid generation methods rely on the 2D Delaunay triangulation technique [9]. In order to start with 2.5D grid generation, a 2D layer is extracted from the input 3D CPG model. Considering an extracted 2D layer as a reference layer, 3D wells are projected onto it. By employing methods proposed to honor features in 2D (described above), a well-path honored mesh of a reference layer is generated, which is then projected in the z-direction generating a 2.5D well-path honoured geological layer aligned mesh. The key steps involved are shown in Fig. 25.
Populating an unstructured grid model Reservoir simulation is commonly performed on structured grids, for which there exist a wide variety of property modelling geostatistical algorithms [15]. There are two approaches to populate unstructured models. Either populate the unstructured computational model using geostatistical algorithms [15] or map geological properties from a structured model [62]. The later approach is commonly employed and involves building geostatistical reservoir models on a fine grid and then average and/or upscale [79] properties to the coarser unstructured grid. Upscaling or homogenization involves assigning heterogeneous properties of a region of structured grid cells to an equivalent homogeneous region comprised of a single coarse-grid cell with an effective property value.
For volumetric properties such as porosity and saturation the effective property representative of a set of geo-cells can be established simply by bulk and pore volume weighted average respectively [81]. Homogenization of intrinsic properties, e.g., absolute permeability does not have an exact analytical solution except for in a few idealized cases, and consequently, there exists a large number of upscaling techniques [79,81]. The quality of effective permeability approximation depends on the complexity of fine-scale permeability distribution and the upscaling method used [81]. Mapping rock properties from geo-cellular grids to unstructured grids is performed in a systematic fashion, honouring layer boundaries, unconformities, and faults.

Pressure equation
The pressure equation arises from mass conservation together with Darcy's law and is written in integral form as: where φ represents field pressure; ∇ is the gradient operator, K is the elliptic symmetric permeability tensor; q is the source term, which is zero away from well sources or sinks.
(a) Well-path aligned primal mesh of a reference CPG layer and associated dual mesh (b) Projecting dual mesh in z-direction The first step of the finite-volume formulation involves use of the Gauss divergence theorem to integrate Eq. 7, over a control-volume . After integration Eq. 7 is then written as where δ corresponds to the boundary of control-volume , d is an element of control-volume surface area, i is the i th face of the control-volume and n f is the number of faces; n i is the outward unit normal to face i as shown in Fig. 26. The resolution of Darcy velocity −K∇φ along the unit normal n i is called the Darcy-flux through face i. Approximation of Darcy-flux is a key step in a finite-volume formulation and many approximations have been proposed. We briefly summarize properties of the finite-volume methods used in this work in the review of CVD-MPFA methods below and discuss aspects of the cell-centred formulation in Appendix D. general structured and unstructured grids. The CVD-MPFA formulations overcome the limitation of the standard reservoir simulation two-point flux approximation (TPFA) which is only consistent on K-orthogonal grids [1,19], while providing a consistent generalization of the TPFA method that depends on exactly the same optimal number of degrees of freedom. Auxiliary interface pressures are introduced on control-volume faces in the CVD-MPFA formulation to locally impose normal flux and pressure continuity conditions across control-volume interfaces. Prior to the solution process, the auxiliary pressures are expressed algebraically in terms of the primal pressures via the local flux continuity conditions and lead to a locally conservative formulation with continuous fluxes only dependent on primal pressures. Both cell-centred and vertex-centred CVD-MPFA approximations are considered in this work [19,20,[22][23][24][25][26][27][28]. A comparison between cell-centred and vertex-centred CVD-MPFA schemes is presented here, in terms of computed flow fields resulting from the respective pressure equation approximations. A comparison between the CVD-MPFA schemes and the reservoir simulation standard two-point flux approximation (TPFA) scheme is also presented. The cell-centred CVD-MPFA formulation and standard TPFA reservoir simulation scheme are summarized in Appendix D.
The cell-centred and cell-vertex CVD-MPFA formulations involve multiple families of schemes defined by the local flux quadrature point parameterization on each control-volume face [25,28]. A single family is parameterized by a dimensionless variable q.
For a given grid type, there are two basic types of CVD-MPFA formulation determined by the choice of basis functions: (a) triangle pressure support (TPS) with linear basis functions over subcell triangles leading to pointwise pressure continuity on control-volume sub-faces [19,20,23,24]; (b) full pressure support (FPS) with subcell bilinear basis functions, leading to full pressure continuity over control-volume sub-faces [22,25,28].
The TPS formulation can yield unstable results with strong spurious oscillations for challenging test cases involving strongly anisotropic full-tensor fields, with nonlinear pressure fields resulting from, e.g., injector/producer wells (or source/sinks). In contrast, FPS methods overcome the TPS limitation and yield well-resolved pressure fields that are free or essentially free of strong spurious oscillations [22,25,27,28], where M-matrix conditions are also presented that subject to their satisfaction, ensure all such linear methods have a local discrete maximum principle (LDMP), hence conditional stability. While both the TPS and FPS formulations are shown to violate M-matrix conditions in such cases (i.e., for strong full-tensor field cases where TPS induces spurious oscillations) and therefore lack a local discrete maximum principle required for stability [22,28], it is also shown that the TPS formulation permits decoupled solution modes of grid level resolution for such cases, and conversely that the FPS formulation prevents such modes and thus prevents their accompanying spurious oscillations, allowing the full multi-family quadrature range for flux approximation [22,25,28]. We also note that the unstructured cell-vertex TPS formulation is an exception that does not suffer from grid level decoupling and associated spurious modes [26], and the results of the work presented verify this property for general unstructured grid examples.
Measures of M-matrix violation are presented in comparative tables including relative measures of diagonal dominance, for the cell-centred and cell-vertex TPS formulations on challenging unstructured grid examples. The tables are seen to correspond with method performance.

Cell-centred versus cell-vertex CVD-MPFA
Traditional reservoir simulation methods are primarily block centred, i.e., cell-centred methods; however, both cell-centred and vertex-centred frameworks continue to be developed. While some comparative studies have been undertaken [16,37], a definitive conclusion has yet to be reached. For structured meshes, the number of primal and dual cells are basically equivalent with an off-set for boundaries. On unstructured grids, the number of cells (or elements) is essentially double the number of vertices in 2-D. Consequently for a given grid, the cell-centred formulation requires approximately twice as many degrees of freedom compared to the vertex centred formulation in 2-D. This is easily illustrated by constructing an unstructured mesh from an underlying structured mesh, in 2-D this requires subdividing each quadrilateral into two triangles, and in 3-D each hexahedron can be split into 5 or 6 tetrahedrons. Consequently the cell-centred formulation can be expected to involve many more degrees of freedom for a given grid of nodes, and will be more computationally expensive, but might be expected to resolve flow fields more accurately.
A comparative study of cell-centred versus cell-vertex control-volume distributed multi-point flux approximation (CVD-MPFA) schemes performance, applied to problems involving variable and distorted geometric boundaries and layers on essentially equivalent primal meshes, has been made possible by the grid generation methods presented here. Previous studies of cell-vertex CVD-MPFA methods performance have only been possible for simple interior geometry due to boundary aligned dual-grid generation limitations. Meshes used for the vertex based formulations employ the median dual(see Fig. 28a) of meshes used in the corresponding cell-centred formulation, with the local exception of internal boundaries. Preservation of internal geological boundaries and features is achieved by the use

Duality of cell-centred and vertex-centred methods
The formulation of CVD-MPFA cell-centred or vertex-centred methods follow essentially analogous principles. The most fundamental being the switching of control-volume from primal to dual or vice versa, which is achieved by switching the control point from primal-cell center to primal-cell vertex or vice versa, together with appropriate boundary conditions. However, on unstructured grids, the vertex-centred method is considerably more efficient and robust, with a relatively simpler data structure compared to the cell-centred method. Figure 27 illustrates important differences between cell-centred and vertex-centred CVD-MPFA schemes. We note that: -In cell-centred CVD-MPFA formulations, primal-cells (triangles and/or quads) act as control-volumes. Interface pressures are eliminated via continuity conditions formed with respect to clusters surrounding vertices in dual-cells. In contrast in vertex-centred formulations, interface pressure elimination takes place in each primal cell. -Cell-centred formulations have larger stencil size compared to vertex-centred formulations, e.g. Fig. 27 where the resulting matrix bandwidth of the cell-centred formulation is larger than that of the vertex-centred CVD-MPFA formulation, here 13 versus 7.
Cell-centred and vertex-centred CVD-MPFA methods require respective feature-protected primal and dual meshes as input. In cell-centred methods, primal-cells are aligned with geological features and act as control volumes. For vertex-centred methods, a feature protected dual mesh is generated such that control-volume boundaries are aligned with geological features as discussed above. If the cellcentred method is designed to cope with general polygonal grid elements, i.e., by adopting the cell-vertex view, then it must inherit the cell-vertex data structure to benefit from an input mesh comprised of voronoi polygons. The cell-centred formulation on such a polygonal (polyhedral in 3-D) mesh would then benefit from the advantages of a vertex-centred formulation, of course in this case such a method would effectively have the cell-vertex framework and data structure implanted. The key advantages of employing vertex-centred methods on dual control-volume meshes over cell-centred methods on primal-meshes are listed below: -For a given primal mesh, there are far fewer vertices than cells, approximately a half in 2D. Thus for a given primal grid a cell-vertex formulation is optimal with respect to degrees of freedom compared to a cell-centred method. -Vertex-centred CVD-MPFA schemes are compact on a given unstructured grid. The vertex centred assembled CVD-MPFA formulation matrix can have a much reduced bandwidth compared to the cell-centred CVD-MPFA formulation, c.f. Fig. 27, which is due to the cell-centred flux assembly being performed over each cluster of triangles attached to each grid vertex. -Vertex-centred CVD-MPFA schemes employ a simpler data structure when compared to cell centred methods as assembly is performed with respect to primal grid cells. -Vertex-centred CVD-MPFA operating on primalmeshes provides freedom of selecting centroid and/or circumcentre (for an acute triangulation) as the dualpoint allowing median and/or voronoi polygons to be used as control-volumes.
Quadrature selection and implementation of CVD-MPFA Control-volume distributed(CVD) flux approximation schemes are implemented here, both in cell-centred and cell-vertex frameworks following the TPS formulations of [19,20,23,24], and the more robust FPS formulations [22,25,27,28]. The families of CVD-MPFA schemes are defined by the choice of quadrature points. In order to assess the comparative performance of cell-centred versus vertexcentred CVD-MPFA formulation, TPS is used both with default (q = 1.0) and the SPD variant (q = 2/3) for arbitrary triangular meshes [24]. Note that the TPS quadrature range is less stable than the FPS quadrature range [22]. For FPS we use (q = 1.0, c = 0.001, p = 0.5). The quadrature point q can be selected anywhere between the cluster vertex (excluding the cluster vertex (q = 0) which is singular) and edge mid-point(q = 1) with 0 < q ≤ 1, where continuity of flux is imposed. The parameter c controls the size of auxiliary dual cell built around the cluster point to impose the divergence condition and p corresponds to the flux quadrature point for this auxiliary dual cell. A full permeability tensor is used in each of the test cases given below and is defined by for choices of θ .

Comparative performance of CVD-MPFA schemes
The comparative performance of cell-centred versus vertexcentred CVD-MPFA formulations is assessed by application to problems involving strongly anisotropic full-tensor permeability fields, presented in the results section below. Note that the grids employed are essentially equivalent in the primal setting, which is achieved by the methods of grid generation presented above. In the cell-vertex formulation, the dual mesh is constructed from the primal mesh using the median dual unless stated otherwise. M-matrix conditions for CVD-MPFA schemes are given in [19,22,25,27,28].

Case 1: highly anisotropic permeability tensor
This case is taken from [49], where a unit square domain [0, 1] × [0, 1] with a square hole (4/9, 5/9), characterizing a source positioned at the centre is considered. The pressure field resulting from Eq. 7 is simulated, with an anisotropic permeability tensor K, derived from Eq. 9 by setting k 1 = 100, k 2 = 1 and rotated through θ = −30 • relative to the domain frame of reference. Dirichlet boundary conditions with pressure φ = 2 and φ = 0 are prescribed at the walls of the inner and outer square regions respectively. Domain discretization is performed by employing the above triangulation technique, the resulting meshes obtained are shown in Fig. 29. The primal (triangle) cell mesh comprised of 4926 number of elements and has been generated for use with cell centred formulations. Whereas its associated dual(median) mesh consists of 2567 control-volumes, equal to the number of points in the primal mesh, and is designed for use with vertex-centred methods. The flow fields produced by the TPFA, TPS and FPS schemes both in cell-centred and cell-vertex modes are shown in Fig. 30. The top row corresponds to contour plots of the respective pressure fields for the cell-centred simulations, and the bottom row displays respective pressure fields for the cell-vertex simulations. The grid employed is not Korthogonal; thus, TPFA does not yield a consistent solution, nevertheless the solution has a discrete maximum principle(DMP), and is bounded between (0, 2) [49]. The TPFA solution presented has a DMP, both in cell and vertex centred configurations. We also note that using TPFA with a voronoi based mesh, i.e., a PEBI grid, yields similar results to TPFA median dual, e.g., see case-2. The CVD-MPFA schemes with triangular and full-pressure support schemes both in cell and vertex centred settings provide consistent approximations of the flow field. The degrees of freedom for cell-vertex simulations are roughly half(2567/4926) compared to those of the cell centred formulation. We note that the cell-vertex-based methods use less computational time and capture the flow field with at least comparable resolution compared to the cell-centred simulations.

Case 2: strong discontinuous full-tensor with imposed source and sink
The following test case has been selected from [22] (where the test is on a structured grid), and involves simulating a pressure field in a unit square domain subdivided into four quadrants, with a piecewise varying permeability tensor. The permeability tensor is defined by setting k 1 = 3000, k 2 = 1 and θ = 25 • in Eq. 9 which is assigned to first and third quadrants respectively, whereas for the second and fourth quadrants the permeability tensor has the same   In cell-centred formulations, the control-point, in general defined as the control-volume centre, can either be selected at the centroid and/or circumcenter of an underlying primal mesh. In vertex-centred formulations control volumes are defined by dual-cells, which are constructed by joining (i) cell centroids or (ii) circumcenters (dual-points) of cells surrounding a vertex. In case (i) a dual cell is constructed by joining cell-centroids to cell edge mid-points of cells surrounding a vertex, leading to a median dual, whereas in case (ii) the dual is obtained by joining circumcentres and is a voronoi(PEBI) dual e.g. see Fig. 28. In general, the control-volume centroid is used as a control (dual) point, this is because an acute BAG can not be guaranteed, unless the geometry is simple and/or a quality mesh generation algorithm in conjunction with mesh cosmetics is employed.
In this case, we study cell-centred versus vertex-centred formulations both with centroid and circumcentre as control (dual) points. The primal-cell acute triangulation shown in Fig. 31a is used as a mesh to simulate the flow field in a cellcentred framework. The comparable median and voronoi dual meshes employed for the counter part cell-vertex formulations are constructed from the same underlying acute  Fig. 31b. The resulting median and voronoi dual meshes are shown in Fig. 31c, d, respectively. The resulting grid generated has 4024 primal-cells (triangles) and 2140 dual cells, the mesh generated respects the interior boundaries in both primal and dual settings. In order to capture flow field resolution more closely, the grid has been generated with increased resolution around points where the source and sink are specified.
Centroid as control and dual point The pressure field simulated by employing TPFA and CVD-MPFA schemes both in a cell-centred and cell-vertex framework with centroid as control point and dual point is shown in Fig. 32 (the first and second rows correspond to the cell-centred and cell-vertex formulations respectively). From Fig. 32a (middle), it is observed that the cell-centred TPS formulation introduces spurious oscillations, due to decoupling [25] and violates the discrete maximum(DMP) principle, the solution presented is in accordance with that computed in [22]. These results provide further evidence that the cell-centred TPS formulation, when used to simulate a flow field governed by a strong anisotropic full permeability tensor together with imposed source/sink, introduces spurious oscillations due to decoupling. In contrast, when the same problem is simulated by employing the cell-vertex TPS formulation on an unstructured grid, the method yields a well resolved  solution, see Fig. 32b (middle). We note that cell-vertex TPS on triangles does not generally lead to decoupling [26,27]. Note these linear methods violate the M-matrix conditions for such cases involving strong full-tensor coefficients and therefore do not satisfy a local discrete maximum principle(LDMP) needed to ensure stability [22,25,28]. The degree of M-matrix violation (the measures are defined above) of cell-centred versus vertex-centred TPS is  Fig. 32c. From the results and the table we conclude that the cell-vertex TPS formulation is far more robust when compared to its cell-centred counterpart on unstructured grids. The FPS formulation also violates the M-matrix conditions for this case, but does not suffer from decoupling and the results verify that FPS yields well resolved solutions for both cell-centred and cell-vertex approximations.
Circumcenter as control and dual point The pressure fields computed with circumcentres as control points (pressure approximation points) in the cell-centred formulation, and circumcentres as dual-points in vertex-centred (voronoi dual) formulations are compared employing both TPFA and CVD-MPFA schemes. The results are displayed in Fig. 33, (a different mode of spurious oscillation is noted for cell-centred TPS (using circumcentres) compared to cell-centroid TPS). While cell-centred TPFA detects more cross-flow than the cell-vertex TPFA method, TPFA is generally inconsistent for this case. The cell-vertex TPS formulations with circumcentre and centroid as control point respectively yield similar stable results, c.f. see Figs. 32 versus 33 and are comparable with FPS results, where FPS is robust for all grid types.

Case 3: strong discontinuous full tensor zigzag field with an embedded well trajectory
This case is motivated by [22] where the test is on a structured grid, where a unit square domain, is divided into three distinct regions with an embedded well trajectory defined along a vertical line across the centre of domain is considered. Dirichlet boundary conditions are prescribed by setting pressure to zero along the boundaries of the square domain. Definition of the zigzag permeability field involves a full-tensor defined in the bottom and top regions by Eq. 9 with k 1 = 3000, k 2 = 1 and θ = 25 • , whereas the middle region has the same k i , with θ = −25 • . Again the grids are generated by employing the above triangulation technique and shown in Fig. 34 both for cell-centred and cell-vertex formulations. The grid generated is boundary and well aligned, both in primal and dual forms. The cell-centred grid has 4412 control volumes(primal cells), whereas the vertex centred grid is comprised of 2329 control volumes(dual cells). Figure 35 shows the numerical pressure solutions computed by the TPFA and CVD-MPFA schemes. Again the cell-centred TPS solution is found to violate the M-matrix conditions [25] and introduces non-physical oscillations consistent with decoupling, while the cell-vertex TPS (as noted previously) is more robust, see Fig. 35a versus b. The FPS-based formulations are found to be consistently far more robust yielding well-resolved solutions for both grid types c.f. right hand figures of Fig. 35a, b, respectively. Measures of M-matrix violation of cell-centred versus vertex-centred TPS schemes are tabulated in Fig. 35c.

SPD variant of cell-centred CVD-MPFA-TPS
Triangular pressure support(TPS) schemes with quadrature q = 2/3 on triangular meshes yields a symmetric positive definite(SPD) variant of the family of schemes [24]. A study using this scheme is also presented in [44]. Here, we compute the pressure field of case-3, using the SPD variant of the TPS family. The resulting flow field for case-3, is shown in Fig. 36. Comparing the result listed in Fig. 36, to that displayed in Fig. 35a (middle), indicates that the SPD variant of the TPS cell-centred approximations, yields an improved flow field relative to the default TPS(q = 1) cell-centred formulation. However, both of the TPS approximations induce significant spurious oscillations.

Case-4: Full tensor zigzag field embedded with a complex well trajectory and deviated layered system
This case contains more challenging deviated internal boundaries than the previous cases and serves to further test the grid generation methods presented and increase the challenge in comparison between the CVD-MPFA cell and vertex centred formulations. A unit square heterogeneous domain contains an embedded layered system, which partitions the computational domain into four distinct regions, Fig. 20 displays a pictorial representation of the test case. A piecewise constant permeability tensor is assumed in each sub-domain with (principal axes) anisotropic ratio of 1000 : 1, and its orientation (θ = ±25) is varied so as to define a zigzag flow field. A vertical well trajectory is located in the middle of domain, that intersects a layer and bifurcates into a multilateral well in the neighbouring region. The well trajectory is considered as a geometrical object, and has prescribed pressure φ = 1. We impose homogeneous Dirichlet boundary conditions(φ = 0) on the outer walls of the computational domain. This problem presents significant challenges to both grid generation methods and numerical schemes.
For cases where layers and well trajectories are nonintersecting, grid generation is relatively easy. The model case at hand requires both boundary and well trajectory aligned meshes. The generation of grids honoring geological features with respect to control volumes (BAG) involves entirely different strategies compared to those requiring control point alignment(WAG). This is further exacerbated, when these intractable features appear while intersecting each other, since they meet conflicting requirements at the point of intersection. For such complex geometries, the proposed feature-based triangulation technique proves versatile. We employ non-uniform triangulation and generate boundary and well-aligned grids, with refinement in the near well regions, using the new meshing technique. The resulting primal and dual-cell boundary and well-aligned meshes thus obtained are shown in Fig. 21a, b, respectively.
The problem poses serious challenges to the numerical schemes, which is mainly due to large anisotropic ratio and local non-aligned orientation of the grid. Figure 37 shows the numerical pressure fields obtained by employing TPFA, and the CVD-MPFA schemes, both in cell-centred(top row) and vertex-centred(bottom row) frameworks. The TPFA schemes are inconsistent and cannot resolve the anisotropy of the problem. Again all of the CVD-MPFA methods violate M-matrix conditions for such cases with strong full-tensors. However, the cell-centred TPS formulation is found to have the maximum M-matrix violation and introduces non-physical oscillations due to grid level decoupling, c.f. Fig. 37a (middle), whereas the vertex-centred TPS(on dual mesh) Fig. 37b (middle) and both cell and vertexcentred FPS (Fig. 37a (right), b (right) respectively), yield consistent solutions that are free of the corresponding spurious modes. We conclude that (i) the cell-vertex unstructured grid TPS formulation does not suffer from decoupling (consistent with [26]) that is inherent in the cell-centred TPS formulation [25], (ii) by design the FPS formulations have   [30,53] time (N represents number of mesh points), which is primarily due to two factors: The proposed grid generation methods involve two novel components, i.e., halo construction and protection circles (Section 4). Halo construction and protection-circles are built around key geological features before field triangulation starts. For field meshing advancing front point placement is used, and field points are introduced during the first level of insertion and filtered by protection circles (e.g., Section 4). Additional computational cost induced by halo construction and protection circles depends on the number of geological features, but is relatively small compared to the time incurred for complete grid generation. Grids generated are employed to simulate pressure fields. When compared to simulation time (dominated by solver), the computational cost of grid generation is far less.

Conclusions
Geological features can be classified into two groups; the first of which deals with domains involving features including layers, fractures and/or faults, and the second treats well-trajectories. In the former, controlvolume aligned grids also known as boundary aligned grids(BAGs) are used, while in the latter, control-point aligned grids are required, which we have called well aligned grids(WAGs). The geological features are honored either with respect to primal-cells, or with respect to dual(voronoi) cells constructed from an underlying simplicial mesh. In terms of features honored, there exists duality between boundary aligned(BAG) grids and well aligned grids(WAG). To honor geological features in primal and dual frameworks, we introduce the ideas of protection circles and quadrilateral halo construction respectively. In order to generate the primal-cell BAG, curve(s) characterizing geological features are embedded in discrete form. To start with triangulation, we generate an empty mesh and then by employing local refinement we honour any missing connections, including interior boundaries. The feature honored empty mesh, generally comprises of low quality simplices, and in order to generate a well resolved mesh comprising of quality simplices, the empty mesh is refined by introducing new(field) points. In order to avoid any swapping which may occur with Delaunay triangulation and to preserve integrity of the features honoured in the empty mesh, we introduce the idea of protection circles (circum and/or diametric circles) around the feature based simplices. In field mesh generation, integrity of features is honored provided new points falling in the protection circles are not accepted.
Dual-cell BAGs are derived from halo protected primal meshes. Beginning with an empty mesh where features are embedded in discretized form, before field mesh generation takes place, halos comprised of quadrilateral cells are constructed such that honoured features can be retrieved as medial lines of halo elements. The proposed triangulation method has been successfully employed to generate boundary and well aligned grids, which honour geological features both in primal and dual settings.
The cell-centred and cell-vertex CVD-MPFA formulations involve analogous steps in flux approximation, and requires switching control-volumes from primal to dual or vice versa together with appropriate data structures and boundary conditions. The cell-centred formulation requires primal cells(triangles and/or quadrilaterals) be boundary aligned with the specified geological features, whereas the vertex-centred formulation requires that the dual-grid, comprised of control-volumes surrounding grid vertices, be boundary aligned with specified geological features, with boundary alignment with respect to the control-volume polygonal element faces. Generation of comparable boundary aligned cell-centred and cell-vertex grids in terms of cell resolution has enabled a comparative study of cell-centred versus vertex-centred CVD-MPFA formulations on essentially equivalent primal meshes. As a result, this work provides the first comparison of the two types of CVD-MPFA approximations on comparable unstructured grids involving challenging geometries and challenging permeability fields. The comparative performance of TPS is summarized in Table 1 below.
Key observations on approximation schemes: -TPFA results verify that the standard reservoir simulation flux approximation (with two-point pressure difference) is inconsistent on non K-orthogonal grids. -Cell-vertex simulation requires much less computational time compared to the cell-centred formulation. This is because the cell-vertex formulation involves approximately half the number of degrees of freedom compared to the cell-centred formulation when using the same primal unstructured grid in two-dimensions, and is thus computationally more efficient. -Cell-vertex CVD-MPFA methods are compact and result in considerable reduction in assembled matrix bandwidth compared to cell-centred CVD-MPFA on a given (primal) unstructured grid making a further contribution to computational efficiency. -The cell-centred TPS-based CVD-MPFA formulation yields results with strong spurious oscillations due to decoupling at the grid level, when applied to problems involving strong full-tensor fields and non-linear solutions. Some improvement is gained by using the SPD variant of the scheme or using the circumcentre as the approximation point, however the results still exhibit significant spurious oscillations. -Cell-vertex CVD-MPFA formulations (point-wise continuous TPS) are computationally more robust than their cell-centred counter part formulations on unstructured grids and free from grid level decoupling, yielding consistent well resolved solutions with approximately half the number of degrees of freedom. While the cellcentred TPS formulation is effective for full-tensor test cases with weaker off-diagonal coefficients, pressure fields with strong spurious oscillations are obtained for some strong full-tensor test cases, even though the cellcentred method uses twice the number of degrees of freedom.  [84], TetGen [83], and Gmsh [82], which can be used to generate primal-cell boundary aligned grids with limited control over mesh quality and constraints (features) being honoured. To the best of our knowledge there is no open source mesh generation tool which can be deployed in its current form to generate dualcell feature honoured grids that produce the grids generated here.
For example in 2D using the open source tools for each feature to be honoured in the dual settings, it would require embedding two PSLGs enclosing the actual feature as its medial line giving rise to a channel. Triangulation would then be performed with a prescribed characteristic length to ensure that no element forms inside the channel, and would be accompanied by a boundary recovery procedure. Assuming that one succeeds in generating such a channel protected mesh, then if the feature to be honoured is a straight line, the dual mesh derived from the triangulated channel protected primal mesh would honour the actual feature approximately with undulation joining centroids of the triangulated channel. Generating dual-cell feature protected mesh is an over constrained problem for an open source tool. However, the meshing performed here is able to guarantee medial alignment and allow primal nodal path alignment crossing the medial alignment. Despite several advances in flux approximation schemes, the two point flux approximation (TPFA) is still regarded as the standard reservoir simulation technique, and uses a two point pressure difference in the flux approximation across the face of a control volume [7]. Figure 38 displays two control volumes A and B, with ab being the interface between them, φ a and φ b represent field pressures at their control points, respectively. The reservoir simulation standard two point flux approximation is given by:

Appendix C: Meshing algorithms
where k a = n • K a n, k b = n • K b n, with n being unit normal to the interface ab ; d a and d b represent the respective distances of control points in control volumes A and B from the interface ab . The expression for the two point flux is obtained by satisfying flux and pressure continuity on a K-orthogonal grid e.g. [64]. While TPFA is locally conservative and always has a discrete maximum principle, this standard approximation is only consistent if the grid is Korthogonal otherwise an O(1) error in flux is introduced by the TPFA flux approximation, e.g. [19].

D.2 Formulation of CVD-MPFA schemes
The main advantage of CVD-MPFA schemes is that they provide a consistent locally conservative flux-continuous  approximation of the pressure equation for any permeability tensor and grid type, while only depending on a single degree of freedom per control-volume. The continuity conditions imposed around every cluster point [24] (e.g. Figs. 40a and 41a), leads to an increased pressure support with wider stencil compared to the standard TPFA scheme as shown in Fig. 39a, b. However, all other methods that rival CVD-MPFA in terms of consistency and flux continuity depend on a much larger number of degrees of freedom and consequently yield much larger matrices. For example mixed finite element methods e.g. [4,68] require three times as many degrees of freedom in 2-D, that include the standard cell-centred pressure in the approximation together with edge values, which increases to a factor of four in 3-D (on structured grids). The mixed hybrid finite element method (MHFEM) [10] and mimetic methods [51] only depend on control-volume face values and have an SPD matrix. While reducing the degrees of freedom compared to the original mixed methods, with the traditional control-volume centred pressures now removed, MHFEM still involves twice the degrees of freedom when compared to the CVD-MPFA formulations in 2-D, while CVD-MPFA depends entirely upon the traditional control-volume centred pressures. CVD-MPFA cell-vertex and cell centred schemes basically divide into two types, namely triangular pressure support(TPS) [19,20,23,24] and full pressure support(FPS) [22,[25][26][27][28] schemes and are discussed further below. Next, we briefly describe the procedure to construct TPS and FPS flux approximations across control-volume faces in a cell-centred setting. The gradient ∇φ(e.g. Eq. 7), is expressed in a (ξ, η) coordinate system (Figs. 40b and 41b) as: where J −1 is the inverse of Jacobian matrix J defined by: In cell-centred formulations primal cells, i.e., triangles and/or quadrilaterals in 2-D act as control volumes. Each control-volume is subdivided to sub-cells(sub-quad) equal to the number of vertices of the underlying control-volume. On each sub-cell or sub-quad, TPS has triangular pressure support, and FPS has full (quadrilateral) pressure support.

Triangular pressure support (TPS) schemes
Triangular pressure support(TPS) schemes are point-wise continuous both in normal pressure and normal flux, and involve the use of triangle pressure support basis functions, e.g., [23,24] for further details.
The resulting ξ and η partial derivative approximations are given by: with analogous approximations for position vector. The discrete Darcy velocity −K∇φ is then resolved along the normal to each primal cell half-edge (each adjoining pair

Flux computation and notion of general tensor in TPS scheme
In the cell-centred formulation, flux approximation is performed across the half edges of the triangles/quadrilaterals cells (control volumes), e.g., referring to Fig. 40a p 1 m 1 and p 1 m 3 form a pair of half-edges. The computation of flux is defined by resolving −K∇φ along the normal (with magnitude equal to half edge length) to these half edges. Let n i be the normal to a representative i th half edge of a sub-cell quad, Darcy velocity is given by: and the flux that results from resolving the Darcy velocity along the normal vector n i is expressed as where i is the local index of half-edge across which −K∇φ has been resolved and T ij corresponds to general tensor component. Since each sub-quad shares two half edges of underlying control volume (p 1 m 1 and p 1 m 3 Fig. 40a), the local indices of these half edges varies i = 1, 2, defining the respective fluxes over two half edges of each sub-quad with n 1 and n 2 being normal to these edges. Note that in d dimensions T is a d ×d tensor, a symmetric scheme is given in [24].

Construction of local system
TPS is point-wise continuous in pressure and normal flux. In order to ensure local flux continuity continuous auxiliary interface pressures are introduced on the interfaces which are then eliminated by enforcing flux continuity, locally around each vertex of a control-volume. The TPS flux is defined at the pressure continuity point which is selected anywhere between the cluster point and associated edge mid-point and is parametrized by quadrature q, e.g., see Fig. 40a. Figure 40c illustrates a three-cell cluster around point p, for each sub-quad edge associated with the cluster point p there are two normal fluxes, one either side. From Eq. 14, each such flux involves three pressures, i.e., one cell-centred pressure and two edge pressures(auxiliary interface pressures). Flux continuity through the half edges of each cluster leads to a local system, expressed as: (15) where φ ν = [φ 1 , φ 2 , φ 3 ] t and φ f = [φ a , φ b , φ c ] t represent cell-centre and interface pressures respectively. Solving the above system for interface pressures yields: where  (17) where flux vector F = [F pa , F pb , F pc ] t , e.g. Fig. 40c.

Quadrature parametrization and family of TPS schemes
Referring to Fig. 40a φ 1 is the pressure at cell-centre (control point), whereas φ a and φ b are the interface pressures with default position at the edge mid points, however, their position can be parametrized anywhere between cell-vertex p 1 and edge mid point m 3 and m 1 with parameters q 1 and q 2 respectively. This parametrization naturally leads to families of TPS schemes [23,28]. Here we choose a symmetric variation of quadrature i.e. q 1 = q 2 = q where 0 < q ≤ 1, with q = 1 corresponds to edge mid point and q = 0 to the excluded cell vertex. Note q = 2/3 corresponds to a symmetric positive definite approximation [24]. Figure 40a displays pressure triangles constructed with different q, note that q = 0 degenerates the pressure triangle to a line segment and is avoided.

Full pressure support (FPS) schemes
When applied to problems involving strongly anisotropic full-tensor fields non-aligned with the local grid, together with source/sink boundary conditions (which result in nonlinear solutions), TPS schemes can violate the discrete maximum principle(DMP) and introduce spurious oscillations in such cases due to decoupling, e.g., [22,25]. Full pressure support(FPS) CVD-MPFA schemes involve imposing full pressure continuity over each face of a control-volume together with continuity of normal flux across each face [25]. In addition to continuous interface pressures, FPS involves introduction of an auxiliary pressure at each cluster point (Fig. 41a), which enables full control-volume face pressure continuity to be imposed. The additional auxiliary pressure support of the FPS formulation leads to a wider and more stable range of quadrature compared to earlier TPS in the presence of high anisotropy ratios enabling the decoupling suffered by TPS to be avoided [28]. The auxiliary pressures are then eliminated by employing flux continuity conditions in conjunction with a divergence condition derived from Eq. 7 over a shrinking volume surrounding the cluster point. An alternative type of method to CVD-MPFA schemes which is designed to prevent spurious oscillations is presented in, e.g., [49,50]; however, these methods are non-linear and involve iteration. In contrast, the FPS formulation is linear and leads to a quasi-M-Matrix [22,28], and has been shown to be far more robust than the earlier linear TPS methods [22,25,28] and does not require iteration. The FPS formulation employs a bilinear variation of pressure over each quadrilateral sub-cell belonging to a cluster, parametrized by (ξ, η). The nomenclature, pressure support and isoparametric mapping involved in formulating FPS are shown in Fig. 41. with a sub-cell quad displayed in Fig. 41b, pressure and position vector have a bilinear variation. φ = (1−ξ)(1−η)φ 1 +ξ(1−η)φ a +ξηφ p +(1−ξ)ηφ b (18) The position vector has an analogous approximation. The local ξ and η derivative approximations are given by: Pressure gradient ∇φ, is expressed in terms of (ξ, η) in Eq. 12. Figure 41b shows quadrature parametrization used in FPS scheme. Here, we employ a symmetric quadrature i.e. q 1 = q 2 = q ∈]0, 1] and p 1 = p 2 = p ∈]0, 1] to define the respective sub-cell and sub-sub-cell quadrature points, where c ∈]0, 1] relates to the size of a sub-sub-cell belonging to the auxiliary cluster vertex control-volume for the auxiliary divergence condition. The quadrature parametrization listed in Fig. 41c assumes that q 1 = q 2 = 1 corresponds to the cell-edge mid-points and q 1 = q 2 = 0 to the cluster point (when the scheme is cell centred the cluster point is cell-vertex, and vice-versa), q = 0 degenerates continuity points to cluster points and is excluded (note anisotropic quadrature is considered in [28]). The quadrature for the edges of sub-sub-cells(used to enforce the divergence free condition) may differ from that of the sub-cell, i.e., in general q = p. Note unlike TPS, FPS has -An additional interface pressure introduced at each cluster point. -Additional divergence free condition to eliminate additional auxiliary pressure introduced at the cluster point.
-Interface pressures are attached to edge mid-points and the respective cluster point, while flux quadrature points define families of schemes [22,25,28].
FPS fluxes are computed by resolving −K∇φ along the local control-volume face outward normal, analogous to TPS, but now using a bilinear sub-cell approximation. The fluxes computed by resolving −K∇φ along the normal to the half edges are used to enforce flux continuity. The fluxes calculated across the faces of the sub-sub-cell are used to construct the auxiliary divergence free condition.
Constructing the local FPS system The local FPS system for a cluster involving three interface continuity conditions and a zero divergence condition for the cluster point is illustrated in Fig. 42b and is written as: (20) where φ ν = [φ 1 , φ 2 , φ 3 ] t is the vector of primary cellcentred pressures in the illustrated cluster, and φ f = [φ a , φ b , φ c , φ p ] t is the vector of interface pressures including the auxiliary pressure at the cell-vertex. From the above system, the interface pressures are expressed in terms of the primary pressure variables as: where B = [B L − B R ] −1 and A = [A R − A L ]. Back substitution into Eq. 20 and after discarding the last row defining the divergence free condition, the expression for fluxes through the half edges is expressed in terms of cell-centre pressures viz: where flux vector F = [F pa , F pb , F pc ] t , e.g. Fig. 42b.