Extending the global-direction stencil with “face-area-weighted centroid” to unstructured finite volume discretization from integral form

Accuracy of unstructured finite volume discretization is greatly influenced by the gradient reconstruction. For the commonly used k-exact reconstruction method, the cell centroid is always chosen as the reference point to formulate the reconstructed function. But in some practical problems, such as the boundary layer, cells in this area are always set with high aspect ratio to improve the local field resolution, and if geometric centroid is still utilized for the spatial discretization, the severe grid skewness cannot be avoided, which is adverse to the numerical performance of unstructured finite volume solver. In previous work [Kong, et al. Chin Phys B 29(10):100203, 2020], we explored a novel global-direction stencil and combined it with the face-area-weighted centroid on unstructured finite volume methods from differential form to realize the skewness reduction and a better reflection of flow anisotropy. Greatly inspired by the differential form, in this research, we demonstrate that it is also feasible to extend this novel method to the unstructured finite volume discretization from integral form on both second and third-order finite volume solver. Numerical examples governed by linear convective, Euler and Laplacian equations are utilized to examine the correctness as well as effectiveness of this extension. Compared with traditional vertex-neighbor and face-neighbor stencils based on the geometric centroid, the grid skewness is almost eliminated and computational accuracy as well as convergence rate is greatly improved by the global-direction stencil with face-area-weighted centroid. As a result, on unstructured finite volume discretization from integral form, the method also has superiorities on both computational accuracy and convergence rate.


Introduction
Compared with the block-structured grid, unstructured grid is highly automated in grid generation [1][2][3] and convenient for the grid adaptation [4,5], while the computational accuracy and stabilities are hard to be guaranteed [6], especially on high-aspect-ratio triangular grids [7][8][9]. Regarding the current difficulties, scholars are continued trying to improve the discretization algorithms to break through the bottleneck of accuracy loss and stability deterioration on highly anisotropic grids, and realize the unification of automated grid generation and accurate numerical simulation [10,11].
As for the research about unstructured finite volume discretization on high-aspect-ratio triangular grids, Diskin and Thomas [12,13] tested the accuracy of gradient reconstruction, and results show that the accuracy of reconstructed gradient on such grid is determined by both solution and grid, and the existing discretization schemes cannot meet the requirements of grid quality and computational accuracy simultaneously. On this basis, Diskin and Thomas et al. [14,15], systematically compared the numerical performance of inviscid and viscous fluxes on different node-centered and cell-centered unstructured finite volume methods. Research suggests that for high-aspect-ratio triangular grids, there are few discretization schemes able to preserve the simulation accuracy, and after adding disturbance to grid points, the magnitude of computational errors obtained by all schemes is proportional to cell aspect ratios. Similar conclusions are obtained in Ref. [16].
On this basis, we notice that although there are numerous discretization approaches for high-aspect-ratio triangular grids, studies related to the stencil selection are quite few, and within the framework of unstructured finite volume method, both inviscid and viscous fluxes are greatly influenced by the gradient reconstruction [17], where different stencils play a crucial role. Commonly used stencils are face-neighbor and vertexneighbor stencils that consist of cells sharing faces or vertices with the central cell. Apart from these two commonly used stencils, an ingenious stencil augmentation method is proposed in Ref. [18], which is named as the Smart Augmentation (SA) stencil. For each vertex, an extra cell closest to geometric centroid is selected from vertexadjacent cells and appended to the face-neighbor stencil. This method improves the stabilities of face-neighbor stencil to a certain extent. In order to further improve stabilities of the face-neighbor stencil, efficient gradient stencils are proposed by Nishikawa [19], where the F-decreasing augmentation in combination with the symmetric augmentation stencil augments a face-neighbor stencil with extra cells to attempt to increase the symmetry of the stencil. On this basis, extra cells are added to decrease the reciprocal of the Frobenius norm of a scaled least-squares matrix to minimize the lower bound of the magnitude of the gradient. This stencil augmentation approach is much more robust as well as efficient, and the instability of face-neighbor stencil is effectively avoided even on extremely irregular grids.
Apart from the mentioned stencil augmentation methods, Xiong et al. [20], proposed a local-direction stencil selection method, where stencil cells are augmented along two local directions. In this method, characteristics of the flowfield are taken into account during the process of determining local directions. On isotropic triangles, two local directions are close to the normal and tangential directions of the wall, while on high-aspect-ratio triangular grids, stencil cells selected by this method will deviate a lot from the boundary normal, and the numerical performance on such grids will get deteriorated [21,22].
Compared with the local-direction stencil selection method, global direction stencil selection method [22,23] effectively overcomes the deviation on grids with high aspect ratios, and cells are always selected along two global directions, that is, normal and tangential directions of the wall. Hence, the flow anisotropy could be well reflected, and after verification, on high-aspect-ratio triangular grids, global-direction stencil has a better numerical performance on both computational accuracy and efficiency than commonly used face-neighbor and vertex-neighbor stencils as well as the local-direction stencil.
Although global-direction stencil preliminary exhibits a better numerical performance, the distribution of reference points is still much more skewed, especially on high-aspectratio triangular grids. You and Mittal et al. [24], first proposed the grid skewness on such grids, and concluded that grid skewness is adverse to the computational accuracy and stabilities of CFD solvers. Regarding the grid skewness, there are various definitions of that, such as the angle between face normal and the vector pointing form face centroid to the cell centroid, the minimal internal angle of grid cell [25], ratio of the max diagonal to the minimum [26], etc. Nevertheless, from different definitions, the same conclusion could be drawn that on high-aspect-ratio triangular grids, the grid skewness is always quite evident [27].
Different from above skewness measures, Nishikawa [28] proposed a novel definition of that, and it is defined at a face shared by two neighbor cells, say, A and B, as the dot product of the unit vector pointing from the centroid of cell A to that of cell B. Therefore, a non-skewed grid has the measure one and highly-skewed grid nearly zero. Specifically, on highly-skewed triangular grids, the reference point distribution is irregular and exhibits deflective phenomenon. Thus, although global-direction stencil cells mentioned above are along the normal and tangential directions of the wall, the irregular distribution of reference points is not changed, and it is hard to say whether flowfield characteristics are well reflected or not. In order to reduce the grid skewness and optimize the reference point distribution, Nishikawa proposed a novel face-areaweighted centroid [28] for the second-order unstructured finite volume discretization from differential form. After verification, the second-order accuracy is also achieved, and the convergence rate is greatly improved. Besides, compared with the traditional geometric centroid, the distribution of this novel reference point is more regular, and the connection of that is almost parallel to the boundary normal vector.
Based on this phenomenon, in previous work, we combined this novel centroid and global-direction stencil for the second-order unstructured finite volume method [23], and it is verified that the global-direction stencil with face-area-weighted centroid has a lower discretization error than stencils with the geometric centroid including face-neighbor and vertex-neighbor stencils as well as the global-direction stencil. What's more, this novel method has a more stable numerical performance on high-mach-number flow. As a result, the current situation related to accuracy loss and stability deterioration on highaspect-ratio triangular grids is greatly ameliorated by the employment of this improved global-direction stencil. This method is designed for the second-order differential finite volume solver [28], where both solution and source term vectors are evaluated as point values, and therefore, the source term integration is totally avoided. Besides, higherorder accuracy of the differential finite volume solver is also able to be achieved, and professor Nishikawa gave an idea in Ref. [28,29] about the primitive function reconstruction of the flux, by which the governing equation is higher-order accurate without the second-order errors.
Greatly inspired by the face-area-weighted centroid and efficient differential finite volume solver, we explore the novel local origin in the finite-volume discretization based on the integral form and investigate its impact not only on a second-order method but also on a third-order method, since on the integral form, the grid skewness is also existing and needs to be eliminated to guarantee the accuracy and convergence rate.
In this research, the second-order accurate integral finite volume discretization is considered at first, where the k-exact reconstruction [30][31][32][33][34][35][36] based on any local origins is derived, and we combine the global-direction with the face-area-weighted centroid to further reduce the grid skewness and capture the variation of flowfield more accurately. In addition, four representative numerical examples are designed to examine the effectiveness of this extension, and on this basis, the third-order k-exact reconstruction and finite volume discretization based on any local origins are also derived to further extend the global-direction stencil with face-area-weighted centroid to higher-order accurate unstructured finite volume method, and the accuracy as well as the discretization errors is examined by a numerical example with the method of manufactured solutions.
The paper is organized as follows. In Section 2, governing equations and spatial discretization from both differential and integral forms are given at first. And then, we theoretically derive the second-order k-exact reconstruction process based on any local origins. Different stencil selection methods and global-direction stencil with face-areaweighted centroid are presented in Section 3. In Section 4, numerical examples governed by linear convective, Euler and Laplacian equations are employed to verify the effectiveness, feasibilities as well as superiorities of the method on the second-order integral unstructured finite volume solver. The third-order k-exact reconstruction and finite volume discretization based on any local origins, as well as the accuracy analysis are given in Section 5. Concluding remarks and future work will be summarized in Section 6.

