On the application of polygonal finite element method for Stokes flow – A comparison between equal order and different order approximation

Journal Pre-proof On the application of polygonal finite element method for Stokes flow – A comparison between equal order and different order approximation

Sundararajan Natarajan

PII:

S0167-8396(19)30122-0

DOI:

https://doi.org/10.1016/j.cagd.2019.101813

Reference:

COMAID 101813

To appear in:

Computer Aided Geometric Design

Please cite this article as: Natarajan, S. On the application of polygonal finite element method for Stokes flow – A comparison between equal order and different order approximation. Comput. Aided Geom. Des. (2020), 101813, doi: https://doi.org/10.1016/j.cagd.2019.101813.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.

Highlights • • • •

Proposes an extension polygonal finite element for Stokes flow with equal and different order approximation for velocity and pressure. For equal order approximation, a pressure stabilized Petrov-Galerkin method is employed to ensure stability. Strain smoothing technique is employed to compute the terms in the bilinear linear form. Yields optimal convergence rates in L 2 norm and H 1 seminorm.

On the application of polygonal finite element method for Stokes flow - a comparison between equal order and different order approximation Sundararajan Natarajana,1,∗ a

Integrated Modelling and Simulation Lab, Department of Mechanical Engineering, Indian Institute of Technology-Madras, Chennai - 600036, India.

Abstract In this work, we discuss the application of polygonal finite elements to the solution of Stokes equations in two dimensions. The proposed framework is based on the mixed finite element discretization that involves velocity and pressure as independent unknown field variables. We present two approaches: (a) same order of approximation for velocity and pressure with Pressure Stabilization Petrov-Galerkin method (PSPG) and (b) different order of approximation for velocity and pressure. The approximation functions over polygons are rational polynomials based on Wachspress coordinates and for numerical integration of the terms in the bilinear and the linear form, we employ the linear smoothing technique. The relative performance between the approaches, the convergence properties and the accuracy is presented for two dimensional numerical examples, which shows that the proposed framework yields accurate results and converges at optimal convergence rate. Keywords: Polygonal finite element method, Wachspress interpolants, Numerical integration, PSPG, Stokes equation, Serendipity shape functions. 1. Introduction Introduction of elements with arbitrary edges/facets into the finite element method (FEM), coined as ‘polygonal FEM (PFEM)’ has revolutionized the FEM. This constitutes a class of approximation method that retain most of the favorable properties of the conventional FEM but without finite-element-like restrictions on the element geometry and topology. Such schemes offer the promise of greatly simplified automatic mesh generation, even on extremely complex and/or evolving domains. Moreover, an adaptive mesh refinement and coarsening are simplified with polytopes since they naturally ∗

Corresponding author Email address: [email protected] (Sundararajan Natarajan) 1 Department of Mechanical Engineering, Indian Institute of Technology-Madras, Chennai - 600036, India.

Preprint submitted to Journal Name

January 6, 2020

address the difficulty in handling hanging nodes [1]. However, the success of the PFEM hinges on two areas: (a) a robust mesh-generation approach that can handle complex geometries and (b) construction of approximation functions over arbitrary polytopes. The polygonal meshes can be generated from: (a) Voronoi tessellations [2–4] or (b) by connecting the centroid of all triangular elements circumventing a particular node in a triangular mesh [5, 6]. Interestingly, the approximation functions over arbitrary polytopes are not unique, some of approaches include: Wachspress basis functions [7], Laplace interpolants [8], mean value coordinates [9], harmonic coordinates [10], blended coordinates [11], to name a few. These different approaches are collectively referred to as ‘Generalized barycentric coordinates’. Interested readers are referred to the recent book by Hormann and Sukumar [12] for more detailed discussion on the construction and applications of the PFEM. Due to the aforementioned advantages, the development of discretization techniques over polytopes has gained considerable attention in the recent past. Some of the popular techniques include: polygonal finite element method [13], scaled boundary FEM [14], the virtual element method [15], strain smoothing technique [16, 17], virtual node method [18], to name a few. Since its inception, the polygonal/polyhedral elements have been used to solve problems involving large deformations [19], incompressibility [20], contact problems [21] and fracture mechanics [22]. Although, a lot of research has been focused on employing the VEM to steady and unsteady Navier-Stokes equations [15, 23–27], to the best of author’s knowledge, the application of polygonal finite element method with Generalized Barycentric coordinates for Stokes flow is scarce in the literature [20, 28]. Talischi et al. [20] presented lower order spaces over convex polygons for incompressible flow problems. Chi et al., [29, 30] in a series of papers studied the performance of polygonal elements for finite elasticity including both lower and higher-order elements. Vu-Huu et al. [28] developed MINI polygonal element, wherein the approximation space of the velocity field is supplemented with a bubble node. The pressure degrees of freedom (dofs) are at the vertices and the velocity dofs are at the vertices and at the center of the element. The main objectives of the paper are to study the convergence properties and the accuracy of similar order and different order approximation for the velocity and the pressure. In case of similar order approximation, a pressure stabilized Petrov-Galerkin method is employed to ensure stability, whilst in the latter case, the stability is ensured by the different order of approximations employed for the velocity and the pressure. The recently introduced quadratic serendipity elements are used for approximating the unknown velocity fields and the linear smoothing technique is employed to integrate the terms in the bilinear and the linear form. The accuracy is demonstrated with two examples for which analytical solutions are available. The paper is organized as follows. The governing equations for Stokes flow and the corresponding weak form is presented in Section 2. Section 3 presents the spatial discretization. An overview of the linear smoothing technique to integrate the terms in the bilinear linear form is presented in Section 4. Two numerical examples are presented in Section 5 to study the convergence properties and accuracy of the proposed 2

