Arbitrary order 2D virtual elements for polygonal meshes: Part I, elastic problem

TL;DR

This paper develops a Virtual Element Method (VEM) for 2D elastic problems, enabling flexible polygonal elements with higher-order approximations, and compares its performance with classical finite element methods.

Contribution

It introduces a standardized procedure for constructing VEM terms for polygonal elements with higher-order interpolation, expanding the applicability of VEM in structural analysis.

Findings

VEM performs comparably to classical FEM in numerical tests.

Higher-order polygonal elements improve approximation accuracy.

The method effectively handles various polygonal shapes and edge counts.

Abstract

The present work deals with the formulation of a Virtual Element Method (VEM) for two dimensional structural problems. The contribution is split in two parts: in part I, the elastic problem is discussed, while in part II [3] the method is extended to material nonlinearity, considering different inelastic responses of the material. In particular, in part I a standardized procedure for the construction of all the terms required for the implementation of the method in a code is explained. The procedure is initially illustrated for the simplest case of quadrilateral virtual elements with linear approximation of displacement variables on the boundary of the element. Then, the case of polygonal elements with quadratic and, even, higher order interpolation is considered. The construction of the method is detailed, deriving the approximation of the consistent term, the required stabilization term and the loading term for all the considered virtual elements. A wide numerical investigation is performed to assess the performances of the developed virtual elements, considering different number of edges describing the elements and different order of approximations of the unknown field. Numerical results are also compared with the one recovered using the classical finite element method.

Tables

Table 1. Table 1 : Indexing s 𝑠 s of vector 𝕧 ~ h subscript ~ 𝕧 ℎ \tilde{{\mathbb{v}}}_{h} components ( u ~ h , s (\tilde{u}_{h,s}   v ~ h , s ) \tilde{v}_{h,s}) for the case with general order k 𝑘 k , and polygon with m 𝑚 m edges.

	${\tilde{u}}_{h,s}$		${\tilde{v}}_{h,s}$
	$\left(s=2i-1\right)$		$\left(s=2i\right)$
Vertex $\left(i\right)$		$i=1,\dots ,m$
Edge $\left(i\right)$		$i=m+1,\dots ,km$
Internal $\left(i\right)$		$i=km+1,\dots ,km+r$

Equations

$$\mathbf{\epsilon }(u)=\mathbb{S}u\text{with}\mathbb{S}=\left[\begin{array}{cc}{\partial }_x& 0\\ 0& {\partial }_y\\ {\partial }_y& {\partial }_x\end{array}\right],$$

$$\mathbf{\sigma }=\mathbb{C}\mathbf{\epsilon }$$