Unstructured finite volume method
As mentioned in introduction, the finite volume method from differential form and the discretization base on the face-area-weighted centroid [28] give us great inspirations, and in this paper, to further extend the applicable scope of this novel local origin, we promote it on unstructured finite volume method from integral form, since when unstructured finite volume discretization from integral form is employed, the grid skewness is existing and needs to be eliminated as well.
As a result, in this section, we give unstructured finite volume discretization from integral form at first. In addition, the k-exact reconstruction method based on any local origins is also derived in this section.

Finite volume discretization from integral form
The finite volume discretization from integral form is utilized in this paper, where the governing equation in integral form could be written as where, u j is a solution vector. F c j is the convective flux and s j is a source (or forcing) term vector. Besides, V j and ∂V j are the volume and boundaries of cell j. According to the divergence theorem, Eq. (1) is transformed as Z where, f c n is the convective flux along the face normal vector. On this basis, Eq. (2) could be discretized as follows.
where, u j is the cell-averaged solution vector and {k j } are faces of cell j. Besides, Φ jk is the numerical flux along the normal vector at the k-th face, and A jk is the k-th face area or length (in 2D). The diagram of unstructured finite volume discretization is shown as follows.
As shown in Fig. 1, where u L and u R are two state vectors at the gauss point obtained by owner and neighbor cells of the common face, and C j , C k are two reference points of cell j and cell k respectively. In fact, for traditional unstructured finite volume discretization, these two reference points are located at the geometric centroid of the control volume, while in this paper, the discretization based on any local origins are derived and it will be specifically described in the following section.
For the second-order unstructured finite volume solver, the gauss point is just located at the face centroid and the weight is equal to 1, while for higher-order discretization, there should be at least two points to guarantee the flux quadrature is higher-order accurate. The gauss point location as well as the corresponding weights will be specifically given in Section 5.2. Besides, r j and r k are two vectors pointing from the cell centroid of owner and neighbor stencil cells to gauss point respectively. Two state vectors could be interpolated by the point value u ori j and u ori k at the local origin obtained by the k-exact reconstruction algorithm [30][31][32][33][34][35][36].
Usually, the reference point is chosen at the geometric centroid, while in this paper, the point-valued solution at any local origins is derived, and the main process will be given in the following section. As a result, u L and u R these two state vectors could be formulated by where point value solution of u ori j and u ori k , as well as the solution gradients are all derived and obtained by the k-exact reconstruction algorithm at the local origins C j and C k , which will be specifically described in Section 2.2, and ðx ori j ; y ori j Þ, ðx ori k ; y ori k Þ are local origins coordinates of cell j and face-adjacent cell k in Cartesian-coordinate system. Besides, (x G , y G ) is the gauss point coordinate at per face. On this basis, the numerical flux along the face normal vector could be formulated as where u n = v ⋅ n jk is computed by the Roe flux [37], and H = (γp/(γ -1) + ρv 2 /2)/ρ is the specific total enthalpy. The above content is the second-order accurate unstructured finite volume discretization from integral form, which is employed in this paper and in Section 5, we will specifically exhibit the higher-order, i.e., the third-order accurate discretization process.
As shown in Fig. 2, when calculating the gradient of cell C i , three stencil cells C 1 , C 2 and C 3 are employed to compose the reconstructed equations.
where ω i1 , ω i2 and ω i3 are weights set to emphasize geometrically nearby data.
and, x i and x j are two position vectors from the origin to geometric centroid of cell i and its stencil cell j respectively. As demonstrated in Eq. (6), solution utilized for the least-square reconstruction is the point value at the geometric centroid, and for second-order finite volume methods, this point value is always written in terms of the cell average: u j →u j , but if the reference point is moved anywhere, the point value used for the least-square reconstruction could not be given by the cell average again, and therefore, the reconstruction process could not be implemented, since the point value u ori j is unknown at any local origins in every time step. Besides, least-square reconstruction is second-order accurate and considering the method extension to higher-order accuracy, the k-exact reconstruction method is employed in this manuscript.
Besides, compared with LSQR, the point value at the reference point could be obtained by the k-exact reconstruction, and this point value is utilized for the flux construction. As a result, although this research is focused on the second-order unstructured finite volume discretization, we also employ the k-exact reconstruction method to solve the solution gradient and obtain the point value at the corresponding reference point. The specific process is shown as follows.
In order to obtain the point value at any local origins, the expansion should be formulated based on any local origins as well, where u ori j is the point valued solution recovered at a local origin from cell averaged solutions as we will describe later. Integrating Eq. (8) over the control volume, we obtain where, x ori j and y ori j are the coordinate of any local origins. Replacing ðxx ori j Þ and ðy − y ori j Þ with ðx − x gc j Þ þ ðx gc j − x ori j Þ and ðy − y gc j Þ þ ðy gc j − y ori j Þ respectively, Eq. (9) could be written as which could be further transformed as where, Δx j and Δy j are abbreviations of ðx gc j − x ori j Þ and ðy gc j − y ori j Þ , and if the reconstructed function is expanded based on any local origins, there are another two terms in its integral average, and the mean constraint is rewritten as On this basis, if the reconstructed function is integrated over a stencil cell k, we obtain Similarly replace ðx − x ori j Þ, ðy − y ori j Þ by ðx − x gc k Þ þ ðx gc k − x ori j Þ and ðy − y gc k Þ þ ðy gc k − y ori j Þ respectively, and reconstructed equation on stencil cell k could be written as This equation is written for every stencil cell, of which there should be more than the number of derivatives to be solved to create an overconstrained system. If we write the novel mean constraint together with the Eq. (12) for each stencil cell, we have where the first row is the mean constraint, and geometric terms could be expressed as and the weights are set to emphasize the geometrically adjacent data, where x j and x k are local origins of central and stencil cells. When solving the Eq. (15), the unconstraint reconstructed equations are solved at first, Assuming that the current reconstructed problem is A * x = b, transposition matrix A * T is multiplied at both sides of the equation as A * T A * x = A * T b at first, and then, since the matrix A * T A * is the square matrix, we obtain On this basis, solution gradients as well as point-valued solution u ori j could be obtained. Note, however, that the mean constraint is not satisfied currently, and therefore, the point-valued solution u ori j should be recalculated as Thus, the mean constraint is always exactly satisfied and the obtained u ori j , as well as ∂u ∂x j j and ∂u ∂y j j is utilized for calculating the pointwise solution at per gauss point as shown in Eq. (4). In short, we perform the solution reconstruction at a Gauss point over a face (i.e., obtain u L and u R ) by first computing the point valued solutions and gradients at the local origins with the LSQR method and then use Eq. (4) or a higher-order variant to be presented later. It sets a stage for the employment of novel reference points on unstructured finite volume discretization from integral form without degrading the design accuracy.
3 Global-direction stencil based on the face-area-weighted centroid In Section 2.2, we theoretically derive the feasibility of k-exact reconstruction based on any local origins. In this section, some commonly used stencil selection methods are introduced at first. And then, we briefly discuss the grid skewness on high-aspect-ratio triangular grids as well as the effect of skewness reduction by the employment of facearea-weighted centroid. Finally, to reduce the grid skewness and achieve a better reflection of flow anisotropy, the global-direction stencil with this novel centroid is introduced. The main content in this section has been analyzed in Ref. [23], while considering the completeness of the article, it is also discussed here.