approaches, followed by conclusions in the last section. 2. Governing equations Consider a domain Ω ∈ R2 with the boundary ∂Ω. In this study, we consider the steady state Stokes equation with homogeneous Dirichlet boundary conditions: find (u, p) such that −νΔu − ∇p =f in Ω div u =0 in Ω u =0 on ∂Ω

(1)

where Δ, div represents the vector Laplacian, the divergence, respectively, whilst, ∇ and ∇ denote the gradient operator for vector and scalar fields, respectively. ν ∈ R, ν > 0 is the viscosity, u, p are the velocity and the pressure fields, respectively and f is the external force. By following the standard Galerkin procedure and considering the spaces: ⎧ ⎫  ⎨ ⎬ V := [H01 (Ω)]2 , Q := L20 (Ω) = q ∈ L2 (Ω)| q dΩ = 0 , (2) ⎩ ⎭ Ω

the corresponding variational formulation for Equation (1) is: find (u, p) ∈ V × Q, such that: a(u, v) + b(v, p) = (f , v) b(u, q) = 0 where, ∀u, v ∈ V and ∀q ∈ Q:  a(u, v) := ν∇u : ∇v dΩ

(3)

 b(v, q) :=



Ω

q div u dΩ, Ω

f · v dΩ

(f , v) :=

∀v ∈ V ∀q ∈ Q

(4)

Ω

It is observed that the finite element method when applied to Equation (3) suffers from instabilities due to incompatibility between the velocity and pressure finite element spaces. For the two-field mixed finite elements, approximations of the additional pressure field is needed. As conforming approximations of Q, either discontinuous or continuous approximations can be adopted. Within the continuous approximation for the pressure field framework, there are two approaches: equal and different order approximation for the velocity and the pressure. When using equal order of approximation for both the 3

velocity and the pressure, it has been observed that the node-to-node pressure shows oscillations. This is because, the equal order approximation fails to satisfy the LBB or inf-sup condition, which requires that the approximation order of the pressure field should be one order lower than the approximation order of the velocity field. The other approach would be to use one order higher approximation for the velocity field when compared to the pressure field, for example, a quadratic variation for the velocity and a linear variation for the pressure. Although this approach yields stable results, the velocity and the pressure fields have different dofs, which makes tracking of these field variables difficult. On the contrary, research has been focused towards using equal order approximation with stabilization techniques, for example, PSPG, Galerkin Least squares, Streamline upwind Petrov-Galerkin (SUPG). A concise discussion on different stabilization approach is presented in [31]. In this paper, we consider both types of polygonal finite elements: equal and different order approximation for the velocity and the pressure field. For the equal order of approximation, we employ the PSPG method for stabilization. Based on this, the modified weak form is:    (5) ν∇u : ∇v dΩ + q div u dΩ + τ ∇q (−νΔu − ∇p − f ) dΩ = 0 Ω

Ω

Ω

where τ ∇q is the perturbation multiplied to the residual of the momentum equation and τ is the stabilization parameter. In this study, we set τ=

h2e 4ν