$$\{\begin{array}{rl}& \text{Find u∈V such that}\\ & a(u,v)=<b,v>\forall v\in \mathbf{V}\end{array}$$

$$a(u,v)$$

$$<b,v>$$

$$\displaystyle V_{h|E}=\big{\{}$$

$$\displaystyle v_{h}|_{e}\in\mathsf{P}_{k}(e)\quad\forall e\in\partial E\big{\}},$$

$${\mathbf{V}}_{h|E}={\left[V_{h|E}\right]}^2.$$

$$n=\text{dim}({\mathbf{V}}_{h|E})=2mk+k(k-1),$$

$${\mathbf{V}}_h=\{v\in \mathbf{V}:v{|}_E\in {\mathbf{V}}_{h|E}\forall E\in {\mathcal{T}}_h\}.$$

$$\{\begin{array}{rl}& \text{mbox}find{u}_h\in {\mathbf{V}}_h\text{mbox}suchthat\\ & a_h(u_h,v_h)=<b_h,v_h>\forall v_h\in {\mathbf{V}}_h,\end{array}$$

$$a_h(u_h,v_h)=\sum _{E\in {\mathcal{T}}_h}{a}_h^E(u_h,v_h)\forall u_h,v_h\in {\mathbf{V}}_h.$$

$$\Pi :$$

$$v_h\mapsto \Pi (v_h),$$

$${\int }_E\Pi (v_h{)}^T{\mathbf{\epsilon }}^P={\int }_E\mathbf{\epsilon }(v_h{)}^T{\mathbf{\epsilon }}^P,\forall {\mathbf{\epsilon }}^P\in {\mathsf{P}}_{k-1}(E{)}_{\text{sym}}^{2\times 2}.$$

$$a_h^E(u_h,v_h)={\int }_E{\left[\Pi (u_h)\right]}^T\mathbb{C}\Pi (v_h)+S^E\left(u_h,v_h\right)$$

$$\xi =\frac{x-x_c}{h_E},\eta =\frac{y-y_c}{h_E},$$

$$n=\text{dim}({\mathbf{V}}_{h|E})=4\times 2=8$$

$$\ell =\text{dim}({\mathsf{P}}_0(E{)}_{\text{sym}}^{2\times 2})=3,$$

$${\mathsf{P}}_0(E{)}_{\text{sym}}^{2\times 2}=\text{span}\left\{\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\right\}.$$

$$\{\begin{array}{rl}& v_h={\mathbb{N}}^V{\tilde{v}}_h,\\ & {\mathbf{\epsilon }}^P={\mathbb{N}}^P\hat{\mathbf{\epsilon }},\end{array}$$

$${\mathbb{N}}^P=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$$

$$\Pi (v_h)={\mathbb{N}}^P{\Pi }^{\text{m}}{\tilde{v}}_h$$

$${\int }_E{\left({\mathbb{N}}^P{\Pi }^{\text{m}}{\tilde{v}}_h\right)}^{\text{T}}{\mathbb{N}}^P\hat{\mathbf{\epsilon }}={\int }_E{\left[\mathbf{\epsilon }({\mathbb{N}}^V{\tilde{v}}_h)\right]}^{\text{T}}{\mathbb{N}}^P\hat{\mathbf{\epsilon }}\forall \hat{\mathbf{\epsilon }}\in {\mathbb{R}}^3.$$

$${\int }_E{\left({\mathbb{N}}^P{\Pi }^{\text{m}}{\tilde{v}}_h\right)}^{\text{T}}{\mathbb{N}}^P\hat{\mathbf{\epsilon }}={\int }_{\partial E}{\left({\mathbb{N}}^V{\tilde{v}}_h\right)}^{\text{T}}\left({\mathbb{N}}_E{\mathbb{N}}^P\hat{\mathbf{\epsilon }}\right)$$

$${\mathbb{N}}_E=\left[\begin{array}{ccc}{n}_x& 0& n_y\\ 0& n_y& n_x\end{array}\right]$$

$${\hat{\mathbf{\epsilon }}}^{\text{T}}\mathcal{G}{\Pi }^{\text{m}}{\tilde{v}}_h={\hat{\mathbf{\epsilon }}}^{\text{T}}\mathcal{B}{\tilde{v}}_h\forall \hat{\mathbf{\epsilon }}\in {\mathbb{R}}^3,$$

$$\mathcal{G}={\int }_E{\left({\mathbb{N}}^P\right)}^{\text{T}}{\mathbb{N}}^P\in {\mathbb{R}}^{3\times 3}$$

$$\mathcal{B}={\int }_{\partial E}{\left({\mathbb{N}}_E{\mathbb{N}}^P\right)}^{\text{T}}{\mathbb{N}}^V\in {\mathbb{R}}^{3\times 8}.$$

$${\Pi }^{\text{m}}={\mathcal{G}}^{-1}\mathcal{B},$$

$${\int }_E{\left(\Pi (v_h)\right)}^{\text{T}}\mathbb{C}\Pi (u_h)={\int }_E{\left({\mathbb{N}}^P{\mathcal{G}}^{-1}\mathcal{B}{\tilde{v}}_h\right)}^{\text{T}}\mathbb{C}\left({\mathbb{N}}^P{\mathcal{G}}^{-1}\mathcal{B}{\tilde{u}}_h\right)={\tilde{v}}_h^{\text{T}}{\mathbb{K}}_{\text{c}}{\tilde{u}}_h$$

$${\mathbb{K}}_{\text{c}}={\mathcal{B}}^{\text{T}}{\mathcal{G}}^{\text{-T}}\left[{\int }_E{\left({\mathbb{N}}^P\right)}^{\text{T}}\mathbb{C}{\mathbb{N}}^P\right]{\mathcal{G}}^{-1}\mathcal{B}.$$

$$({\mathsf{P}}_1(E){)}^2=\text{span}\left\{\left(\begin{array}{c}1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\end{array}\right),\left(\begin{array}{c}\xi \\ 0\end{array}\right),\left(\begin{array}{c}0\\ \xi \end{array}\right),\left(\begin{array}{c}\eta \\ 0\end{array}\right),\left(\begin{array}{c}0\\ \eta \end{array}\right)\right\},$$

Full text

1 Introduction

The Virtual Element Method (VEM) was recently introduced in volley ; hitchhikers as a generalization of the Finite Element Method (e.g. Bathe1996 ; Hughes:book ; Wriggers:book ) with the additional capability of dealing with very general polygonal/polyhedral meshes and the possibility to easily implement highly regular discrete spaces (see for instance Brezzi:Marini:plates ). By making use also of non-polynomial shape functions, but avoiding the explicit construction of the local basis functions, the VEM can easily handle general polygons/polyhedrons without complex integrations on the element. Polytopal meshes can be very useful for a wide range of reasons, including meshing of the domain (such as cracks) and data (such as inclusions) features, automatic use of hanging nodes, use of moving meshes, adaptivity. Recently, the interest in numerical methods that can make use of general polytopal meshes has undergone a significant growth in the mathematical and engineering literature. Here, we limit to citing some work in the area of polygonal Finite Elements ST04 applied to advanced problems in structural mechanics, that is Paulino-nonlinear-polygonal ; TPPM10 ; Wriggers-2014 ; Paulino-cracks , see also DiPietro-Ern-1 .

The Virtual Element Method has strongly developed and is now able to cope with a wide range of topics and problems, see for instance CH-VEM and references therein. In the more specific framework of structural mechanics, VEM has been introduced in VEM-elasticity for (possibly incompressible) two dimensional linear elasticity and general “polynomial” order, in Paulino-VEM for three dimensional linear elasticity and lowest order, in BeiraoLovaMora for general two dimensional elastic and inelastic problems under small deformations (lowest order). More specific applications have been also considered, such as contact problems wriggers , topology optimization Topology-VEM and geomechanics Andersen-geo .

The present contribution is the first part of a work made of two papers, that take the steps from BeiraoLovaMora ; VEM-elasticity and extend the VEM to the case of arbitrary order of accuracy (or “polynomial” order), in the framework of general elastic and inelastic small deformation 2D problems. In brief, we are able to develop a numerical scheme of arbitrary order of accuracy that works for nonlinear deformation problems, general polygonal meshes and still yields a conforming solution. The key idea of the method is to use a projection of the strains on a polynomial space, apply the constitutive law to this approximated strain at specific points (which may be seen as the analogous of the Gauss points in standard FEM), and finally stabilize the ensuing discrete stiffness operator by standard techniques in the VEM literature.

This first part of the work is focused on linear elasticity and sets all the basic tools necessary for the more advanced problems considered in part two. Namely, we introduce the Virtual Element space, the associated degrees of freedom and projectors, and the approximation of the elasticity stiffness matrix. When compared to the VEM in VEM-elasticity , the present scheme has strong similarities (both handle arbitrary order and linear elasticity). Nevertheless, in this paper there are three important additional issues:

Differently from VEM-elasticity , the presented method is already suitable for generalization to nonlinear problems (treated in Part II), without any change in the foundations of the scheme. 2. 2.

An important contribution of this paper is that it gives extensive guidelines on the construction and coding of the proposed scheme, starting from the simple case of lowest order on quadrilaterals and going up (step by step) to general order and general polygons. Although it already exists a well written coding guide for VEM hitchhikers , that work is more focused on the diffusion problem rather than elasticity; moreover, here we propose a coding approach that is more suitable for engineers. The reader interested in the coding of VEM will surely benefit from reading both works; 3. 3.

In this paper we also present, for the first time, a deep numerical testing of a VEM method of general order in the realm of linear elasticity, with highly encouraging results.

The paper is organized as follows. In Section 2 we present the continuous problem and the general structure of the Virtual Element formulation. In Section 3 we thoroughly describe the construction of the method. To help the reader understand all the details, we start with the lowest order case (linear accuracy) on quadrilaterals in Section 3.1, then consider the quadratic case on general polygons in Section 3.2, and only as a last step we consider the fully general case, see Section 3.3. Finally, in Section 4 we present a set of numerical tests for different orders and different families of meshes.

2 General structure of the Virtual Element formulation

In the present section we introduce the continuous problem and give an overview on the proposed Virtual Element discretization. More details and specific explanations will be found in the following sections.

2.1 The continuous problem

In this section we present the linear elastic problem in the two-dimensional (2D) framework, under the assumption of infinitesimal strain and small displacements. In the following the Voigt notation is adopted, so that stress and strain tensors are represented as $3-$component vectors, and the fourth-order constitutive tensor is represented as a $3\times 3$ matrix.

Let $\Omega$ be a continuous body occupying a region of the two-dimensional space $\mathscr{R}^{2}$ in which the Cartesian coordinate system $(\textrm{O},x,y)$ is introduced. The displacement field is denoted by the vector ${\mathbb{u}}(x,y)=\left(u\;v\right)^{\textrm{T}}$ and the associated strain defined as:

[TABLE]

where $\partial_{(\bullet)}$ denotes the partial derivative operator with respect to the $(\bullet)$-coordinate.

The linear elastic constitutive law is considered for the body $\Omega$, so that the stress is given by the relationship:

[TABLE]

where $\mathbb{C}=\mathbb{C}(x)$ is the $3\times 3$ (uniformly positive) elastic matrix, possibly depending on the position vector ${\bf x}=(x\;y)^{\textrm{T}}\in\Omega$.

The body $\Omega$ is subjected to distributed volume forces $\mathbb{b}$. Without loosing in generality, for simplicity it is assumed that the displacements are vanishing on the whole boundary of $\Omega$ (see Remark 1 at the end of this section). The variational formulation of the elastostatic problem is provided by the virtual work principle:

[TABLE]

where ${\bf V}:=(H^{1}_{0}(\Omega))^{2}$ is the space of the admissible displacement field and

[TABLE]

The form $a(\cdot,\cdot)$ is symmetric, continuous and coercive on ${\bf V}$, so that problem (3) is well posed.

2.2 The Virtual Element formulation

We start by presenting the discrete (virtual) space of admissible displacements, that is the same one introduced in VEM-elasticity . See also volley ; hitchhikers for a scalar version of the space used for the diffusion problem.

Let ${\mathscr{T}}_{h}$ be a simple polygonal mesh on $\Omega$. This can be any decomposition of $\Omega$ in non overlapping polygons $E$ with straight edges. Let ${\cal E}_{h}$ denote the set of the edges of ${\mathscr{T}}_{h}$. The symbol $m$ represents the number of edges of a polygon $E$,and the typical edge of the polygon $E$ is indicated by $e$, (i.e. $e\in\partial E$). A sample element $E$, with an ordering of vertexes and edges, is depicted in figure 1. The space ${\bf V}_{h}$ will be defined element-wise, by introducing local spaces ${{\bf V}_{h|E}}$ and the associated local degrees of freedom, as in standard Finite Element (FE) analysis. Instead, differently from standard FE, the definition of the local spaces ${{\bf V}_{h|E}}$ is not fully explicit.

Let $k$ be a positive integer, representing the “degree of accuracy” of the method. Then, given an element $E\in{\mathscr{T}}_{h}$, we start by defining the scalar space:

[TABLE]

where $\mathsf{P}_{k}(E)$ denotes the space of polynomials of degree (up to) $k$ defined on $E$, with the agreement that $\mathsf{P}_{\!-1}(E)=\{0\}$. The proposed virtual displacement space is

[TABLE]

The space ${{\bf V}_{h|E}}$ is made of vector valued functions ${\mathbb{v}}_{h}$ such that VEM-elasticity ; volley :

• •

${\mathbb{v}}_{h}$ is a polynomial of degree $\leq k$ on each edge $e$ of $E$, i.e. ${\mathbb{v}}_{h|e}\in(\mathsf{P}_{k}(e))^{2}$;

• •

${\mathbb{v}}_{h}$ is globally continuous on $\partial E$;

• •

{\bf div}\big{(}\mathbb{C}\boldsymbol{\varepsilon}({\mathbb{v}}_{h})\big{)} is a polynomial of degree $\leq k-2$ in $E$.

The following important observations hold:

• •

the functions ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$ are explicitly known on $\partial E$;

• •

the functions ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$ are not explicitly known inside the element $E$;

• •

it holds $(\mathsf{P}_{k}(E))^{2}\subseteq{{\bf V}_{h|E}}$ (that is important for approximation).

A sample polygon with $m=5$ edges is represented in Figure 1 with indication of vertexes and edges. Taking into account that the unknown vector ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$ has two components, for any $E\in{\mathscr{T}}_{h}$, the degrees of freedom for ${{\bf V}_{h|E}}$ will be:

• •

$2m$ pointwise values ${\mathbb{v}}_{h}(\nu_{i})$ at corners $\nu_{i}$ of $E$, $i=1,2,..,m$,

• •

$2m(k-1)$ pointwise values ${\mathbb{v}}_{h}(\boldsymbol{x}^{e}_{j})$ at edge internal nodes $\{\boldsymbol{x}^{e}_{j}\}$ (i.e. not comprising extrema of edge $e$), $j=1,...,k-1$, being $m$ the number of edges of a polygon;

• •

$2\frac{k(k-1)}{2}$ scalar moments of the unknown field over the element, not associated with a specific location over $E$ (these degrees of freedom will be specified in the following).

Sample figures for the degrees of freedom will be shown in the following sections, where more details are given. The dimension of the space ${{\bf V}_{h|E}}$ results

[TABLE]

where $2m$ degrees of freedom are associated with the $m$ corners of $\partial E$, $2m(k-1)$ degrees of freedom are associated with the $m$ edges (with $(k-1)$ internal nodal points per edge), and $k(k-1)$ degrees of freedom are the moments. A local “Lagrange-type” basis for the space ${{\bf V}_{h|E}}$ is chosen as follows: the $i-th$ function of the local basis is given by the unique function in ${{\bf V}_{h|E}}$ that has unit value on the $i-th$ degree of freedom and vanishes for the remaining ones.

As in standard FE methods, the global space ${\bf V}_{h}\subseteq{\bf V}$ is built by assembling the local spaces ${{\bf V}_{h|E}}$ as usual:

[TABLE]

A virtual element method (VEM) for the above problem is constructed following a procedure that closely resembles the Finite Element approach, i.e. by restricting the original variational formulation (3) to the discrete space ${\bf V}_{h}$ and approximating the ensuing terms:

[TABLE]

where

• •

${\bf V}_{h}\subset{\bf V}$ is the virtual space introduced above;

• •

$a_{h}(\cdot,\cdot)\>:\>{\bf V}_{h}\times{\bf V}_{h}\rightarrow{\mathbb{R}}$ is a discrete bilinear form approximating the continuous form $a(\cdot,\cdot)$;

• •

$<\mathbb{b}_{h},{\mathbb{v}}_{h}>$ is the term approximating the virtual work of the external load.

The discrete bilinear form is built element by element, assuming:

[TABLE]

On the other hand, since the functions of the space ${\bf V}_{h}$ are not know explicitly inside the elements, the bilinear forms $a_{h}^{E}(\cdot,\cdot)$ cannot be evaluated by standard Gauss integration. The construction of the local discrete form $a_{h}^{E}(\cdot,\cdot)$, a key point in the VEM, is performed following the procedure described in the following.

Since neither the function ${\mathbb{v}}_{h}$ nor its gradient are explicitly computable in the element interior points, the method proceeds by introducing a projection operator $\Pi$, representing the approximated strain associated with the virtual displacement, defined as:

[TABLE]

where, albeit the adopted notation, the $2\times 2$ symmetry character of the strain tensor is remarked. It is noted that $\ell=\textrm{dim}(\mathsf{P}_{k-1}(E)^{2\times 2}_{\textrm{sym}})=2\frac{k(k+1)}{2}$.

Given ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, such an operator $\Pi$ is defined as the unique function $\Pi({\mathbb{v}}_{h})\in\mathsf{P}_{k-1}(E)^{2\times 2}_{\textrm{sym}}$ that satisfies the condition:

[TABLE]

This operator represents the best approximation of the strains (in the square integral norm) in the space of piecewise polynomials of degree $(k-1)$. Although the functions in ${{\bf V}_{h|E}}$ are virtual, the right hand side in (11) (and thus the operator $\Pi$) turns out to be computable with simple calculations, as we will show in the next section.

Once computed, the $\Pi$ operator approximates the internal work part associated to the equilibrium weak formulation as follows:

[TABLE]

with $S^{E}(\cdot,\cdot)$ a suitable stabilizing term, needed to preserve the coercivity of the system, and described in the sequel. The numerical approximation of the loading term is presented in the next section.

Remark 1

Since on the boundary of the elements the functions of ${\bf V}_{h}$ are piecewise polynomials, the imposition of more general boundary conditions follows exactly the same procedure as in standard finite elements. For instance, non-homogeneous Dirichlet boundary conditions are enforced by standard interpolation of the prescribed displacements, while Neumann boundary conditions are implemented by adding the usual boundary integral to the formulation.

3 Construction of the method

In the present section we detail the construction of the proposed virtual scheme. In order to make things clearer, our presentation will be subdivided into three parts. Initially, we start describing the simple case $k=1$ for quadrilateral elements $E$. Afterwards, we describe the case $k=2$ for general polygons $E$ with $m$ edges. Last, we address the case of general polynomial order $k$.

Remark 2

In the following, we adopt local (to the current element $E$) scaled space coordinates, such that: (1) the origin of the coordinate axes is the centroid of $E$; (2) all the coordinate values are scaled by the diameter $h_{E}$ of the element $E$. For instance, if the position of a point $P$ is $(x\;y)$ in standard cartesian coordinates, then its local scaled coordinates $(\xi\;\eta)$ are given by

[TABLE]

where $(x_{c}\;y_{c})$ denotes the cartesian coordinates of the centroid of $E$.

3.1 The case $k=1$ for quadrilaterals

We start considering the case of $E$ being a quadrilateral polygon, i.e. $m=4$, and assuming $k=1$; in this case, the degrees of freedom for ${{\bf V}_{h|E}}$ are the only displacement field values at the corners of the element, with no edge degrees of freedom and no moments. The degrees of freedom are depicted in figure 2.

Therefore, the virtual element space and the polynomial space of approximated strains, on the basis of Eq. (10), satisfy the following conditions:

[TABLE]

so that it is set

[TABLE]

Formally, by representing the virtual displacement field and the approximated strains in terms of the basis functions, we can write

[TABLE]

where ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, $\tilde{{\mathbb{v}}}_{h}\in\mathbb{R}^{8}$, $\boldsymbol{\varepsilon}^{P}\in\mathsf{P}_{0}(E)^{2\times 2}_{\textrm{sym}}$, $\hat{\boldsymbol{\varepsilon}}\in\mathbb{R}^{3}$. The matrix $\mathbb{N}^{P}$ is given by

[TABLE]

Matrix $\mathbb{N}^{V}$ can be viewed as the shape function operator of a typical displacement-based 2D finite element, even if its entries are not analytical over the element interior. In particular, it has the same structure, i.e. odd-index columns are associated with $u_{h}$ and have zero second entries, whereas even-index columns are associated with $v_{h}$ and have zero first entries. Since the only degrees of freedom are the vertex values and recalling that the displacements have two components, $\mathbb{N}^{V}$ is a $2\times 8$ matrix operator. Ordering the four vertexes of $E$ with $i=1,...4$, for $s=2i-1$, the first entry of the $s^{th}$-column of $\mathbb{N}^{V}$ is formally given by the unique function of $V_{h|E}$ that takes the value $1$ at vertex $i$ and vanishes at all other vertexes, the second entry being [math]. Analogously, for $s=2i$, the $s^{th}$-column first entry is zero; the second one is formally given by the unique function of $V_{h|E}$ that takes the value $1$ at vertex $i$ and vanishes at all other vertexes. The operator $\mathbb{N}^{V}$ is key for understanding the present construction, but is actually not computable in the interior of the element. On the other hand, since the restrictions of the functions of ${{\bf V}_{h|E}}$ to $\partial E$ are piecewise linear and continuous, the computation of $\mathbb{N}^{V}$ can be carried out easily on the element boundary.

3.1.1 Consistent term

In the present section we detail the construction of the consistent part of the bilinear form $a^{E}_{h}(\cdot,\cdot)$, i.e the first term appearing in the right hand side of (12). We start by the practical construction of the projection operator $\Pi$. We denote by $\Pi^{\textrm{m}}$ the operator $\Pi$ expressed as a matrix, in terms of the basis of ${{\bf V}_{h|E}}$ and $\mathsf{P}_{0}(E)^{2\times 2}_{\textrm{sym}}$. In other words, the projected strain results in

[TABLE]

with $\Pi^{\textrm{m}}\in\mathbb{R}^{3\times 8}$. The following steps lead to the computation of the matrix operator $\Pi^{\textrm{m}}$.

With the above definitions, (11) becomes:

[TABLE]

Taking into account Eq. (1) and noting that the internal contribution vanishes since $\boldsymbol{\varepsilon}^{P}\!=\!\mathbb{N}^{P}\!\!\hat{\boldsymbol{\varepsilon}}$ is constant on the element, an integration by parts on the right hand side yields

[TABLE]

for any $\hat{\boldsymbol{\varepsilon}}$ in $\mathbb{R}^{3}$, where the matrix

[TABLE]

is built by the components of the outward normal ${\bf n}_{E}=(n_{x},n_{y})^{T}$ to the element edges.

We can re-write Eq. (20) more compactly as:

[TABLE]

where:

[TABLE]

Remark 3

Note that, in light of the previous observations, also the matrix $\mathscr{B}$ is computable, since it only requires knowledge of $\mathbb{N}^{V}$ on the boundary of the element.

Eq. (21) provides the required operator $\Pi^{\textrm{m}}$ as the solution of a linear system:

[TABLE]

where $\mathscr{G}$ is clearly symmetric and positive definite.

Combining eqs. (18) and (24), from (12), one obtains:

[TABLE]

which defines explicitly the first part of the stiffness matrix:

[TABLE]

As already noted in (12), the above term is not sufficient for the stability of the method, and needs to be completed with a stabilization term that is discussed in the next section.

3.1.2 Stabilization term

The discrete bilinear form cannot be composed of the consistent term alone, otherwise the coercivity of the system may be lost and hourglass modes could arise. The presence of a stability term is therefore standard for the Virtual Element Method, as illustrated for instance volley ; Brezzi:Marini:plates ; VEM-elasticity ; hitchhikers ; BeiraoLovaMora ; CH-VEM for a few examples. Herein, we first detail the construction of the stability term, commenting its motivations at the end of the section.

Let a basis for the space $(\mathsf{P}_{1}(E))^{2}$ be chosen as:

[TABLE]

The elements $\mathbb{p}\in(\mathsf{P}_{1}(E))^{2}$ can be written in terms of the above basis as vectors $\hat{\mathbb{p}}\in\mathbb{R}^{6}$. Moreover, since $(\mathsf{P}_{1}(E))^{2}\subseteq{{\bf V}_{h|E}}$, the elements $\mathbb{p}\in(\mathsf{P}_{1}(E))^{2}$ can be written also in terms of the basis of ${{\bf V}_{h|E}}$. We call ${\bf D}$ the matrix associated to such change of basis, that can be simply calculated by evaluating the six functions in (27) through the degrees of freedom of ${{\bf V}_{h|E}}$. Specifically, one has

[TABLE]

where $\tilde{\mathbb{p}}\in\mathbb{R}^{8}$ represents $\mathbb{p}$ written in terms of the ${{\bf V}_{h|E}}$ basis. The matrix ${\bf D}$ is given by

[TABLE]

where the $i-th$ vertex of the polygon has scaled coordinates $(\xi_{i},\eta_{i})$, $i=1,...,4$.

A stabilization strategy is obtained by considering the form $S^{E}({\mathbb{u}}_{h},{\mathbb{v}}_{h})=\tilde{{\mathbb{u}}}_{h}^{T}\mathbb{K}_{S}\tilde{{\mathbb{v}}}_{h}$ for all ${\mathbb{u}}_{h}\;,{\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, taking:

[TABLE]

with $\tau$ positive real number (see below). The total stiffness operator results

[TABLE]

The role of the stabilization matrix $\mathbb{K}_{S}$ is to keep the positivity of the discrete (local) energy form (up to the constant functions that should still give zero energy). Indeed, it can be checked that, thanks to the contribution of the matrix $\mathbb{K}_{S}$, the discrete bilinear form satisfies the condition:

[TABLE]

for any ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$ such that $a^{E}_{h}({\mathbb{v}}_{h},{\mathbb{v}}_{h})=0$. Thus, a vanishing discrete elastic energy implies a rigid body motion.

The trace term in the definition of $\mathbb{K}_{S}$ is added in order to guarantee the correct scaling of the energy with respect to the element size and material constants. The factor $\tau\in\mathbb{R}$, $\tau>0$, can be discarded but is included in order to allow for generalizations. Indeed, a typical good choice is simply $\tau=1$, but other choices can be found in the literature (see for instance Paulino-VEM ). As shown later in the numerical tests, the method is not much sensible to the parameter $\tau$. We refer to the initial paper volley for more details on the motivations of the stabilization term, while a comparison with the standard “single Gauss point” quadrilateral FEM element and the related hourglass mode control can be found in hourglass .

Remark 4

It is straightforward to check that, if $\tilde{\mathbb{p}}\in{{\bf V}_{h|E}}$ represents a polynomial (in other words satisfies (28) for some $\hat{\mathbb{p}}\in\mathbb{R}^{6}$) then $\mathbb{K}_{S}\>\tilde{\mathbb{p}}$ vanishes. Therefore, from the definition of $\mathbb{K}$, it can be derived that (in the case of constant elastic coefficients) the local patch test is satisfied by the discrete bilinear form, i.e.:

[TABLE]

a property that in the Virtual Element terminology is called $k$-consistency volley ; hitchhikers .

3.1.3 Loading term

The loading term for the case $k=1$ can be approximated by applying an integration rule based on vertexes. As in standard FE technology, the wirtual work exerted by the external loads is the element-wise sum of local contributions:

[TABLE]

where we remind that $<\mathbb{b}_{h},{\mathbb{v}}_{h}>$, see also (8), represents the VEM approximation of the right hand side, expressed below (for varying ${\mathbb{v}}_{h}$) as a linear operator on ${\bf V}_{h}$.

Let the vector $\overline{\mathbb{b}}\in\mathbb{R}^{2}$ be defined component-wise

[TABLE]

with $b_{1},b_{2}$ the two components of the volume force $\mathbb{b}$. Then, we can approximate the local loading term $\int_{E}\mathbb{b}^{\textrm{T}}{\mathbb{v}}_{h}$, for all ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, with the following simple integration rule

[TABLE]

where $(\tilde{{\mathbb{v}}}_{h})_{s}$, $s=1,...,8$, is the $s-th$ component of the vector $\tilde{{\mathbb{v}}}_{h}$, evaluated at the four corners $\nu_{i}$, $i=1,...,4$, of $E$.

The loading term approximation reported in Eq. (33) is exact whenever the load $\mathbb{b}$ is constant over the element, and the displacement field ${\mathbb{v}}_{h}$ is a first order polynomial.

3.2 The case $k=2$ for general polygons

The step from $k=1$ to $k=2$ is substantial, since in the case $k=2$ edge degrees of freedom and moments come into play in the approximation of the displacement field. We now consider the case of a general polygon with $m\geq 3$ edges.

By referring to the list of degrees of freedom in Section 2, it is immediate to check that for $k=2$ there are

• •

2 (since the functions are vector valued) degrees of freedom per vertex, representing the displacement value at the vertex

• •

2 degrees of freedom per edge (representing the displacement value at the midpoint of the edge)

• •

2 scalar moments over the element

The degrees of freedom are depicted in figure 3. Note that the edge midpoints are denoted by ${\boldsymbol{x}}_{1}^{e}$, $e\in\partial E$.

The moments degrees of freedom are the component-wise average of the function ${\mathbb{v}}_{h}$

[TABLE]

The virtual element space and the polynomial space of approximated strains, from (6) with $k=2$, satisfy:

[TABLE]

and we set

[TABLE]

Representation (16) still applies: ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, $\tilde{{\mathbb{v}}}_{h}\in\mathbb{R}^{4m+2}$, $\boldsymbol{\varepsilon}^{P}\in\mathsf{P}_{1}(E)^{2\times 2}_{\textrm{sym}}$, $\hat{\boldsymbol{\varepsilon}}\in\mathbb{R}^{9}$. The matrix $\mathbb{N}^{P}$ is given by

[TABLE]

The matrix operator $\mathbb{N}^{V}$ representing the basis functions of ${{\bf V}_{h|E}}$ is a $2\times(4m+2)$ operator. Ordering the $m$ vertexes and the $m$ edge midpoint nodes of $E$ with $i=1,..,2m$, setting $s=2i-1$, the $s^{th}$-column first entry is given by the unique function of $V_{h|E}$ that takes the value $1$ at $i-th$ location, vanishes at all other vertexes and midside nodes and has null average over the polygon (cf. (34)); the second entry is zero. Analogously, setting $s=2i$, the $s^{th}$-column has null first entry, being the second entry a basis function with the above mentioned properties. The subsequent $(4m+1)-th,\;(4m+2)-th$ columns contain basis functions associated with moments degrees of freedom, hence they vanish at vertexes and edge nodes, and have unit average over the element. Similarly to the $k=1$ case, we note that for $k=2$ the functions of ${{\bf V}_{h|E}}$ are piecewise quadratic and continuous on the boundary of $E$; therefore the ff $\mathbb{N}^{V}$ can be easily computed on the element boundary. Moreover, also the integral of those functions can be evaluated. Indeed, from the definition of the internal degrees of freedom one gets that

[TABLE]

3.2.1 Consistent term

In the present section we describe the consistent part of the bilinear form $a^{E}_{h}(\cdot,\cdot)$, starting by the practical construction of the projection operator $\Pi$ through the definition of the associated matrix $\Pi^{\textrm{m}}$

[TABLE]

with $\Pi^{\textrm{m}}\in\mathbb{R}^{9\times(4m+2)}$. By inserting the above positions into Eq. (11) one again obtains (20), that now holds for all $\hat{\boldsymbol{\varepsilon}}$ in $\mathbb{R}^{9}$. Integrating by parts, identity (20) yields

[TABLE]

with the divergence matrix $\boldsymbol{\partial}={\mathbb{S}}^{\textrm{T}}$.

Eq. (3.2.1) can be written as

[TABLE]

where

[TABLE]

Therefore, the required operator $\Pi^{\textrm{m}}$ is given by the solution of a linear system:

[TABLE]

Remark 5

Note that matrix $\mathscr{B}$ is also computable. Indeed, the first term in the right hand side of (43) can be calculated: recall that the columns of $\mathbb{N}^{V}$, representing the ${{\bf V}_{h|E}}$ basis functions, are piecewise quadratic functions explicitly known on the boundary, see also Remark 6 below. Regarding the second term in the right hand side of (43), it can be immediately computed from the observation that $\boldsymbol{\partial}\mathbb{N}^{P}$ has constant entries

[TABLE]

and using (38).

Combining eqs. (39) and (44), from (12), one obtains:

[TABLE]

which defines the consistent symmetric stiffness matrix:

[TABLE]

Remark 6

The boundary term in (43) is split as a sum of edge integrals involving the product of a polynomial of degree $\leq 2$ and a polynomial of degree $\leq 1$. Therefore, a Gauss-Lobatto integration rule with three nodes on the edge (that is, the endpoints and the midpoint) is sufficient to compute exactly such integrals. This means that in practice one needs to compute the value of the matrix $\mathbb{N}^{V}$ only at the element vertexes and element edge midpoints. Due to our choice of the basis functions, and since such points are associated to the degrees of freedom of the space ${{\bf V}_{h|E}}$, the matrix $\mathbb{N}^{V}$ always takes either the value [math] or the value $1$ at the vertexes and the edge midpoints. This observation clearly simplifies the computation of matrix $\mathscr{B}$.

3.2.2 Stabilization term

Consider the following basis for the space $(\mathsf{P}_{2}(E))^{2}$:

[TABLE]

where we recall that we are using local scaled coordinates (cf. Eq. (13)). Similarly to the case $k=1$, the elements $\mathbb{p}\in(\mathsf{P}_{2}(E))^{2}$ can be written in terms of the above basis as vectors $\hat{\mathbb{p}}\in\mathbb{R}^{12}$ and, since $(\mathsf{P}_{2}(E))^{2}\subseteq{{\bf V}_{h|E}}$, also in terms of the basis of ${{\bf V}_{h|E}}$. We call ${\bf D}$ the matrix associated to such change of basis, that can be simply calculated by evaluating the functions in (47) through the degrees of freedom of ${{\bf V}_{h|E}}$. One has

[TABLE]

where $\tilde{\mathbb{p}}\in\mathbb{R}^{(4m+2)}$ represents $\mathbb{p}$ written in terms of the ${{\bf V}_{h|E}}$ basis. The matrix ${\bf D}$ is given by

[TABLE]

where the $i-th$ vertex of the polygon is supposed to have scaled coordinates $(\xi_{i}\;\eta_{i})$, $i=1,...,m$, the $i-th$ edge midpoint is assumed to have coordinates $(\zeta_{i}\;\omega_{i})$, $i=1,...,m$, and we denote with the short symbol $\int\!\!\!\!\!\!-\>p$ the integral average of any scalar polynomial $p$ on the polygon $E$.

Once the matrix ${\bf D}$ is built, the rest of the procedure follows exactly the same line of the case $k=1$, namely Eq. (30) for the construction of the stabilizing matrix, and Eq. (31) for the final stiffness matrix.

Remark 7

Analogously to Remark 4, it is easy to check that in this case the following “patch test” condition holds (in the case of constant elastic coefficients)

[TABLE]

3.2.3 Loading term

The loading term for the case $k=2$ is simpler than $k=1$, due to the presence of the moment degrees of freedom. Indeed, after the usual element-wise splitting (32), we approximate the local loading term by

[TABLE]

that corresponds to an approximation of the local load by its average. Following the notation of Section 3.1.3, this yields

[TABLE]

that can be immediately computed using (38).

3.3 The case of general order $k$

The degrees of freedom for the case of general order $k$ are given in Section 2.2. In particular, for any edge of a polygon, edge degrees of freedom are interpolatory displacements at points $\{\boldsymbol{x}^{e}_{i}\}$, $i=1,...,k-1$, chosen as the internal nodes of the Gauss-Lobatto integration rule with $k+1$ nodes. For instance, for $k=3$ the Gauss-Lobatto integration rule leads to $k+1=4$ nodes, i.e. $\{\boldsymbol{x}^{e}_{i}\}$, $i=1,2$ internal nodes, plus $2$ extrema. In passing, we note that, for $k=2$, this scheme provides the edge midpoint of the previous section. The definition of the $k(k-1)$ internal degrees of freedom is set by considering the following basis for $\mathsf{P}_{k-2}(E)$

[TABLE]

The monomials of basis (50) are ordered as the exponent couples $(\alpha\;\beta)$ of Pascal’s triangle and named accordingly $q_{1},q_{2},...,q_{r}$, with $r=k(k-1)/2$ the dimension of $\mathsf{P}_{k-2}(E)$. This ordered basis results:

[TABLE]

The moments of the displacement components, playing the role of internal degrees of freedom, are the $2r$ linear operators from ${{\bf V}_{h|E}}$ into $\mathbb{R}$ defined as:

[TABLE]

for all ${\mathbb{v}}_{h}=(u_{h}\;v_{h})^{T}\in{{\bf V}_{h|E}}$. The internal degrees of freedom are simply the scaled moments of the function ${\mathbb{v}}_{h}$ (for the first and second component) up to order $k-2$. For example, for $k=3$, the $2r=6$ internal degrees of freedom are:

[TABLE]

The dimensions of the virtual element space and of the polynomial space of approximated strains (cf. Eqs. (6), (10)) are, respectively:

[TABLE]

We order the basis for $\mathsf{P}_{k-1}(E)^{2\times 2}_{\textrm{sym}}$ coherently with the ordering of $q_{i}$’s:

[TABLE]

Representation (16) still applies, where ${\mathbb{v}}_{h}\in{{\bf V}_{h|E}}$, $\tilde{{\mathbb{v}}}_{h}\in\mathbb{R}^{n}$, $\boldsymbol{\varepsilon}^{P}\in\mathsf{P}_{k-1}(E)^{2\times 2}_{\textrm{sym}}$, $\hat{\boldsymbol{\varepsilon}}\in\mathbb{R}^{\ell}$. The matrix $\mathbb{N}^{P}\in\mathbb{R}^{3\times\ell}$ is immediately derived from the basis in (54)

[TABLE]

$\mathbb{N}^{V}$ is a $2\times n$ matrix operator representing the basis function of ${{\bf V}_{h|E}}$; odd columns have zero second entry, even columns have the reverse property. Each shape function takes the value $1$ for the corresponding degree of freedom and vanishes for all the others. The ordering of the degrees of freedom in vector $\tilde{{\mathbb{v}}}_{h}$ (and associated basis functions) is the same one for the case $k=2$, for $\tilde{u}_{h}$ and $\tilde{v}_{h}$ displacement components, respectively. Such indexing depends on the usual indexing of vertexes, edge nodes, and on the ordering of basis functions (51). It is reported in Table 1 for compactness.

Remark 8

The basis functions associated to the internal degrees of freedom clearly vanish on the boundary, since by definition these functions take zero value at all the boundary nodes. The internal basis functions are ordered in the same way as the polynomials $\{q_{1},q_{2},...,q_{r}\}$ introduced above. Namely, for $j=1,...,r$ and $i=km+1,...,km+r$ (that is the index range associated to the internal basis functions accordingly to the ordering in Table 1) we have for the first component

[TABLE]

with $\delta_{ij}$ representing the Kronecker symbol. The internal basis functions associated to the second component follow an analogous rule. For $j=1,...,r$ and $i=km+1,...,km+r$, it holds:

[TABLE]

Similarly to the previous cases, the functions of ${{\bf V}_{h|E}}$ on the boundary of $E$ are continuous and piecewise polynomials of degree $k$. Therefore the evaluation of $\mathbb{N}^{V}$ can be easily computed on the element boundary. Moreover, the integral on $E$ of the basis functions multiplied by any polynomial $q_{j}$ is computable immediately by using Eqs. (56)-(57). This allows to deduce that, for all $j=1,2,..,r$, the integral matrix

[TABLE]

has only two non-vanishing entries, with both value $|E|$, at position $(1,2km+2j-1)$ and $(2,2km+2j)$, respectively. Note that the first $2km$ columns of zeros are associated with the boundary degrees of freedom. For instance, in the case $k=3$ (recalling that $q_{1}=1,q_{2}=\xi,q_{3}=\eta$) one gets

[TABLE]

3.3.1 Consistent term

The practical construction of the projection operator $\Pi$, through the definition of the associated matrix $\Pi^{\textrm{m}}$

[TABLE]

follows the same steps shown for the $k=2$ case in Section 3.2. Indeed, by inserting the above definitions into Eq. (11), one again obtains (20), that now holds for all $\hat{\boldsymbol{\varepsilon}}$ in $\mathbb{R}^{\ell}$. Integrating by parts, identity (20) again gives (3.2.1), that in turn still can be written as in (41) and thus yields the fundamental Eq. (44). Clearly, now the dimensions involved are different:

[TABLE]

Finally, the consistent stiffness matrix $\mathbb{K}_{\textrm{c}}$ can be computed as in (46).

We close this section with important observations regarding the construction of matrix $\mathscr{B}$.

The first term on the right hand side of (61) can be calculated recalling that the columns of $\mathbb{N}^{V}$ (representing the basis functions of ${{\bf V}_{h|E}}$) are piecewise polynomial functions of degree $\leq k$ that are explicitly known on the boundary. Moreover, such calculation can be greatly simplified if one adopts as edge degrees of freedom the $(k-1)$ internal points associated to the Gauss-Lobatto integration rule, as already suggested. Indeed, the boundary term in (61) is split as a sum of edge integrals involving the product of a polynomial of degree $\leq k$ and a polynomial of degree $(k-1)$. Therefore, a Gauss-Lobatto integration rule with $(k+1)$ nodes on the edge (that is, the vertex nodes and the $(k-1)$ internal edge nodes) is sufficient to compute exactly such integrals. This means that in practice one needs to compute the value of the matrix $\mathbb{N}^{V}$ only at the element vertexes and element edge nodes $\{\boldsymbol{x}^{e}_{i}\}$, $i=1,...,(k-1)$. Due to our choice of the basis functions, and since such points are associated with the degrees of freedom of the space ${{\bf V}_{h|E}}$, the matrix $\mathbb{N}^{V}$ always takes either the value [math] or the value $1$ at the vertexes and the edge nodes. This observation clearly simplifies the computation of matrix $\mathscr{B}$.

Regarding the second term in the right hand side of (61), it can be computed by noting that the entries of $\boldsymbol{\partial}\mathbb{N}^{P}$ are polynomials of degree $\leq(k-2)$ and using (58). First of all, one needs to express the components of the polynomial matrix $\boldsymbol{\partial}\mathbb{N}^{P}$ in terms of the $q_{j}$ basis, introduced previously in the present section. Indeed, one can write

[TABLE]

where the matrices $M^{j}\in\mathbb{R}^{2\times\ell}$ have constant entries. Note that, given $k$, the matrices $M_{j}$ can be computed offline once and for all with symbolic calculus tools, since (62) is a polynomial identity and does not depend on the element $E$ (see also Remark 9 below). Having the above matrices, one can simply write

[TABLE]

and make use of (58).

Remark 9

Note that the matrix $\int_{E}q_{j}\mathbb{N}^{V}$ in (58) is a sparse matrix (only two components are different from zero) and therefore the calculations above can be done at minimal computational cost. Moreover, the matrix $\int_{E}q_{j}\mathbb{N}^{V}$ can be immediately written, for every element $E$, as a sparse matrix with only two unit entries, multiplied by the factor $|E|$. Therefore, the product matrices $(M^{j})^{T}\int_{E}q_{j}\mathbb{N}^{V}$ have the following properties: (1) only two columns are different from zero; (2) they can be computed and stored offline very easily, since the only factor depending on the element is the scalar factor $|E|$, which can be included later on, during each element stiffness matrix construction.

3.3.2 Stabilization term

The guidelines of the procedure are the same as in the cases $k=1$ and $k=2$. Therefore, since the construction of the stabilization term has been already explained with full details in Section 3.2.2, we here limit to a compact description of the involved terms.

We start by introducing a basis $\{{\mathsf{p}}_{j}\}$, $j=1,...,(k+1)(k+2)$, for the space $(\mathsf{P}_{k}(E))^{2}$. We set

[TABLE]

where the scaled monomials $q_{j}$ have been introduced at the beginning of Section 3.3 (cf. (50)). Note that for $k=2$ we obtain the same identical basis shown in (47).

Similarly to the previous cases, the elements $\mathbb{p}\in(\mathsf{P}_{k}(E))^{2}$ can be written in terms of the above basis as vectors $\hat{\mathbb{p}}\in\mathbb{R}^{(k+1)(k+2)}$ and, since $(\mathsf{P}_{k}(E))^{2}\subseteq{{\bf V}_{h|E}}$, also in terms of the basis of ${{\bf V}_{h|E}}$. As usual we call ${\bf D}$ the matrix associated to such change of basis, that can be simply computed by evaluating the polynomials $\{{\mathsf{p}}_{j}\}$, $j=1,...,(k+1)(k+2)$, through the degrees of freedom of ${{\bf V}_{h|E}}$. It holds:

[TABLE]

where $\tilde{\mathbb{p}}\in\mathbb{R}^{n}$ represents $\mathbb{p}$ written in terms of the ${{\bf V}_{h|E}}$ basis. The components of the matrix ${\bf D}$ are given by

[TABLE]

and can be calculated as follows.

Let $\{\nu_{i}\}$, $i=1,...,m$ be the set of the polygon vertexes, and $\{\boldsymbol{x}^{e}_{i}\}$, $i=(m+1),...,km$ the edge nodes, associated with the space ${{\bf V}_{h|E}}$. These are assumed to be ordered as described in Section 3.3, i.e. first the vertex nodes and then the edge nodes (given some ordering of the edges and the local ordering of the Gauss-Lobatto nodes along each edge). Thus, the first $2km$ rows of the matrix ${\bf D}$ are given by

[TABLE]

and all the remaining components of the first $2km$ rows are zero. It may be helpful for the reader to give a look at the explicit matrix representation given in (49) for the case $k=2$.

The remaining $2r=k(k-1)$ rows of ${\bf D}$ are given by

[TABLE]

and all the remaining components of the last $2r$ rows are zero.

Once the matrix ${\bf D}$ is built, the rest of the construction is identical to the previous cases, namely Eq. (30) is used for the construction of the stabilizing matrix, and Eq. (31) for the final stiffness matrix. Analogously to Remark 4, it is easy to check that in this case the following “patch test” condition holds (for constant elastic coefficients)

[TABLE]

3.3.3 Loading term

As usual, we start by the element-wise splitting (32), and we approximate the local loading term $\int_{E}{\mathbb{v}}^{T}\mathbb{b}$. The first step is to introduce $\mathbb{b}_{h}\in(\mathsf{P}_{k-2}(E))^{2}$ as the $L^{2}$ projection of the load $\mathbb{b}|_{E}$ onto the polynomials of degree $(k-2)$. The second step is simply to define

[TABLE]

We now describe such steps more in detail. Since $\mathbb{b}_{h}\in(\mathsf{P}_{k-2}(E))^{2}$, $\mathbb{b}_{h}$ can be expressed in terms of the basis (51):

[TABLE]

The collection of the above real numbers leads to the following vector $\widetilde{\mathbb{b}}_{h}\in\mathbb{R}^{2r}$:

[TABLE]

By definition of $L^{2}$ projection, and using (65), it is easy to check that the coefficient vectors in (65) must satisfy

[TABLE]

The above linear system can be explicitly written for the coefficient vector (66), as

[TABLE]

where

[TABLE]

Finally, once the $\widetilde{\mathbb{b}}_{h}$ vectors is computed as shown above, using (64) and (65), we have

[TABLE]

that is immediately computed recalling (58). The same observations of Remark 9 still apply.

4 Numerical tests

The proposed VEM methods are implemented into in-house Matlab matlab:2015 codes for validation and analysis. In this section we assess accuracy and robustness through a number of representative boundary value problems on a unit square domain plus a classical 2D problem taken from the literature. They are organized as follows:

• $\bullet$

Tests $1a$, $1b$: 2D plane strain elasticity linear patch tests on the unit square $\Omega=[0,1]^{2}$;

• $\bullet$

Tests $2a$, $2b$: 2D plane strain convergence tests with known analytical solution on the unit square $\Omega=[0,1]^{2}$;

• $\bullet$

Test $3$: assessment of the influence of the stabilization parameter;

• $\bullet$

Test $4$: Cook’s Membrane.

For the sake of simplicity, the geometric and physical quantities are set omitting the system of measurement; it is nonetheless assumed that the relevant units are taken consistently. In all the numerical tests we use the choice $\tau=1/2$ for the stabilization parameter, see (30). In section 4.3 we show that taking other choices has limited influence on the results.

4.1 Test $1a$ and $1b$

We consider two boundary value problems on the unit square domain $\Omega=[0,1]^{2}$ which induce constant stress states examined in Bathe1996 , under plain strain assumption. The material parameters are introduced setting the Young modulus $E=7000$ and the Poisson ratio $\nu=0.3$. Geometry, boundary conditions and loading cases for Test $1a$ [resp. Test $1b$] are reported in Fig. 4(a) [resp. Fig. 4(b)]. The former case represents a constant tensile stress state induced by the normal traction $q=2000$ on the right edge of the unit square, the latter represents a constant shear stress state induced by the tangential traction $t=400$ assigned on the whole boundary.

Linear VEMs are sufficient for the problem under investigation. Test $1a$ is analyzed with a very coarse mesh made of $4$ equal concave quadrilaterals and one square; the latter is centered in the center of $\Omega$ and has a side length $1/3$, and is rotated counter-clockwise of $\pi/3$ as can be seen in Fig. 5(a), where the meshes for the undeformed (dotted line) and deformed (continuous line) configurations are reported. Test $1b$ is analyzed with a very coarse mesh made of $4$ equal concave pentagons and one star-shaped octagon; the latter is centered in the center of $\Omega$ and has is derived by the square element of Test $1a$ by punching each midside point through the inside of a quarter of its side length, plus a counter-clockwise of $\pi/3$ as can be seen in Fig. 5(b), where the meshes for the undeformed (dotted line) and deformed (continuous line) configurations are reported.

Qualitative color plots111For brevity, contour plots are interpolated through vertex values using command fill matlab:2015 . for the displacement components stemming from the two test solutions are shown in Fig. 6: as expected a linear pattern can be appreciated in both cases, being horizontal displacement for Test $1b$ obviously zero. The order of magnitude of the maximum euclidean norm of the difference between the analytical and the numerical stresses over a $100\times 100$ uniform grid, and over the whole domain, is $10^{-16}$ in both cases, ascertaining that the method is capable of exactly representing a simple stress state.

4.2 Test $2a$ and $2b$

We consider two boundary value problems on the unit square domain $\Omega=[0,1]^{2}$, with known analytical solution, discussed in BeiraoLovaMora . In this case, the material parameters are assigned in terms of Lamé constants $\lambda=1$, $\mu=1$, and plain strain regime is still invoked. The tests are defined by choosing a required solution and deriving the corresponding load $\mathbb{b}$, as synthetically indicated in the following:

• $\bullet$

Test $2a$

[TABLE]

• $\bullet$

Test $2b$

[TABLE]

As it can be immediately observed, Test $2a$ is a problem with Dirichlet non-homogeneous boundary conditions, zero loading and a polynomial solution; whereas Test $2b$ has homogeneous Dirichlet boundary conditions, trigonometric distributed loads with a trigonometric solution.

For the purpose of comparing the accuracy level and the convergence rate of the proposed VEM to standard FEM, meshes with only quadrilaterals will be considered for Test $2a$. Linear and quadratic VEMs are compared to Lagrangian linear and quadratic quadrilateral finite elements, indicated in the following as $Q4$ and $Q9$, respectively Zienckiewicz_Taylor2000 . Test $2b$ is instead analyzed with general polygonal meshes only testing VEM solutions.

In the following, such accuracy and convergence rate assessment is carried out using the following error norms:

• $\bullet$

Energy-type error norm:

[TABLE]

where $\boldsymbol{\boldsymbol{\varepsilon}}_{\textrm{ex}}$ is the exact strain, and $\boldsymbol{\boldsymbol{\varepsilon}}({\mathbb{v}}_{h})$ is the numerical strain according to either a FEM or a VEM solution. In particular, for the VEM solution, such a numerical strain is computed by means of (10) and (11) while for the FEM solution it is obtained by standard derivation of the displacements.

• $\bullet$

Discrete $H^{1}$-type error norm:

[TABLE]

where $h_{f}$ denotes the length of any edge $f$ of the reticulation ${\cal T}_{h}$, and ${\mathbb{t}}_{f}$ denotes the unit tangent vector to the edge $f$ chosen once and for all.

Sample meshes for Test $2a$ are represented in Fig. 7, comprising evenly distributed squares (a), rhombic and concave quadrilaterals (b), and rectangular trapezoids almost collapsing into triangles (c). Correspondingly, Figure 8(a)-(f) show the $h-$convergence plots in terms of error $D_{1}$ and $D_{2}$ for the $Q4$, $Q9$ FEM, and $k=1$, $k=2$ VEM solutions. It is observed that, the tested methods present the expected convergence rates, namely linear and quadratic, respectively. The comparison in terms of error $D_{1}$ indicates that FEM methods slightly outperform the VEM of same order; while the opposite patterns is encountered in terms of error $D_{2}$. This seems to indicate that the VEM solution on the mesh skeleton is slightly better than the FEM solution in the presence of distorted meshes (see $D_{2}$ errors). On the other hand the VEM strain, since it is computer through the “filter” of the projection operator instead of being directly derived from the local displacements, suffers from a slight loss of accuracy (see $D_{1}$ errors). Possibly making use of more sofisticated post-processing procedures to compute the VEM strains could yield results that are more accurate than FEM. Nevertheless, it is important to remark that all the errors shown in Figure 8 are comparable; this result fits exactly into the scope of the present test $2a$, that is to show than even for conforming quadrilaterals (that is, playing on a ground where FEM are very strong) the VEM is able to obtain competitive results.

Sample meshes for Test $2b$ are represented in Fig. 9, comprising non-uniform mesh of quadrilaterals (a), non-uniform mesh of triangles (b), uniform mesh of convex hexagons (c), and a centroid-based Voronoi tessellation (d). Correspondingly, Figures 10(a)-(d) show the $h-$convergence plots in terms of error $D_{1}$ for the $k=1$, and $k=2$ VEM solutions. It is again observed that VEM show the expected convergence rates, which, as for Test $2a$, are maintained regardless of mesh non-uniformity and distortion. Similar considerations may be drawn in relation to the $h-$convergence pattern for error $D_{2}$, as can be appreciated in Fig. 11.

4.3 Test $3$

The present section is devoted to the analysis of the influence on solution accuracy of the stabilization parameter, appearing in (30) pre-multiplying the stabilization stiffness, taken by default as $\tau=1/2$. Reference is made to Test $2a$, solved on the four different meshes reported in Fig. 12, comprising: non-uniform mesh of triangles (a), distorted convex quadrilaterals (b), a centroid-based Voronoi tessellation, and a random-based Voronoi tessellation. Correspondingly, we plot error $D_{1}$ (cf. (77)) in Fig. 13 against a parameter $\alpha_{0}$ ranging from $10^{-2}$ to $10^{2}$, with the corresponding stabilization parameter taken as $\tau=\alpha_{0}/2$, cf. equation (30). Note that for $\alpha_{0}=1$ the standard choice $\tau=1/2$ adopted in the tests of this paper is recovered, while for different values of $\alpha_{0}$ stronger or weaker stabilizations are introduced. We consider both $k=1$ and $k=2$ VEMs. It can be observed that each examined case seems to be quite insensitive to the variation of $\tau$: the curves are quite flat and the overall variation of error with respect to the scaling is limited to within less than one order of magnitude. The proposed method thus seems to be quite robust with respect to the stabilization parameter.

4.4 Test $4$

The present section deals with the classical Cook’s membrane 2D problem Zienckiewicz_Taylor2000 . Geometry of the domain $\Omega$ is presented in Fig. 14 with length data $H_{1}=44$, $H_{2}=16$, $L=48$. Material parameters are Young modulus $E=70$, and Poisson ration $\nu=1/3$. The loading is given by a constant tangential traction $q=6.25$ on the right edge of the domain. The problem is solved using two types of meshes: an evenly distributed quad mesh (see fig. 14(a)), and centroid based [resp. random-based] Voronoi tessellations (see fig. 14(b)). The former type is used for linear and quadratic VEMs in conjunction with $Q4$ and $Q9$ FEM for comparison purposes; whereas the latter types are adopted in conjunction to VEMs only. Convergence results are reported in terms of mesh refinement monitoring $v_{A}$, the vertical displacement of point A (see Fig. 14), as can be appreciated in Fig. 16(a)-(b), where a reference solution is indicated with a dotted black line corresponding to an overkilling accurate solution obtained with the hybrid-mixed CPE4I element ABAQUS:2011 . It is observed that quadrilateral VEMs have a slight edge in terms of accuracy with respect to Lagrangian finite elements of the corresponding order, while the polygonal VEMs from the centroid-based Voronoi tessellation slightly outperform the random-based VEMs in terms of accuracy.

In Fig. 17 the deformed configuration corresponding to quadratic quadrilaterals (a) and centroidal Voronoi polygons (b) is reported, showing a perfect agreement between the two solutions and hence no sensitivity to mesh type, distortion and non-uniformity. Finally, in Fig. 18, the qualitative color plots of the displacement components stemming from the latter solution are displayed showing the expected smoothness of the displacement field over the mesh.

5 Conclusion

A Virtual Element method for plane elasticity problems has been presented. The newly developed methodology, together with implementation details, has been illustrated for the case of quadrilateral elements with linear approximation of the displacement field, and for general polygonal elements with quadratic and higher order interpolation. The stiffness matrix and loading vector have been explicitly derived proving the straightforward implementability of the method into into a computer code for structural analysis. An extensive campaign of numerical applications illustrate accuracy and convergence patterns of the method, and its performances in comparison with a classical finite element approach characterized by the same order of approximation of the unknown fields. An accurate investigation on the convergence for structured and non structured discretizations has been performed, showing the superior robustness of the VEM method with respect to a FEM approach. It can be stated that VEM is almost insensitive to the mesh distortion. The assessed reliability and the good performances of the VEM formulation indicate that the virtual element approach can rapidly become a numerical strategy largely diffused and used for structural computations.