Stencil selection methods
Commonly used stencils are vertex-neighbor and face-neighbor stencils. As Fig. 3 exhibits, face-neighbor stencil includes entire neighbor cells that share faces with the central cell, and vertex-neighbor stencil is similarly constructed by cells that share vertices with the central cell. Besides, both of them are topological-dependent, and therefore are limited by the original mesh topology. Furtherly, the stencil size of them is hard to accurately control, especially for vertex-neighbor stencil, and characteristics of flowfield cannot be well reflected.
Apart from two commonly used stencils, in 2018, Xiong et al. [20], put forward the local-direction stencil selection method, by which selected stencil cells are along two local directions. As shown in Fig. 4a, on isotropic grid, two local directions are close to the normal and tangential directions of the wall, while on high-aspect-ratio triangular grids, as Fig. 4b displays, one of the local directions has severely deviated from the normal direction of the wall, and flow anisotropy is not well reflected. Besides, it is verified that on this grid type, accuracy loss and stability deterioration cannot be avoided [22], and the implementation process of this stencil selection method is quite complicated.
In previous work, based on the existing problems on local-direction stencil, a novel global-direction stencil selection method [22,23] was proposed. Compared with the localdirection stencil, problems mentioned above are well solved by this novel stencil. Specifically, for this method, two global directions, that is normal and tangential directions of the wall, are determined at first. And then, for each central cell, two lines which are parallel to global directions respectively and pass the cell centroid are generated. Finally, cells in a given set, such as the vertex-adjacent cells that intersect with these two lines are selected to construct the new stencil, and the stencil size is governed by layer of vertex-adjacent cells.
As Fig. 5 demonstrates, stencil cells selected by this method are always along two global directions no matter the grid with high aspect ratio or not. Therefore, flow anisotropy can (a) Vertex-neighbor stencil (b) Face-neighbor stencil  [23], global-direction stencil has a better numerical performance than commonly used vertex-neighbor and faceneighbor stencils as well as local-direction stencil. Besides, on cases with simple shapes, two global directions could be easily determined, while for complex surface, there are no analytical expressions to obtain the normal and tangential directions of the wall. Hence, in this situation, we can refer to the method of computing wall distance to get the corresponding normal vector [42], and construct the global-direction stencil. But after analysis, although global-direction stencil cells are always along the normal and tangential directions of the wall, no matter on grid with high aspect ratio or not, the only data obtained by k-exact reconstruction and required for the flux evaluation are solution vectors at the reference point rather than stencil cells themselves, and it is hard to guarantee whether flow anisotropy is well reflected or not. In the next section, we will focus on different locations of reference points.

Face-area-weighted centroid and global-direction stencil
In this section, we give a brief analysis about grid skewness and face-area-weighted centroid on triangular mesh, and explain the reason as well as the idea of the combination of global-direction stencil and this novel reference point in detail.