where he is the maximum size of the circle that can be inscribed inside the polygonal element. For the current study, −νΔu is set to zero. Figure 1 shows a schematic representation of the degrees of freedom in a typical polygonal element for the velocity and the pressure. In case of an equal order approximation, both the velocity and the pressure dofs are at the vertices and in case of different order approximation, the pressure dofs are at the vertices and the velocity dofs are at the vertices and at the mid-point of the edges. 3. Spatial discretization and construction of basis functions In the proposed approach, the domain Ω is discretized into a sequence of decompositions Ωh , where the decomposition consists of arbitrary polygonal elements E with hE := diameter(E) and h := sup hE . These elements satisfy: (a) star-convexity and E∈Ωh

(b) the distance between any two vertices of E ≥ ρhE , where ρ is a positive constant. The construction of basis functions over arbitrary polygons are not unique and are collectively termed as ‘Generalized Barycentric coordinates’ [12]. In this study, we employ Wachspress coordinates for linear polygons and its serendipity extension for quadratic polygons. For a convex polygon with n vertices (see Figure 2), Wachspress [7], by using the principles of perspective geometry, proposed rational basis functions on polygonal 4

Pressure dof Velocity dof

(a)

(b)

Figure 1: Schematic representation of an arbitrary polygonal element with pressure and displacement Dofs: (a) equal order approximation and (b) different order approximation.

elements, in which the algebraic equations of the edges are used to ensure nodal interpolation and linearity on the boundaries. The generalization of Wachspress shape functions to simplex convex polyhedra was given by Warren [32, 33]. The construction of the coordinates is as follows: Let P ⊂ IR2 be a simple convex polygon with edges F and vertices V . For each edge f ∈ F , let nf be the unit outward normal and for any x ∈ P , let hf (x) denote the perpendicular distance of x to f , which is given by hf (x) = (ve − x) · nf

(6)

for any vertex ve ∈ V that belongs to f . For each vertex ve ∈ V , let f1 , f2 be the two edges incident to ve and for x ∈ P , let wve (x) =

det(pf1 , pf2 ) . hf1 (x)hf2 (x)

(7)

where, pf := nf /hf (x) is the scaled normal vector, f1 , f2 , · · · , fd are the d edges adjacent to ve listed in a counter-clockwise ordering around ve as seen from outside P (see Figure 2) and det denotes the regular vector determinant in Rd . The shape functions for x ∈ P is then given by wve (x) λve (x) = . (8) wue (x) ue∈V

The Wachspress shape functions are the lowest order shape functions that satisfy boundedness, linearity and linear consistency on convex polytopes [32, 33]. A simple MATLAB implementation is given in [34]. In this work, the quadratic serendipity shape 5

Figure 2: Barycentric coordinates: Wachspress interpolants, where x is any point in P , a convex polygon.

3

3

3

4

3

7

7 4

6

7

4

6

4

6

9 2

8

2

8

2

5 1

λij

Linear

1

1

μij

Quadratic

8

2

5 ξij

Serendipity

5 1 ij

Lagrange type

Figure 3: Construction of quadratic serendipity shape functions based on generalized barycentric coordinates.

6

functions are constructed from pairwise product of set of barycentric coordinates [35]. These shape functions are C ∞ continuous inside the domain Ω while C 0 continuous at the boundary Γ . The pairwise product of a set of barycentric coordinates (λi ), yields a total of n(n + 1)/2 functions with mid-nodes on the edges and nodes inside the domain. Then by appropriate linear transformation technique, the n(n + 1)/2 set of functions are reduced to 2n set of functions that satisfies Lagrange property. The essential steps involved in the construction of quadratic serendipity shape functions are (also see Figure 3): 1. Select a set of barycentric coordinates λi , i = 1, · · · , n, where n is the number of vertices of the polygon. 2. Compute pairwise functions μab := λa λb . This construction yields a total of n(n + 1)/2 functions. 3. Apply a linear transformation A to μab . The linear transformation A reduces the set μab to 2n set of functions ξij indexed over vertices and edge midpoints of the polygon. 4. Apply another linear transformation B that converts ξij into a basis ψij which satisfies the “Lagrange property.” Interested readers are referred to [35] for a detailed discussion on the construction of quadratic serendipity elements. Recently, the construction of serendipity elements has been extended to arbitrary polyhedra in three dimensions [36]. The approximation functions discussed earlier are used to define the discrete trial and test spaces by constructing the shape functions over the union of all nel ∈ Ω. The trial functions and test are written as a linear combination of the shape functions: uh = φ I uI ; p h = φI pI , a