Grid skewness measure and face-area-weighted centroid
In introduction, we have analyzed that although there are various definitions of grid skewness, the same conclusion could be drawn that on high-aspect-ratio triangles, the grid skewness is always evident. Here, a typical skewness measure is taken into account.
As shown in Fig. 6, the grid skewness is defined at common face shared by two neighbor cells, whereê jk is a unit vector pointing from centroid of cell j to that of cell (a) Minor-aspect-ratio (b) High-aspect-ratio k,n jk is the unit outward normal vector of common face. The grid skewness is measured by jê jk Án jk j, and from Fig. 6, it could be easily found that the non-skewed grid has the measure one, and with the increase of cell aspect ratio, the skewness measure of highly-skewed grid is close to zero. On this basis, the most direct intuition is using isotropic or minor-aspect-ratio triangles to ensure the computing process is carried out on non-skewed grids, while for some typical flows, such as boundary-layer-type flow, solutions in boundary layer are changed dramatically, particularly along the normal direction of the wall. Therefore, to enhance the resolution in this local field, highlyanisotropic grids cannot be avoided, and we can only rely on a novel local origin to reduce the grid skewness and improve the numerical performance.
For high-aspect-ratio triangular grids, a novel face-area-weighted centroid was proposed by Nishikawa [28] on the second-order differential unstructured finite volume solver. By the employment of this reference point, the grid skewness is almost eliminated. In the following analysis, we will briefly summarize the face-area-weighted formula and the effect about skewness reduction. The specific derivation has been given in Ref. [28], and the main results are listed in this section.
A typical choice for the local origin is the geometric centroid that could be written for a triangle by the arithmetic average of face midpoints, where, (x mk , y mk ) is coordinate of the k-th face centroid. Besides, we generalize Eq. (20) to the face-area-weighted centroid formula [28]: where, A jk is the area or length (in 2D) of the face shared by two neighbor cells, and p (> 0) is a real value that controls the skewness degree. Note when p = 0, face-areaweighted centroid is just consistent with the geometric centroid. As shown in Fig. 7, where h is the grid spacing in y-direction and R is cell aspect ratio. For a typical triangular grid in Cartesian-coordinate system, skewness measure along the x-direction approaches one, while in y-direction is nearly zero. Thus, in the Fig. 7 A typical high-aspect-ratio triangular grid in Cartesian-coordinate system following analysis, we will focus on common faces shared by cell 1 and cell 2 as well as cell 2 and cell 3.
For these two common faces, if the geometric centroid is utilized, the grid skewness measure is [28] With the increase of cell aspect ratio, both measures nearly equal to zero. Hence the grid is highly skewed. But if we employ the face-area-weighted centroid and when parameter p is equal to 2, the skewness measure becomes [28] As a result, for high-aspect-ratio triangular grid, the grid skewness could be eliminated, and note, moreover, that with the increase of cell aspect ratio, the grid skewness, at the face like between cell 2 and cell 3, could be further reduced. Besides, it is demonstrated in Fig. 8 when the face-area-weighted centroid is employed, line connecting the novel local origins is almost parallel to the normal direction of the wall, and the serrated phenomenon exhibited on geometric centroid is effectively avoided.
This special distribution just coincides with our original motivation of designing the global-direction stencil, and in previous work, we combined the global-direction stencil and face-area-weighted centroid on the second-order unstructured finite volume solver in differential form to capture the flow anisotropy more accurately. After verification, better computational accuracy and stabilities are obtained by this novel method. In this work, we further investigate feasibilities of extension to the integral form.

Combination of global-direction stencil and face-area-weighted centroid
In Section 3.1, we have demonstrated that compared with commonly used vertexneighbor and face-neighbor stencils, although global-direction stencil cells are along the normal and tangential directions of the wall, the only data obtained by k-exact reconstruction and required for the flux evaluation are solution vectors at cells local origins, rather than stencil cells themselves, and if we still use the geometric centroid, grid skewness is unable to be eliminated. Therefore, it is hard to guarantee whether the flowfield characteristics are well captured or not. But when face-area-weighted centroid is employed, the distribution of reference points is more regular, and the line connecting them is almost along the normal direction of the wall.
(a) Geometric centroid (b) Face-area-weighted centriod Besides, in Section 2.2, we have derived the k-exact reconstruction process based on any local origins, and therefore, in this section, we give the method of combining the global-direction stencil and face-area-weighted centroid in detail to realize the unified direction of both stencil cells and local origins. Global-direction stencils with geometric and face-area-weighted centroids on grids with straight and curved boundaries are displayed in Figs. 9 and 10.
From these two figures, we can easily find although same stencil cells are selected, geometric centroids are deflective and exhibit serrated phenomenon. Particularly on high-aspect-ratio triangular grids, the mentioned phenomenon is much more evident. By comparison, when face-area-weighted centroids are employed, line connecting them is close to the normal direction of the wall no matter on grid with straight or curved boundaries, and it is consistent with one of the global directions. As a result, the flow anisotropy could be well reflected, and the grid skewness is reduced.
On unstructured finite volume method from differential form, both computational accuracy and stabilities are greatly improved by this novel method, and in the next section, four representative numerical examples are designed to verify the effectiveness and superiorities of this novel method on integral unstructured finite volume solver.

Numerical examples
To examine the effectiveness of global-direction stencil with face-area-weighted centroid on unstructured finite volume method from integral form, in this section, numerical examples governed by linear convective, Euler and Laplacian equations are utilized. For comparison, these numerical examples are simulated with four different stencils, including vertex-neighbor and face-neighbor stencils, as well as the global direction stencil with geometric centroid and face-area-weighted centroid. Note that from Section 3.2.1, we can find the face-area-weighted centroid could be further distinguished by different p values, and when p = 2, skewness has been almost eliminated. Hence, the face-area-weighted centroid in this section refers to the result of p = 2. In addition, to simplify the presentation in the following analysis, different stencils are abbreviated as follows (Table 1).
More importantly, in this section, we make a discussion about the result exhibition of unstructured finite volume method from integral form, and notice that the point-(a) Geometric centroid (b) Face-area-weighted centriod valued solution u ori j at the face-area-weighted centroid obtained by the mean constraint of k-exact reconstruction algorithm with global-direction stencil is more accurate than the cell-averaged solution u j . The specific discussion will be given in Section 4.1.3.

Manufactured boundary layer (governed by linear convective equation)
In this section, we first use the Method of Manufactured Solutions (MMS) [43][44][45][46] on linear convective equation to verify feasibilities of the employment of face-areaweighted centroid on unstructured finite volume discretization from integral form, and examine the numerical performance of global-direction stencil with this novel reference point. The model equation can be written as where, α ¼ ðcos π 16 ; sin π 16 Þ is the constant convective velocity, and to simulate characteristics of boundary layer, the manufactured solution is where c = 0.59 is a constant, and parameter μ is utilized to control the thickness of boundary layer. Flowfields corresponding to different μ values are shown in Fig. 11, and in the following test, μ is set as 10 − 6 . By bringing the manufactured solution to Eq. (24), the modified equation with source term could be written as (a) Geometric centroid (b) Face-area-weighted centriod On this basis, the manufactured solution represents the analytical solution of this modified governing equation, and we can calculate both L 2 and L ∞ errors of different stencils.
whereũ j , u analy and A j are numerical and analytical solutions and area of cell j respectively. Note, however, that errors are computed for the point valued solutions recovered at the local origin from the cell averaged solutions stored at cells rather than for the cell-averaged solutions for the reason explained in Section 4.1.3. Therefore, the exact solution used to compute the error is a point value evaluated at the local origin. Besides, as shown in Fig. 12, both regular and randomly perturbed triangular grids are used in this numerical example, and two levels of grid stretching, including 10 3 and 10 4 these two wall cell aspect ratios (AR), are tested. In each level, five sets of triangular grids from the coarsest to finest are generated within x, y ∈ [0.05, 1.05] × [0, 0.001].
During the grid generation process, nodes in x-direction are equidistantly distributed, while the y-coordinates of different nodes are determined by whereĥ y is the first layer vertical spacing, and β is the stretching factor that could be (a) µ = 10 -6 (b) µ = 10 -8  computed by the known condition y N = 10 −3 . Besides, distribution of five sets of background quadrilateral grids from the coarsest to finest is shown in Table 2.
On this basis, computational accuracy and discretization errors of four different stencils on integral unstructured finite volume solver are counted and given in Figs. 13 and 14.

Computational results on regular grids
From Figs. 13 and 14, we can easily find the second-order accuracy is achieved by all stencils we tested for both L 2 and L ∞ errors. Therefore, the effectiveness of employing face-areaweighted centroid on unstructured finite volume method from integral form is verified. In addition, it is also proved that the derivation of k-exact reconstruction based on any local origins is feasible, and its correctness has been well demonstrated by computational results.
Besides, combining the specific data listed in Table 3, we find both L 2 and L ∞ errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested, and when AR = 10 4 , computational accuracy of L ∞ errors could be obviously improved. What's more, from Table 3, we notice that G-Stencil requires the least number of stencil cells, and the efficiency can also be improved.

Computational results on randomly perturbed grids
On randomly perturbed triangular grids, we also test the case of AR = 10 3 and 10 4 . For simplicity, results of AR = 10 4 are given here.
From Fig. 15, we find errors on randomly perturbed grids demonstrate the similar trends to results on regular grids, and when face-area-weighted centroid is employed on the global-direction stencil, both L 2 and L ∞ errors are greatly reduced. As a result, on linear convective governing equation, the effectiveness as well as superiorities of this novel method is well verified.

Discussion about results exhibition of finite volume discretization from integral form
In this work, we also realize that no matter what numerical scheme is adopted, the main target is obtaining a more accurate result or flowfield, and as mentioned above, for unstructured finite volume method from integral form, the obtained solution vector is the cell-averaged value. Consequently, the result is also exhibited by the cell average, and there is no relevant research to further find the point-valued solution and compare it with the cell average, since the point-valued solution could not be obtained by solving the integral governing equation directly. But we notice that by k-exact reconstruction, the solution at the local origin could be obtained by the mean constraint shown as follows, Besides, when counting the cell-averaged errors, the source term and analytical solution need to be integrated over the control volume, and in this section, three-point quadrature rule is employed [34,35]. For a triangular cell, the integral point is just per face centroid and the corresponding weights are all equal to 1 3 . For the convenience of analysis, the cell average and point value at the geometric centroid and face-area-weighted centroid are all based on the identical global-direction stencil, i.e., G-Stencil. Besides, "Avg" and "Ref" these two labels in Figs. 16 and 17 represent the cell-averaged and point-valued errors respectively.
Combining the results shown in the two figures and Table 4, we notice that no matter the spatial discretization is based on the geometric centroid or face-area-weighted centroid, cell averages are almost identical and both of them are close to point value at the geometric centroid.
In addition, both cell-averaged errors and errors at the geometric centroid are always higher than that at the face-area-weighted centroid. On this basis, we wonder whether the computational result could be exhibited by the value at a certain local origin, that is u ori j , which could be obtained by the k-exact reconstruction. By further consideration, we realize there may be three requirements that should be satisfied if the computational result exhibited by the value at a certain local origin. The first (a) L 2 errors (b) L ∞ errors Therefore, combining the analysis and numerical verification, we notice that compared with the traditional finite volume scheme from integral form, result exhibited at the specific face-area-weighted centroid with the global-direction stencil is more accurate, and the designed order of accuracy could also be achieved.
In short, although the spatial discretization is from integral form, it is not necessary to exhibit the result by the cell average, and the point value at the face-area-weighted centroid is more appropriate to reflect the computational result. Hence, in the following test, errors are all computed for the point-valued solutions recovered at the local origin.

Supersonic vortex flow (governed by Euler equations)
To further test the effectiveness and feasibilities of this novel method on Euler equations, in this section, the supersonic vortex flow is introduced. Computational domain is two concentric circular arcs with radius r i = 1 and r 0 = 1.384 located in the first quadrant. These two circular arcs represent the inviscid wall boundary, and the flow at both inlet and outlet is supersonic. Analytical solution [47] of this numerical example could be derived by isentropic relation and is given as follows, where the value of Mach number at the inner radius is M i = 2.25 and the density ρ i = 1. Besides, the sound speed is calculated as Flow structure of this numerical case is shown in Fig. 18. For this numerical example, both regular and randomly perturbed triangular grids are utilized. As shown in Fig. 19, where the regular grid is generated by splitting the background quadrilateral grid with right diagonals, and randomly perturbed grid is generated by introducing the random node perturbation to the regular grid with topology and the number of cells unchanged. Besides, two levels of grid aspect ratios are employed, and in each level, five sets of grids from the coarsest to finest are generated.
Considering that the aspect ratio of grid with the curved boundary is not a fixed value, similarly, the wall cell aspect ratio is utilized, and two aspect ratios are approximately equal to 4 and 8 respectively. The distribution of background quadrilateral grids in the radical and circumferential is shown in Table 5.