I

where φ are the approximations functions, depending on the choice, they could be • Wachspress interpolants (φ ≡ λ) for equal order approximation; or • the serendipity shape functions (φ = ψ) for node I of the element for different order approximation, with Wachspress interpolants for pressure and serendipity shape functions for velocity. We close this section by discussing the numerical integration employed in this paper to compute the terms in the bilinear and linear forms. The shape functions over arbitrary polygons are typically rational functions (polynomials over polynomials) and hence conventional quadrature rules cannot be employed [37]. Due to the nonpolynomial nature of the shape functions, typically higher order quadrature rules are employed to improve the accuracy, however, this leads to large amount of integration points. The numerical integration over arbitrary polygons has attracted considerable attention in the recent past. They include: Green-Gauss quadrature [38], complex mapping [39], generalized quadrature [40], homogeneous numerical integration [41] and 7

strain smoothing [16, 17, 42, 43], to name a few. In this paper, we employ the recently introduced linear smoothing technique [17], an overview of which is presented in the next section. 4. Overview of the linear smoothing technique In this section, a numerical integration scheme known as linear Smoothing (LS)[44] for polygons is discussed. With LS technique the accuracy of the polygonal element is obtained at lower computational expense in comparison to the conventional PFEM, which typically requires ‘more’ integration points to compute the terms in the bilinear/linear form (see Equation (3)). The procedure for numerical integration with LS technique is by sub-dividing the polygonal element into smoothing domain (i.e triangular smoothing domain, see Figure 4). Note that the process of sub-division is only for the purpose of numerical integration and does not introduce any additional degrees of freedom. 5 4



3

0; 1)

y

6

c

x1

c

c x3

x2 x

1

c

0; 0)

2

1; 0)



Figure 4: Discretization of arbitrary polygon into triangular smoothing domain ΩC using virtual point shown by ‘open’ circle. The smoothed strain is obtained at open square.

h is evaluated from the standard gradient Now, the smoothed gradient matrix ∇u h operator ∇u (x) as given below:  h  ∇uh (x) q(x) dV (9) ∇u = h ΩC

where q(x) is the smoothing function. On writing Equation (9) at the basis functions derivative level and invoking Gauss-Ostrogradsky theorem, we get:    φI,j q(x) dV = φI q(x)nj dS − φI q,x (x) dV (10) h ΩC

h ΓC

h ΩC

8

where, j = 1, 2 represents the x and y coordinates, the linear smoothing function q(x) = {1, x, y} in two dimensions. Note that the domain integral in Equation (10) is evaluated at the points shown by open square within the smoothing domain, xm = (xm , ym ) (see Figure 4), where m = {1, 2, 3}. Whilst, the boundary integral is evaluated along the boundary of the smoothing domain (the location of integration point on the boundary is represented by ‘filled’ square in Figure 4). On expanding Equation (10), we get  h ΩC

 



h ΩC

h ΩC

φI,1 dV = 

φI,1 x1 dV =  φI,1 x2 dV =

h ΓC

h ΓC

h ΓC

φI n1 dS,

(11a) 

φI x1 n1 dS −

h ΩC

φI dV,

(11b)

φI x2 n1 dS,

(11c)

φI n2 dS,

(11d)

φI x1 n2 dS,

(11e)

for φI,1 , and   

 h ΩC

h ΩC

h ΩC

φI,2 dV = 

φI,2 x1 dV =  φI,2 x2 dV =

h ΓC

h ΓC

h ΓC

 φI x2 n2 dS −

h ΩC

φI dV

(11f)

for φI,2 . The integration is performed within smoothing domain ΩC , where φ is the basis functions for node I of a polygonal element and n = (n1 , n2 ) is the outward normal along the boundary ΓC of the smoothing domain. Let the coordinates of the m-th interior subcell Gauß point be defined as m r = (m r1 , m r2 ) and its associated Gauß weight as m w; the coordinates and the Gauß weight of the g th Gauß point that is located on the k th edge of the subcell is gk s = (gk s1 , gk s2 ) and gk v, respectively; and the unit outward normal to the k th edge of the subcell is denoted as k n = (k n1 , k n2 ). In two dimensions, three interior Gauß points (ngp = 3) per subcell is required to compute the smoothed strain displacement matrix. On using numerical integration in Equation (11) leads to the following system of linear equations: Wdj = f j , j = 1, 2 where