Computational errors on regular grids
As shown in Figs. 20 and 21, for three stencils with the geometric centroid, both L 2 and L ∞ errors of global-direction stencil are the lowest. On this basis, we employ the face-area-weighted centroid, and discretization errors are further reduced. Besides, according to the specific data listed in Table 6, higher-order accuracy on both L 2 and L ∞ errors is achieved by G-Stencil (F-a-w centroid). Particularly, L 2 errors can reach 2.418 and 2.472 order accuracy. As a result, the employment of face-area-weighted centroid on unstructured finite volume method from integral form will not deteriorate the designed order of accuracy but greatly improve the computational accuracy of the finite volume solver after being introduced in global-direction stencil.

Computational errors on randomly perturbed grids
In order to adequately illustrate the numerical performance of different methods, errors on randomly perturbed grids are also counted. For simplicity, results of AR ≈ 8 are given in Fig. 22.  Fig. 22, errors exhibit the similar trends to that of regular grids, and both average and max errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested. Therefore, correctness of k-exact reconstruction based on any local origins as well as the effectiveness and superiorities of the global-direction stencil with face-area-weighted centroid is verified on Euler equations as well.

Subsonic flow over a NACA0012 airfoil (governed by Euler equations)
In Section 3.1, we give a brief introduction about the determination of global directions based on wall distance [42], and in this section, the subsonic flow over a NACA0012 airfoil governed by Euler equations is simulated to examine the effectiveness of the mentioned method. The angle of attack is α = 0 ∘ , and the initial condition is Ma = 0.5.
Regular and randomly perturbed triangular grids with both O-type and C-type topologies are utilized in this numerical example. For O-type regular grid, the first layer spacing in normal direction is 10 −3 , and there are 201 and 71 grid points distributed in circumferential and normal directions respectively. Besides, for C-type regular grid, the background quadrilateral grid is split with random diagonals, and the distribution is 226 × 66 (129 points on airfoil surface). Regular and randomly perturbed grids with both O-type and C-type grids near the airfoil are shown in Figs. 23 and 24.   Besides, to intuitively contrast the effect of geometric centroid and face-areaweighted centroid on cells adjacent to the airfoil, we display these two reference points on randomly perturbed grid with O-type respectively.
From Fig. 25, we can easily find that compared with the geometric centroid, the distribution of face-area-weighted centroid is more regular, and the reference points are almost spread along the normal direction of the wall, especially on grid adjacent to the airfoil surface. On this basis, we count lift and drag coefficients as well as the residual of different stencils with O-type and C-type grids respectively. Before calculating these two coefficients, pressures at the per gauss point of airfoil boundary face are needed at first, and they are interpolated by the point value at the local origin, and the corresponding formula is given as follows, where, p iG is the pressure at the gauss point of wall-adjacent face. Similarly, for the second-order accurate unstructured finite volume discretization, the gauss point is just the face centroid. Besides, (x iG , y iG ) is the coordinate of the current gauss point, and ðx ori i ; y ori i Þ is the coordinate of the local origin C i .   Similar phenomenon could be concluded on drag coefficient, where the result of G-Stencil (F-a-w centroid) is lower than V-Stencil and G-Stencil (Cell centroid) but slightly higher than F-Stencil. Nevertheless, oscillations of F-Stencil are quite obvious, and the convergence property of this stencil is much poorer than another three stencils.

Results on regular and randomly perturbed O-type grid
Conclusions related to lift and drag coefficients mentioned above can be better illustrated by residuals of different stencils. As exhibited in Fig. 28, we can easily find that for three stencils with geometric centroid, the convergence speed of G-Stencil is close to the commonly used V-Stencil, and in the same time steps, residuals of these two stencils are decreased by 7 orders of magnitude, while only 4 orders of magnitude is decreased on F-Stencil. On this basis, with the introduction of face-area-weighted centroid on global-direction stencil, the residual is further decreased, and therefore, the convergence property of integral unstructured finite volume solver is greatly improved by global-direction stencil with face-area-weighed centroid. In the following section, we will further examine the numerical performance on C-type grids.

Results on regular and randomly perturbed C-type grid
From the results shown in Figs. 29 and 30, we can easily find on C-type grid, more accurate lift and drag coefficients are obtained by the global-direction stencil with face-areaweighted centroid as well. However, oscillations of the face-neighbor stencil are extremly evident compared with another three stencils we tested, and the result does not converge at all. It could be well illustrated by the residual shown in Fig. 31.
As shown in Fig. 31, vertex-neighbor and global-direction stencils have a similar convergence rate, and it is improved by the global-direction stencil with face-area-weighted

Dissipative term evaluation (governed by Laplacian equation)
From Section 4.1 to 4.3, three numerical examples govern by linear convective and Euler equations are simulated, and the correctness of k-exact reconstruction based on any local origins is verified. Moreover, on unstructured finite volume method from integral form, the global-direction stencil with face-area-weighted centroid also has a better computational accuracy, efficiency as well as convergence rate for convective flux evaluation.
But, the numerical performance on dissipative term also needs to be further tested to set stage for the extension to viscous problems. In this section, the MMS method is also used on Laplacian equation to examine discretization errors and convergence speed of different stencils. Laplacian model equation could be formulated as where ν is a constant equaling to 1. Integrating this model equation over the control volume, we obtain (a) Geometric centroid (b) Face-area-weighted centroid According to the divergence theorem, Eq. (34) could be transformed as From Eq. (35), we can clearly find the difference between convective term and dissipative term. For dissipative term computation, what we need to compute at the face integral point is not two state vectors obtained by owner and neighbor cells, but solution gradient, which could be evaluated by the arithmetic average of left and right cells for the second-order unstructured finite volume method, On this basis, in order to improve the stability and reduce the truncation errors, the average face gradient is always added a solution jump term [28,48,49]. Here are many different schemes for adding the jump term, such as the edge-normal (EN) and face-tangent (FT) schemes [50]. Based on the FT scheme, a novel α-damping scheme is proposed by Nishikawa [51][52][53][54]. The accuracy as well as robustness of second-order unstructured finite volume solver is greatly improved by this scheme, and it could be formulated as (a) C d on regular grid (b) C d on perturbed grid  where, α is a damping coefficient, and here, a special value is α = 4/3, which has been known to provide superior accuracy and robustness [28,48,51].n jk is the unit face normal vector, andê jk is the unit vector between two cell centroids. L jk is the length from one centroid to another.
Besides, ϕ Ã L and ϕ Ã R are solutions linearly reconstructed at integral point, On the basis of determining the discretization method of dissipative term, the manufactured solution of Eq. (33) is where ω = 2π is a constant, and the source term could be derived as As shown in Fig. 32, the simulation is carried out on five sets of triangular grids generated by splitting the quadrilateral grids with regular or random diagonals, and the cell aspect ratio is AR = 10. The distribution of these five sets of grids is from 10 × 10 to 80 × 80, and the computational domain is x, y ∈ [−0.5, 0.5] × [−5 × 10 −2 , 5 × 10 −2 ].
(a) C l on regular grid (b) C l on perturbed grid  From Eq. (37), we notice that the damping term is greatly influenced by the grid skewness measure. On this basis, we first count it on both two grid types, and list the results in Table 7.
As a result, superiorities of face-area-weighted centroid are illustrated again, and with the employment of this novel centroid, the grid skewness is almost eliminated. In order to reflect the relative position of two different centroids on grid with regular and random diagonals, the distribution of these two local origins is shown in Fig. 33.
From Fig. 33, we can easily find that face-area-weighted centroids are almost distributed along y-axis no matter on grid with regular or random diagonals, but the distribution of geometric centroid is irregular, and according to the data listed in Table 7, it has been demonstrated that the evident skewness will be introduced by this traditional centroid. Based on above analysis, we count discretization errors of four different stencils, and corresponding results are shown in Figs. 34 and 35 and Table 8.
From Figs. 34 and 35 and combining the specific data listed in Table 8, we find that discretization errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested no matter on grid with regular or random diagonals, and the computational accuracy is greatly improved by this novel method.
On this basis, residuals on the coarsest and finest regular grids of different stencils are shown in Fig. 36.
From Fig. 36, we can easily find that among three stencils with the geometric centroid, the G-Stencil (Cell centroid) has a faster convergence speed than two commonly used V-Stencil and F-Stencil. What's more, when face-area-weighted centroid is employed, the convergence rate is further promoted. Hence, except for convective term, on dissipative term, a better computational accuracy and faster convergence rate  are also realized by this novel method. In short, feasibilities and superiorities of employing face-area-weighted centroid on unstructured finite volume discretization from integral form are verified again on Laplacian model equation.

Accuracy analysis of the method on the third-order discretization
In Section 2.2, the k-exact reconstruction based on any local origins is derived, and from Section 4.1 to 4.4, we notice that on the second-order finite volume solver, the result of global-direction stencil with face-area-weighted centroid is much more accurate than two commonly used vertex-neighbor and face-neighbor stencils, as well as the global-direction stencil with the geometric centroid. In this section, the k-exact reconstruction and spatial discretization based on any local origins are derived at first. On this basis, a typical numerical example based on the Method of Manufactured Solutions (MMS) [43][44][45][46] is employed for the accuracy analysis.

Derivation of higher-order k-exact reconstruction based on any local origins
For higher-order k-exact reconstruction, e.g., the third-order, the reconstructed polynomial based on any local origins could be formulated as where ðx ori j ; y ori j Þ is the coordinate of the current local origin, which is not a necessary geometric centroid, and if the ðx − x ori j Þ and ðy − y ori j Þ are replaced by ðx − x gc j Þ þ ðx For k-exact reconstruction, no matter the reconstructed function is formulated based on which reference point, the mean constraint requires [34,35] Besides, integrating the reconstructed polynomial over the control volume, we have Computational errors on the finest regular grid and L 2 , L ∞ accuracy between the last two grids of different stencils This mean constraint could be simplified as which is utilized for calculating the point value u ori j , and after that, the mean constraint is always satisfied. If we take T j to represent the sum of some terms we have To achieve kth-order accuracy requires that we compute the kth derivatives by minimizing the error in predicting the mean value of the reconstructed function for the stencil k, that is, by minimizing the difference between the actual cell average u k and the average of R j ðx − x ori j Þ over stencil cell k. The mean value for a single cell V k of the reconstructed function R j is Likewise, ðx − x ori j Þ and ðy − y ori j Þ could be written as ðx − x gc k Þ þ ðx gc k − x ori j Þ and ðy − y gc k Þ þ ðy gc k − y ori j Þ respectively, Combining the mean constraint given in Eq. (48) and writing Eq. (50) for every stencil cell, we obtain 1 x j y j x 2 j xy j y 2 where, the first row is the mean constraint, and the unconstraint equations are solved at first. After that, the point value at any local origins should be computed by the mean constraint to ensure Eq. (48) is exactly satisfied, which is utilized to calculate the solution at per gauss point of cell face. Eventually, the k-exact reconstruction based on any local origins is derived, and in the next section, higher-order spatial discretization will be specifically described.

Higher-order finite volume spatial discretization
In Section 2.1, the second-order finite volume discretization has been introduced, and as shown in Fig. 1, we mentioned that only one gauss point, i.e., the face centroid, is required for the second-order accurate discretization, while there should be at least two gauss points to guarantee the flux integration is higher-order accurate. The required order of accuracy, gauss point location as well as corresponding weights are given in Table 9, where A and B represent two edge points of a certain cell face. As a result, for the third-order accurate discretization, two gauss points are employed, and the left and right state vectors at per gauss point could be interpolated by the point value u ori j and u ori k , and both of them are obtained by the mean constraint of k-exact reconstruction.
Accordingly, the numerical flux could also be obtained by the Roe flux [37].