⎡

1

w W = ⎣ 1 w 1 x1 1 1 w x2

2

w 2 2 w x1 2 2 w x2

9

⎤ w 3 3 w x1 ⎦ , 3 3 w x2

(12a)

3

(12b)

⎡

3

⎤ φI (gk s) k n1 gk v

⎢ ⎢ k=1 g=1 ⎢ 3 ⎢ 3 φ (g s) g s n g v − φI (m r) m w f1 = ⎢ I k 1 k 1 k k ⎢ k=1 g=1 m=1 ⎢ 3 ⎣ g g g φI (k s) k s2 k n1 k v ⎡

k=1 g=1 2 3

⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎦

(12c)

⎤ φI (gk s) k n2 gk v

⎢ ⎢ k=1 g=1 ⎢ 2 ⎢ 3 φI (gk s) gk s1 k n2 gk v f2 = ⎢ ⎢ k=1 g=1 ⎢ 3 2 3 ⎣ φI (gk s) gk s2 k n2 gk v − φI (m r) m w k=1 g=1

⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎦

(12d)

m=1

and the solution vector of the j th basis function derivative evaluated at the three interior subcell Gauß points is T    . (12e) dj = 1 dj 2 dj 3 dj = φI,j (1 r) φI,j (2 r) φI,j (3 r) The smoothed derivatives given in Equation (12e) are then substituted in the weak form (see Equations (3) and (5)) and the nodal velocity and pressure are computed using standard FE procedures. 5. Numerical examples The accuracy and the convergence properties of the proposed scheme is studied with two examples in this section. The computational domain considered is a square Ω = [0, 1]2 . The domain is discretized with structured quadrilateral elements and with arbitrary polygonal elements. The two dimensional polygonal meshes are generated using the MATLAB functions in Polytop [45]. The arbitrary polygons are based on centroid Voronoi tessellation. For the purpose of error estimation and convergence studies, the L2 −norm of the velocity and the pressure are used that is plotted against the size of the minimum circle that can be embedded in a polygon. The following convention is employed whilst discussing the results: • Q8/Q4 - second order approximation for velocity and first order approximation for pressure with quadrilateral elements. • Q4/Q4+PSPG - bilinear approximation over quadrilateral elements for both the fields, ie., velocity and pressure with pressure stabilization. • k2/k1 - second and first order approximations over arbitrary polygons for velocity and pressure, respectively. • k1/k1+PSPG - first order approximation for both velocity and pressure over arbitrary polygons with pressure stabilization. 10

Figure 5 shows an example of adopted meshes. In all the test cases, we set ν = 1, unless mentioned otherwise.

(a)

(b)

Figure 5: Representative finite element mesh: (a) structured four noded quadrilateral elements and (b) arbitrary polygonal elements.

Example 1. In this example, the body load f and the Dirichlet boundary conditions are chosen such that the exact solution takes the following form [46]:   sin(2πx) cos(2πy) u(x) = − cos(2πx) sin(2πy) p(x) = x2 + y 2

(13)

Example 2. In this example, the external load f is chosen such that the analytical solution takes the following form [47]:   sin(x) cos(y) u(x) = − sin(y) sin(2πx) cos(x) p(x) = 3x2 + sin(xy) − 1.239811742000564725943866

(14)

In both the examples, the exact solution for the velocity given by Equation (14) is prescribed as Dirichlet boundary condition on all four sides of the domain Ω and the pressure at the bottom left corner is set to zero. Figures 6 -8 illustrate the convergence of the relative error in the L2 −norm of the velocity and the pressure for different discretization types for example 1 and example 2, respectively. From the convergence plots, it can be verified that the same order with PSPG and with different order of 11

10−1 10−2

1

1

10−3

10

−4

1

2

10−2

10−1 h Q8/Q4;

k2/k1;

100

Relative error in the L2 −norm-velocity

Relative error in the L2 −norm-pressure

100

100 10−1 10−2 10

−3

1

1.9

10−4 10

Q4/Q4+PSPG;

−5

1

10−2

3 10−1 h

100

k1/k1+PSPG

100 10−1

10

−2

10

−3

10−4 −2 10

1

1

1

2 10−1 h

Q8/Q4;

k2/k1;

100

Relative error in the L2 −norm-velocity

Relative error in the L2 −norm-pressure

Figure 6: Convergence of the relative error in the L2 − norm of velocity and pressure for example 1.