Accuracy analysis
In the first two sections, the third-order k-exact reconstruction as well as the spatial discretization is presented in detail. In this section, a numerical example, which is governed by the Euler equations, based on the Method of Manufactured Solutions (MMS) [43][44][45][46] is employed for the accuracy analysis. The manufactured solution [28] is shown as follows, Taking the analytic solution to Euler equations, the source term is obtained where H = (γp/(γ − 1) + ρv 2 /2)/ρ is the specific total enthalpy. For higher-order unstructured finite volume discretization, the source term should be integrated over the control volume rather than being evaluated at the reference point to avoid the second-order error. Three-point quadrature [34,35] is utilized again for the source term integration. Besides, regular and randomly perturbed triangular grids with 5 × 10 2 and 10 3 these two aspect ratios are employed, and in each cell aspect ratio, five sets of triangular grids from the coarsest to finest are utilized to count the discretization errors and computational accuracy. The coarsest regular and randomly perturbed triangular grids and the distribution of background quadrilateral grids in x and y directions are shown in Fig. 37 and Table 10.
Computational results on both regular and randomly perturbed triangular grids are shown as follows.

Computational results on regular triangular grids
As shown in Figs. 38 and 39, the third-order accuracy is achieved by employing the globaldirection stencil with face-area-weighted centroid. Therefore, correctness as well as feasibilities of k-exact reconstruction and spatial discretization based on any local origins is demonstrated. Besides, combining the result given in Table 11, we can easily find that discretization errors of global-direction stencil with face-area-weighted centroid are always the lowest among all methods we tested, no matter on grids with AR = 5 × 10 2 or 10 3 .
But note that errors are also counted at the local origin rather than the cell-averaged result, and for higher-order accurate finite volume method, the comparison between the cell-averaged result and that counted at the local origin is listed in Section 5.3.3.

Computational results on randomly perturbed triangular grids
From Figs. 40 and 41, we find that on randomly perturbed triangular grids, there is a similar trend to results on regular grids, that is, the global-direction stencil with facearea-weighted centroid is lower than another three stencils. As a result, computational accuracy is improved, which is beneficial for obtaining a more accurate flowfield. In short, feasibilities and superiorities of extending the face-area-weighted centroid and the combination of global-direction stencil with this novel local origin to higher-order unstructured finite volume method from integral form are verified.

Discussion about results exhibition of higher-order discretization from integral form
In Section 4.1.3, it is verified by analysis and numerical example that for the second-order unstructured finite volume discretization, the point value at the face-area-weighted centroid is more appropriate for the result exhibition, since errors at this certain local origin (a) Regular grid (a) Randomly perturbed grid According to the three requirements mentioned in Section 4.1.3, firstly, the point value at any local origins u ori j could be obtained by the mean constraint. In addition, as shown in Sections 5.3.1 and 5.3.2, the third-order accuracy is always achieved by the point-value solution at both geometric centroid and the face-area-weighted centroid, and more importantly, as shown in Fig. 42, errors counted at the face-area-weighted centroid are obviously lower than that counted at the geometric centroid as well as the cell-averaged value.
In summary, for higher-order accurate unstructured finite volume method from integral form, when the global-direction stencil is utilized, we find the result counted at the facearea-weighted centroid is more accurate than the cell-averaged solution. Hence, feasibilities as well as the superiorities of the face-area-weighted centroid are demonstrated not only on the second-order finite volume method, but on discretization with higher-order accuracy.

Concluding remarks
Inspired by the novel face-area-weighted centroid utilized for unstructured finite volume discretization from differential form, in this paper, we have shown that the reference point used in the integral form does not have to be the geometric centroid as (a) L 2 errors (b) L ∞ errors well, and the face-area-weighted centroid is a better choice, which effectively reduces the grid skewness. Besides, we combine the global-direction stencil with this novel reference point trying to improve the computational accuracy and convergence rate of unstructured finite volume solver. Specifically, we first illustrate in detail that the traditional LSQR reconstruction method is unable to obtain the point-value solution located at any local origins, and that is the reason why we focus on k-exact reconstruction algorithm. On this basis, the k-exact reconstruction based on any local origins is analytically derived, and during the reconstruction process, the mean constraint can always be satisfied. Besides, we extend the global-direction stencil with novel face-area-weighted centroid from unstructured finite volume method in differential form to the integral form for the grid skewness reduction and a better reflection of flow anisotropy. More importantly, it sets stage for the promotion on higher-order accuracy.
Four representative numerical examples governed by linear convective, Euler and Laplacian equations are simulated to examine the correctness of k-exact derivation and the effectiveness of global-direction stencil with face-area-weighted centroid on the unstructured finite volume discretization from integral form. After verification, in numerical cases governed by linear convective and Euler equations, when the face-area-  weighted centroid is utilized on finite volume discretization, the second-order accuracy also could be achieved, and the global-direction stencil with this novel reference point has the lowest discretization errors among all stencils we tested, where, the result is exhibited by the point value rather than the cell average, since we realize when the global-direction stencil is utilized, the point-value solution u ori j obtained by the mean constraint of k-exact reconstruction at the face-area-weighted centroid is more accurate than the cell average, and this conclusion is also demonstrated on the third-order finite volume discretization. On this basis, the result exhibition of unstructured finite volume method could be replaced by this point value at the face-area-weighted centroid. Besides, in subsonic flow over NACA0012 airfoil, pressure values utilized to compute the lift and drag coefficients at the gauss point of wall-adjacent face are also interpolated by the point-valued solution at the local origin, and we find this novel method has a better numerical performance on calculating lift and drag coefficients, and has a faster convergence speed than commonly used stencils based on the geometric centroid. On this basis, a Laplacian model equation is utilized to further test its effectiveness on dissipative term evaluation. From the result, similar conclusions are demonstrated that errors of the global-direction stencil with face-area-weighted centroid are still lower than another three stencils, and the convergence rate is also greatly promoted by this novel method. As a result, this conclusion sets stage for the further extension to viscous flows.
Finally, the accuracy analysis is implemented on the third-order unstructured finite volume discretization from integral form as well, where the third-order k-exact reconstruction as well as the spatial discretization based on any local origins is derived. After verification, it is illustrated that the third-order accuracy is also achieved by the globaldirection stencil with face-area-weighted centroid, and discretization errors of this method are lower than another three stencils we tested.
In conclusion, the computational accuracy of the second and third-order unstructured finite volume method from integral form will not get deteriorated by the employment of face-area-weighted centroid, and it is feasible to extend the global-direction stencil with face-area-weighted centroid from differential finite volume solver to the integral form. Future work will be carried out from two aspects. Firstly, we will not only further compare the computational accuracy of point-valued solution at the certain local origin with the abstract cell-averaged solution, but figure out the specific approach to display the flowfield by this point value. In addition, extensions of the global-direction stencil selection method and face-area-weighted formula to three dimensions are also quite necessary for practical applications.