Q4/Q4+PSPG;

10−1 10−2 10−3 10

−4

1

1.7

10−5 10−6 10−7 −2 10

1

2.7 10−1 h

k1/k1+PSPG

Figure 7: Convergence of the relative error in the L2 − norm of velocity and pressure for example 2.

12

100

Error in the L2 −norm-velocity gradient

100

Q8/Q4 k2/k1 Q4/Q4+PSPG k1/k1+PSPG

1

10−1

1

10−2

2 1

10−3 10−4 −2 10

10−1 h

100

Figure 8: Convergence of the relative error in the L2 − norm of gradient of the smoothed velocity for example 2.

approximation converges at optimal convergence rate in both the velocity and pressure. To further study the robustness and the accuracy of the framework, the numerical infsup test is done for Example 2. The numerical inf-sup test guarantees stability of the mixed formulation that requires the following condition to be satisfied [48–50]  |ph divuh | dΩ Ω inf sup = γh ≥ 0 (15) ph uh ||ph ||L2 ||uh ||H 1 where γh is a positive constant independent of the mesh size h. In this study, we employ the numerical test proposed by Chapelle and Bathe [49] to demonstrate the stability of the mixed formulation. For the numerical test, matrices Gh and Sh are defined by [49]:  ph divuh dΩ = WhT Gh Vh Ω

||ph ||L2 = WhT Gh Vh uh ||H 1 = VhT Sh Vh

(16)

where Wh and Vh are the vectors corresponding to the nodal quantities. Then, the smallest non-zero eigenvalue λ of the following eigenvalue problem is calculated to √ obtain the numerical inf-sup value γh = λ: Gh φ = λSh φ

(17)

Figure 9 shows the numerical inf-sup value γh , evaluated for different mesh discretizations, h. It can be inferred from Figure 9 that the numerical inf-sup value converges to a 13

value that is bounded away from zero with mesh refinement and is almost independent of the mesh size h. The pressure distribution is shown in Figure 10 with and without stabilization when the domain is discretized with arbitrary polygonal elements. The pressure contours are shown for a mesh consisting of 320 polygons. It can be seen that without stabilization, the pressure distribution has unphysical oscillations, whilst, this is alleviated when stabilization is employed.

γh

100

Q8/Q4 k2/k1

10−1 −2 10

10−1 Maximum edge size h

100

Figure 9: Example 2: stability parameter as a function of mesh size for different order approximation for velocity and pressure. The domain is discretized with quadrilateral and arbitrary polygonal elements.

6. Conclusions In this paper, we studied the performance of arbitrary polygons in solving the Stokes equation. Both equal and different order of approximation for the velocity and pressure were considered. In case of same order of approximation, a pressure stabilized Petrov Galerkin approach was adopted to alleviate the nodal pressure oscillations. The convergence properties and the accuracy were studied for two example problems for which analytical solutions are available. From the numerical study, it can be inferred that: (a) equal order approximation without stabilization yields unphysical oscillations in pressure distribution, (b) equal order approximation with stabilization and different order approximation yield stable results and converges at optimal convergence rates in the respective norms. In case of equal order approximation, the effort involved is in the computation of stabilization parameter, whilst, in case of different order approximation, 14

(a)

(b)

(c)

Figure 10: Pressure plot for example 2: (a) analytical pressure (p = x2 + y 2 ), (b) without stabilization and (c) with stabilization. The domain is discretized with arbitrary polygons and equal order approximation for velocity and pressure.

15

additional effort is required: (a) to keep track of the velocity and pressure dofs and (b) to compute the serendipity barycentric coordinates. The serendipity coordinates need to be computed for each polygon. Extension of this work to the steady and the unsteady Navier-Stokes equation and its extensions to three dimensions will be a scope of future work. References References [1] N. Sukumar, A. Tabarrei, Polygonal interpolants: Construction and adaptive computations on quadtree meshes, European Congress on Computational methods in applied sciences and engineering (2004) 24–28. [2] C. Talischi, G. H. Paulino, Addressing integration error for polygonal finite elements through polynomial projections: A patch test connection, Mathematical Models and Methods in Applied Sciences 24 (2014) 1701–1727. [3] Y. Liu, W. Wang, B. L´evy, F. Sun, D. Yan, L. Lu, C. Yang, On centroidal voronoi tessellation - energy smoothness and fast computation, ACM Trans. Graph. 28 (2009) 101:1–101:17. [4] D. Yan, W. Wang, B. L´evy, Y. Liu, Efficient computation of clipped voronoi diagram for mesh generation, Computer-Aided Design 45 (2013) 843–852. [5] H. Ledoux, Computing the 3d voronoi diagram robustly: An easy explanation, in: Voronoi Diagrams in Science and Engineering, ISVD 2009. [6] R. Garimella, J. Kim, M. Berndt, Polyhedral mesh generation and optimization for non-manifold domains, in: Proceedings of the 22nd International Meshing Roundtable, pp. 313–330. [7] E. Wachspress, A Rational Finite Element Basis, volume 114, Academic Press, 1975. [8] N. Sukumar, Construction of polygonal interpolants: a maximum entropy approach, International Journal for Numerical Methods in Engineering 61 (2004) 2159–2181. [9] M. S. Floater, Mean value coordinates, Computer Aided Geometric Design 20 (2003) 19–27. [10] J. Bishop, A displacement based finite element formulation for general polyhedra using harmonic shape function, International Journal for Numerical Methods in Engineering 97 (2013) 1–31.

16

[11] D. Anisimov, D. Panozzo, K. Hormann, Blended barycentric coordinates, Computer Aided Geometric Design 52–53 (2017) 205–216. [12] Generalized Barycentric coordinates in Computer Graphics and Computational Mechanics, CRC Press, 2017. [13] N. Sukumar, A. Tabarrei, Conforming polygonal finite elements, International Journal for Numerical Methods in Engineering 61 (2004) 2045–2066. [14] S. Natarajan, E. T. Ooi, I. Chiong, C. Song, Convergence and accuracy of displacement based finite element formulation over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation, Finite Elements in Analysis and Design 85 (2014) 101–122. [15] A. Cangiani, V. Gyrya, G. Manzini, The nonconforming virtual element method for the Stokes equations, SIAM Journal on Numerical Analysis 54 (2016) 3411–3435. [16] K. Dai, G. Liu, T. Nguyen, An n-sided polygonal smoothed finite element method (nsfem) for solid mechanics, Finite Elements in Analysis and Design 43 (2007) 847–860. [17] A. Francis, A.Ortiz-Bernardin, S. P. A. Bordas, S. Natarajan, Linear smoothed polygonal and polyhedral finite elements, International Journal for Numerical Methods in Engineering 109 (2017) 1263–1288. [18] X. hai Tang, S.-C. Wu, C. Zheng, J. hai Zhang, A novel virtual node method for polygonal elements, Applied Mathematics and Mechanics 30 (2009) 1233–1246. [19] S. O. R. Biabanaki, A. R. Khoei, A polygonal finite element method for modeling arbitrary interfaces in large deformation problems, Computational Mechanics 50 (2012) 19–33. [20] C. Talischi, A. Pereira, G. H. Paulino, I. F. M. Menezes, M. S. Carvalho, Polygonal finite elements for incompressible fluid flow, International Journal for Numerical Methods in Engineering 74 (2014) 134–151. [21] S. Biabanaki, A. R. Khoei, P. Wriggers, Polygonal finite element methods for contact-impact problems on non-conformal meshes, Computer Methods in Applied Mechanics and Engineering 269 (2014) 198–221. [22] A. R. Khoei, R. Yasbolaghi, S. Biabanaki, A polygonal finite element method for modeling crack propogation with minimum remeshing, International Journal of Fracture 194 (2015) 123–148. [23] P. Antonietti, L Beir˜ao da Veiga, D. Mora, M. Verani, A stream virtual element formulation of the Stokes problem on polygonal meshes, SIAM Journal on Numerical Analysis 52 (2014) 386–404. 17

[24] L Beir˜ao da Veiga, C. Lovadina, G. Vacca, Divergence free virtual elements for the stokes problem on polygonal meshes, M2AN 51 (2017) 509–537. [25] E. C´aceres, G. N. Gatica, F. A. Sequeira, A Mixed Virtual Element Method for Quasi-Newtonian Stokes Flows, SIAM Journal on Numerical Analysis 56 (2018) 1210–1242. [26] L Beir˜ao da Veiga, C. Lovadina, G. Vacca, Virtual elements for the Navier-Stokes problem on polygonal meshes, SIAM Journal on Numerical Analysis 56 (2018) 317–343. [27] L. Chen, F. Wang, A divergence free weak virtual element method for the Stokes problem on polytopal meshes, Journal of Scientific Computing 78 (2019) 864–886. [28] T. Vu-Huu, C. Le-Thanh, H. Nguyen-Xuan, M. Abdel-Wahab, Incompressible fluid computation based on polygonal finite element, in: Proceedings of the 1st International Conference on Numerical Modelling in Engineering. [29] H. Chi, C. Talischi, O. Lopez-Pamies, G. H. Paulino, Polygonal finite elements for finite elasticity, International Journal for Numerical Methods in Engineering 101 (2015) 305–328. [30] H. Chi, C. Talischi, O. Lopez-Pamies, G. H. Paulino, A paradigm for higherorder polygonal elements in finite elasticity using a gradient correction scheme, Computer Methods in Applied Mechanics and Engineering 306 (2016) 216–251. [31] T.-P. Fries, H. G. Matthies, A Review of Petrov-Galerkin Stabilization Approaches and an Extension to Meshfree Methods, Technical Report, Institute of Scientific Computing, Technical University Braunschweigh, Brunswick, Germany, 2004. [32] J. Warren, On the uniqueness of barycentric coordinates, in: Proceedings of AGGM02, pp. 93–99. [33] J. Warren, S. Schaefer, A. Hirani, M. Desbrun, Barycentric coordinates for convex sets, Advances in Computational Mechanics 27 (2007) 319–338. [34] M. S. Floater, A. Gillette, N. Sukumar, Gradient bounds for Wachspress coordinates in polytopes, SIAM J. Numer. Anal. 52 (2014) 515–532. [35] A. Rand, A. Gillette, C. Bajaj, Quadratic serendipity finite elements on polygons using generalized barycentric coordinates, Mathematics of Computation 83 (2014) 2691–2716. [36] A. Sinu, S. Natarajan, K. Shankar, Quadratic serendipity finite elements over convex polyhedra, International Journal for Numerical Methods in Engineering 113 (2018) 109–129.

18

[37] N. Sukumar, E. A. Malsch, Recent advances in the construction of polygonal finite element interpolants, Archives of Computational Methods in Engineering 13 (2006) 129–163. [38] A. Sommariva, M. Vianello, Product gauss cubature over polygons based on green’s integration formula, BIT Numerical Mathematics 47 (2007) 441–453. [39] S. Natarajan, S. P. A. Bordas, D. R. Mahapatra, Numerical integration over arbitrary polygonal domains based on schwarz-christoffel conforming mapping, International Journal for Numerical Methods in Engineering 80 (2009) 103–134. [40] S. Mousavi, H. Xiao, N. Sukumar, Generalized Gaussian quadrature rules on arbitrary polygons, International Journal for Numerical Methods in Engineering 82 (2010) 99–113. [41] S. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons, Computational Mechanics 47 (2011) 535–554. [42] C. Lee, H. Kim, S. Im, Polyhedral elements by means of node/edge-based smoothed finite element method, International Journal for Numerical Methods in Engineering (2016). [43] S. Natarajan, S. P. A. Bordas, E. T. Ooi, Virtual and smoothed finite elements: a connection and its application to polygonal/polyhedral finite element methods, International Journal for Numerical Methods in Engineering 104 (2015) 1173–1199. [44] A. Francis, S. Natarajan, E. Atroshchenko, B. L´evy, S. P. Bordas, A one point integration rule over star convex polytopes, Computers and Structures 215 (2019) 43–64. [45] C. Talischi, G. H. Paulino, A. Pereira, I. F. Menezes, Polytop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes, Struct. Multidisc Optim 45 (2012) 329–357. [46] Z. Cai, C. Tong, P. S. Vassilevski, C. Wang, Mixed finite element methods for incompressible flow: stationary stokes equations, Numerical Methods for Partial Differential equations 26 (2009) 957–978. [47] B. S. Hosseini, M. M¨oller, S. Turek, Isogeometric analysis of the navier-stokes equations with taylor-hood b-spline elements, Computer Methods in Applied Mechanics and Engineering 267 (2015) 264–281. [48] F. Brezzi, K. Bathe, A discourse on the stability conditions for mixed finite element formulations, Computer Methods in Applied Mechanics and Engineering 82 (1990) 27–57. 19

[49] D. Chapelle, K. Bathe, The inf-sup test, Computers & Structures 47 (1993) 537–545. [50] K. J. Bathe, The inf-sup condition and its evaluation for mixed finite element methods, Computers and Structures 79 (2001) 243–252.

20