A Dual Hybrid Virtual Element Method for Plane Elasticity Problems

TL;DR

This paper introduces a dual hybrid Virtual Element Method for plane elasticity problems, demonstrating stability, convergence, and effectiveness through numerical tests on benchmarks and microelectromechanical systems.

Contribution

It presents a novel dual hybrid Virtual Element scheme with proven stability and convergence, offering an effective alternative to traditional displacement-based methods.

Findings

First order convergence demonstrated

Numerical tests confirm method's accuracy

Validates dual hybrid scheme as an alternative

Abstract

A dual hybrid Virtual Element scheme for plane linear elastic problems is presented and analysed. In particular, stability and convergence results have been established. The method, which is first order convergent, has been numerically tested on two benchmarks with closed form solution, and on a typical microelectromechanical system. The numerical outcomes have proved that the dual hybrid scheme represents a valid alternative to the more classical low-order displacement-based Virtual Element Method.

Equations

$$\left\{{\begin{aligned} &\mbox{Find }({\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt},\mathbf{u})~{}\mbox{such that}\\ &-\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}=\mathbf{f}\qquad\mbox{in $\Omega$}\\ &{\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}=\mathbb{C}{{\kern 0.20004pt\hbox{$\varepsilon$}\kern-4.66318pt\kern-0.20004pt\hbox{$\varepsilon$}\kern-4.66318pt\raise 0.29999pt\hbox{$\varepsilon$}\kern 0.20004pt}}(\mathbf{u})\qquad\mbox{in $\Omega$}\\ &\mathbf{u}_{|\partial\Omega}={\bf 0},\end{aligned}}\right.$$

$$h_E:=\mathrm{diameter}(E),h:=\underset{E\in {\mathcal{T}}_h}{\sup }{h}_E.$$

$$\mathcal{E}(\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v})=-\frac{1}{2}\sum _{E\in {\mathcal{T}}_h}{a}_E(\tau \tau \text{raise}0.29999pt\tau ,\tau \tau \text{raise}0.29999pt\tau )+\sum _E{\int }_{\partial E}\tau \tau \text{raise}0.29999pt\tau \mathbf{n}\cdot \mathbf{v}\tau \tau \text{raise}0.29999pt\tau \in {\Sigma }^f,\mathbf{v}\in U^0.$$

$$\{\begin{array}{rl}& {\Sigma }^f=\prod _{E\in {\mathcal{T}}_h}{\Sigma }^f(E),\\ & U^0=H_0^1(\Omega {)}^2,\end{array}$$

$$\Sigma^{f}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in H(\mathop{\bf div}\nolimits;E)_{s}\,:\,\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}+\mathbf{f}={\bf 0}\Big{\}}$$

$$H(\mathop{\bf div}\nolimits;E)_{s}=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in L^{2}(E)^{2\times 2}\,:\,\mbox{${\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}$ is symmetric,}\,\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in L^{2}(E)^{2}\Big{\}}.$$

$$||\tau \tau \text{raise}0.29999pt\tau |{|}_{\Sigma }:={\left(\sum _{E\in {\mathcal{T}}_h}(||\tau \tau \text{raise}0.29999pt\tau |{|}_{0,E}^2+||\mathrm{div}\tau \tau \text{raise}0.29999pt\tau |{|}_{0,E}^2)\right)}^{1/2}.$$

$${\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}={\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}^{0}+\widehat{{\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}}_{f}\qquad\mbox{with ${\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}^{0}\in\Sigma^{0}$}.$$

$$\{\begin{array}{rl}& \text{mbox}Find({\sigma \sigma \text{raise}0.29999pt\sigma }^0,\mathbf{u})\in {\Sigma }^0\times U^0\text{mbox}suchthat\\ & a({\sigma \sigma \text{raise}0.29999pt\sigma }^0,{\tau \tau \text{raise}0.29999pt\tau }^0)+b({\tau \tau \text{raise}0.29999pt\tau }^0,\mathbf{u})=F({\tau \tau \text{raise}0.29999pt\tau }^0)\forall {\tau \tau \text{raise}0.29999pt\tau }^0\in {\Sigma }^0\\ & b({\sigma \sigma \text{raise}0.29999pt\sigma }^0,\mathbf{v})=G(\mathbf{v})\forall \mathbf{v}\in U^0,\end{array}$$

$$\{\begin{array}{rl}& a({\sigma \sigma \text{raise}0.29999pt\sigma }^0,{\tau \tau \text{raise}0.29999pt\tau }^0)=\sum _{E\in {\mathcal{T}}_h}{a}_E({\sigma \sigma \text{raise}0.29999pt\sigma }^0,{\tau \tau \text{raise}0.29999pt\tau }^0),\\ & b({\tau \tau \text{raise}0.29999pt\tau }^0,\mathbf{v})=-\sum _{E\in {\mathcal{T}}_h}{\int }_{\partial E}{\tau \tau \text{raise}0.29999pt\tau }^0\mathbf{n}\cdot \mathbf{v},\\ & F({\tau \tau \text{raise}0.29999pt\tau }^0)=-\sum _{E\in {\mathcal{T}}_h}{a}_E({\widehat{\sigma \sigma \text{raise}0.29999pt\sigma }}_f,{\tau \tau \text{raise}0.29999pt\tau }^0)\\ & G(\mathbf{v})=\sum _{E\in {\mathcal{T}}_h}{\int }_{\partial E}{\widehat{\sigma \sigma \text{raise}0.29999pt\sigma }}_f\mathbf{n}\cdot \mathbf{v}.\end{array}$$

$$H=\Big{\{}\mathbf{v}\in U^{0}\,:\,b({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt},\mathbf{v})=0\quad\forall\,{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in\Sigma^{0}\Big{\}}=\Big{\{}\mathbf{v}\in U^{0}\,:\,\mathbf{v}_{|{\mathcal{E}}_{h}}={\bf 0}\Big{\}}.$$

$$\underset{\tau \tau \text{raise}0.20999pt\tau \in {\Sigma }^0}{\sup }\frac{b(\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v})}{||\tau \tau \text{raise}0.29999pt\tau |{|}_{\Sigma }}\gtrsim ||\mathbf{v}|{|}_{U^0/H}\forall \mathbf{v}\in U^0.$$

$$K=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in\Sigma^{0}\,:\,b({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt},\mathbf{v})=0\quad\forall\,\mathbf{v}\in U^{0}\Big{\}}=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}\in H(\mathop{\bf div}\nolimits;\Omega)\,:\,\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}={\bf 0}\Big{\}}.$$

$$a(\tau \tau \text{raise}0.29999pt\tau ,\tau \tau \text{raise}0.29999pt\tau )\gtrsim ||\tau \tau \text{raise}0.29999pt\tau |{|}_{\Sigma }^2\forall \tau \tau \text{raise}0.29999pt\tau \in K.$$

$$\{\begin{array}{rl}& \text{mbox}Find(\sigma \sigma \text{raise}0.29999pt\sigma ,\mathbf{u})\in {\Sigma }^f\times U^0\text{mbox}suchthat\\ & a(\sigma \sigma \text{raise}0.29999pt\sigma ,{\tau \tau \text{raise}0.29999pt\tau }^0)+b({\tau \tau \text{raise}0.29999pt\tau }^0,\mathbf{u})=0\forall {\tau \tau \text{raise}0.29999pt\tau }^0\in {\Sigma }^0\\ & b(\sigma \sigma \text{raise}0.29999pt\sigma ,\mathbf{v})=0\forall \mathbf{v}\in U^0,\end{array}$$

$$RM(E)=\left\{\mathbf{r}(\mathbf{x})=\mathbf{a}+b(\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }\mathbf{a}\in {\mathbb{R}}^2,b\in \mathbb{R}\right\}.$$

$$R(e)=\left\{\mathbf{t}(s)=\mathbf{c}+ds\mathbf{n}\mathbf{c}\in {\mathbb{R}}^2,d\in \mathbb{R},s\in [-1/2,1/2]\right\}.$$

$$\displaystyle\Sigma_{h}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in$$

$$\displaystyle({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}\in R(e)\quad\forall e\in\partial E;\quad\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in RM(E)\Big{\}}.$$

$$\displaystyle\Sigma_{h}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in$$

$$\displaystyle({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}\in R(e)\quad\forall e\in\partial E;\quad\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in RM(E)\Big{\}}.$$

$${\int }_E\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h\cdot \mathbf{r}={\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi \cdot \mathbf{r}\forall \mathbf{r}\in RM(E)$$

$$\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h={\alpha \alpha \text{raise}0.29999pt\alpha }_E+{\beta }_E(\mathbf{x}-{\mathbf{x}}_C{)}^{\perp },$$

$$\{\begin{array}{rl}& {\alpha \alpha \text{raise}0.29999pt\alpha }_E=\frac{1}{|E|}{\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi =\frac{1}{|E|}\sum _{e\in \partial E}{\int }_e{\mathbf{c}}_e\\ & {\beta }_E=\frac{1}{{\int }_E|\mathbf{x}-{\mathbf{x}}_C{|}^2}{\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi \cdot (\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }=\frac{1}{{\int }_E|\mathbf{x}-{\mathbf{x}}_C{|}^2}\sum _{e\in \partial E}{\int }_e({\mathbf{c}}_e+d_{e}s\mathbf{n})\cdot (\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }.\end{array}$$

$$\Sigma^{f}_{h}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in\Sigma_{h}(E)\,:\,\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}+\mathbf{f}={\bf 0}\Big{\}}.$$

$$U_{h}(E)=\Big{\{}\mathbf{v}_{h}\in H^{1}(E)^{2}\ :\ \mathbf{v}_{h|e}\in{\mathcal{P}}_{1}(e)^{2}\quad\forall e\in\partial E\Big{\}}.$$

$${\Sigma }_h^f(E)$$

$$\displaystyle d_{n_{E}}=-\frac{1}{\int_{e_{n_{E}}}s\,\mathbf{n}\cdot(\mathbf{x}-\mathbf{x}_{C})^{\perp}}\Big{(}\sum_{i=1}^{n_{E}}\int_{e_{i}}\mathbf{c}_{i}\cdot(\mathbf{x}-\mathbf{x}_{C})^{\perp}$$

$$\displaystyle\quad+\sum_{i=1}^{n_{E}-1}\int_{e_{i}}d_{i}s\,\mathbf{n}\cdot(\mathbf{x}-\mathbf{x}_{C})^{\perp}\Big{)}\Big{\}},$$

$${\int }_{\partial E}{\tau \tau \text{raise}0.29999pt\tau }_h\mathbf{n}\cdot {\mathbf{v}}_h$$

$$a_E({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h)={\int }_E\mathbb{D}{\sigma \sigma \text{raise}0.29999pt\sigma }_h:{\tau \tau \text{raise}0.29999pt\tau }_h$$

$$\{\begin{array}{rl}& {\Pi }_E:{\Sigma }_h(E)\to {\mathcal{P}}_{\ast }(E{)}_s^{2\times 2}\\ & {\tau \tau \text{raise}0.29999pt\tau }_h\mapsto {\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h\\ & a_E({\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h,\pi \pi \text{raise}0.29999pt\pi )=a_E({\tau \tau \text{raise}0.29999pt\tau }_h,\pi \pi \text{raise}0.29999pt\pi )\forall \pi \pi \text{raise}0.29999pt\pi \in {\mathcal{P}}_{\ast }(E{)}_s^{2\times 2},\end{array}$$

Full text

1 Introduction

The Virtual Element Method (VEM) is a recent methodology to approximate partial differential equation problems. Introduced in [10], it is a Galerkin method which can be considered as an evolution of the Finite Element Method (FEM). In contrast to FEM, VEM is able to naturally manage several mesh complexities, such as polytopal shapes or hanging nodes. In addition to this flexibility, it has been realized that VEM is also able to efficiently deal with other non-trivial situations, for instance problems with internal constraints (incompressibility for solids and fluids is an example), or problems with high-continuity requirements (the fourth-order Kirchhoff plate is an example). The price to pay is that the shape functions are not explicitly known and thus they are called virtual. However, the available information on them is sufficient to form the stiffness matrix and the right-hand side of the discretized problem. For the analysis of the VEM technique for the basic second-order elliptic problems, we refer to [10, 8, 14].

Focusing on the linear elasticity problem, the VEM philosophy has been already conjugated in several ways. In the easiest framework, the displacement-based variational formulation, VEMs have been considered in [7, 3, 22], and a procedure to recover an accurate stress field has been proposed in [6]. For nearly-incompressible materials, the VEM schemes proposed in [11] for the Stokes problem has been applied in connection with the displacement/pressure formulation, see for example [23, 30]. Furhermore, the Hellinger-Reissner variational principle has been recently employed to develop Virtual Element Methods with approximated stresses exhibiting, a priori, symmetry and inter-element traction continuity features, see [4, 5]. Among other works regarding linear elasticity problems, also concerning incompressible materials or plate structures, as well as non-conforming schemes or different variational frameworks, we here mention [1, 16, 33, 17, 9, 25, 34, 18, 31, 32].

In the present contribution, we explore the possibility to develop VEM schemes in the framework of the so-called dual hybrid formulation, see [12] for instance. Thus, we are concerned with a variational setting where the unknowns are both the stress and the displacement fields. Once the computational domain is partitioned into (polytopal) elements, the stress field is assumed to a priori satisfy the equlibrium equation locally on each element. It is then required to locally maximize the complementary energy; the displacement field enters into play only on the interelement boundaries and plays the rôle of the Lagrange multiplier for the traction continuity constraint. For Finite Elements, this variational approach has been used in [28, 27], for example. In our VEM scheme, the local stress space is borrowed from the one introduced in [4], while the standard low-order nodal Virtual Element space is essentially used for the displacement field. As usual for VEMs, a local polynomial projection is introduced as a basic ingredient to form the stiffness matrix of the method. Specifically, here, given a virtual stress, the VEM projection returns a suitable computable polynomial stress. In this paper we present and numerically investigate two different projections, the second of which is a low-cost improvement of the first one (already used in [4]). We remark that hybrid formulations for the linear elastic problem are not only interesting per se, but also they may be used as a building block for other more complex situations, for instance plate problems (see [21]).

A brief outline of the paper is as follows. In Section 2 we present the (2D) elastic problem together with its dual hybrid variational formulation. Section 3 describes the Virtual Element approximation we propose. In Section 4 we develop the stability and convergence analysis, by using a suitable mesh-dependent norm for the stress field. The mumerical results, which confirm the theoretical predictions, are detailed in Section 5, while some concluding remarks are drawn in Section 6.

Throughout the paper, given two quantities $a$ and $b$, we use the notation $a\lesssim b$ to mean: there exists a constant $C$, independent of the mesh-size, such that $a\leq C\,b$. If $a\lesssim b$ and $b\lesssim a$, we will write $a\approx b$. Moreover, we use standard notations for Sobolev spaces, norms and semi-norms (cf. [24], for example). Finally, given a subset $\omega\subseteq\mathbb{R}^{2}$, we will denote with ${\mathcal{P}}_{k}(\omega)$ ($k\geq 0$) the space of polynomials of degree up to $k$ and defined on $\omega$.

2 The elasticity problem in the dual hybrid form

In this section we briefly present the linear elasticity problem in dual hybrid form. More details about the dual hybrid formulations of second-order problems can be found in [12]. We start by considering the strong form of the problem we are interested in:

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is a polygonal domain, \sigma$$\sigma$$\sigma and $\mathbf{u}$ are the unknown stress and displacement fields, respectively. Moreover, $\mathbf{f}$ represents the body force density, $\mathbb{C}$ is the elasticity tensor and {{\kern 0.20004pt\hbox{\varepsilon}\kern-4.66318pt\kern-0.20004pt\hbox{\varepsilon}\kern-4.66318pt\raise 0.29999pt\hbox{\varepsilon}\kern 0.20004pt}}(\cdot) is the usual symmetric gradient operator. We consider only the clamped case along the whole boundary, but other boundary conditions can be treated in standard ways (see [12], for instance).

Let now $\{\mathcal{T}_{h}\}$ be a sequence of decompositions of $\Omega$ into general polygonal elements $E$ with

[TABLE]

Given ${\mathcal{T}}_{h}$, let us denote with ${\mathcal{E}}_{h}=\bigcup_{E\in{\mathcal{T}}_{h}}\partial E$ the skeleton of ${\mathcal{T}}_{h}$. The dual hybrid formulation of Problem (1) is a variation approach for which:

• •

the stress field \sigma$$\sigma$$\sigma a-priori satisfies the equilibrium equation in each element $E\in{\mathcal{T}}_{h}$;

• •

the displacement field $\mathbf{u}$ enters into play essentially only on the skeleton ${\mathcal{E}}_{h}$, where it acts as a Lagrange multiplier for the continuity of the tractions {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}\mathbf{n}.

More precisely, defining with $(\cdot,\cdot)_{E}$ the scalar product in $L^{2}(E)$, and a_{E}({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}):=(\mathbb{D}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt})_{E}, with $\mathbb{D}$ the compliance tensor (i.e. the inverse of the tensor $\mathbb{C}$), the dual hybrid formulation stems from considering the critical points of the following functional:

[TABLE]

Above, the spaces $\Sigma^{f}$ and $U^{0}$ are defined by:

[TABLE]

where

[TABLE]

and

[TABLE]

In $\Sigma^{f}$ we introduce the obvious norm:

[TABLE]

Selecting a particular locally self-equilibrated symmetric stress solution \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}\in\Sigma^{f}, i.e. such that (\mathop{\bf div}\nolimits\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}+\mathbf{f})_{|E}={\bf 0} for every $E\in{\mathcal{T}}_{h}$, the stress solution \sigma$$\sigma$$\sigma can be decomposed as

[TABLE]

Consequently, stationarity of functional (2) leads to the variational problem for the unknowns {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0} and $\mathbf{u}$:

[TABLE]

where $\Sigma^{0}$ is defined according with (3) and (4), choosing $\mathbf{f}={\bf 0}$. Above, for {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}^{0}\in\Sigma^{0} and $\mathbf{v}\in U^{0}$, we have set:

[TABLE]

Once {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0} has been found, the stress solution \sigma$$\sigma$$\sigma is simply recovered by using (6). We remark that $\mathbf{u}$, part of the solution to Problem (7), is not unique, but it is defined up to an element of the subspace

[TABLE]

Moreover, the following inf-sup condition holds (see [12]).

[TABLE]

In the sequel, we will also consider the following important subspace of $\Sigma^{0}$:

[TABLE]

We remark that, on $K$, we have ||{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}||_{\Sigma}=||{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}||_{0,\Omega} (cf. (5)), and it holds:

[TABLE]

Obviously, Problem (7) is equivalent to:

[TABLE]

From (10) and (12), the general theory of mixed methods (see [12], for instance) gives that Problem (7) has a unique solution ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0},\mathbf{u})\in\Sigma^{0}\times U^{0}/H (with an abuse of notation, we here use $\mathbf{u}$ to denote the equivalence class of the function $\mathbf{u}$ in $U^{0}/H$). Therefore, we infer that also Problem (13) has a unique solution ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u})\in\Sigma^{0}\times U^{0}/H. The quotient space $U^{0}/H$ essentially means that in Problems (7) and (13) the displacement field $\mathbf{u}$ enters only through its trace $\mathbf{u}_{|{\mathcal{E}}_{h}}$ on the skeleton ${\mathcal{E}}_{h}$, see (9). The following result can be deduced using the discussion in [29], and justifies the variation framework (7) (or (13)) for the elasticity Problem (1).

Proposition 2.1**.**

Let ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u}) be a sufficiently smooth solution to Problem (1). Then ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u})\in\Sigma^{0}\times U^{0}/H is the solution of Problem (13). ∎

3 The Virtual Element Method

We now outline the Virtual Element Method we propose for the discretization of Problem (7). To develop the theoretical analysis of the scheme, we suppose that for all $h$, each element $E$ in $\mathcal{T}_{h}$ fulfils the following assumptions:

• •

$\mathbf{(A1)}$ $E$ is star-shaped with respect to a ball of radius $\geq\,\gamma\,h_{E}$,

• •

$\mathbf{(A2)}$ the distance between any two vertexes of $E$ is $\geq c\,h_{E}$,

where $\gamma$ and $c$ are positive constants. We remark that the hypotheses above, though not too restrictive in many practical cases, can be further relaxed, as noted in  [10].

We suppose that the compliance tensor $\mathbb{D}$ is piecewise constant with respect to the underlying mesh ${\mathcal{T}}_{h}$. The analysis for a general (sufficiently smooth) tensor $\mathbb{D}$ follows from this case using the same arguments of [4]. In addition, we assume that the load term $\mathbf{f}$ is piecewise constant with respect to the underlying mesh ${\mathcal{T}}_{h}$. If $\mathbf{f}$ is indeed smooth, this modification introduces an $O(h)$ perturbation on the solution which does not spoil the convergence rate of our scheme.

3.1 The local spaces

Given a polygon $E\in{\mathcal{T}}_{h}$ with $n_{E}$ edges, we first introduce the space of local infinitesimal rigid body motions:

[TABLE]

Here above, given $\mathbf{c}=(c_{1},c_{2})^{T}\in\mathbb{R}^{2}$, $\mathbf{c}^{\perp}$ is the clock-wise rotated vector $\mathbf{c}^{\perp}=(c_{2},-c_{1})^{T}$, and $\mathbf{x}_{C}$ is the baricenter of $E$. For each edge $e$ of $\partial E$, we introduce the space

[TABLE]

Here above, $s$ is a local linear coordinate on $e$, such that $s=0$ corresponds to the edge midpoint. Furthermore, $\mathbf{n}$ is the outward unit normal to the edge $e$. Hence, $R(e)$ consists of vectorial functions which are constant in the edge tangential direction, while they are linear along the edge normal direction. Then, we set:

[TABLE]

Remark 1**.**

Alternatively, the space (16) can be defined as follows.

[TABLE]

Here above, the equation \mathop{\rm curl}\nolimits\mathop{\bf curl}\nolimits(\mathbb{D}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h})=0 is to be intended in the distribution sense.

We remark that, once ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}=\mathbf{c}_{e}+d_{e}s\,\mathbf{n} is given for all $e\in\partial E$, cf. (15), the quantity \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in RM(E) is determined. Indeed, denoting with {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}:\partial E\to\mathbb{R}^{2} the function such that {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}_{|e}:=\mathbf{c}_{e}+d_{e}s\,\mathbf{n}, the integration by parts formula

[TABLE]

allows to compute \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h} using the $\mathbf{c}_{e}$’s and the $d_{e}$’s. More precisely, setting (cf (14))

[TABLE]

from (18) we infer

[TABLE]

We now define the affine space for the local approximation of the stress field:

[TABLE]

The space $\Sigma^{0}_{h}(E)$ is defined according with (21), choosing $\mathbf{f}={\bf 0}$.

The local approximation space for the displacement field is defined as follows:

[TABLE]

We notice that $U_{h}(E)$ is an infinite dimensional space. However, this will not lead to any computational trouble, since only $\mathbf{v}_{h|\partial E}$, the trace of $\mathbf{v}_{h}\in U_{h}(E)$ on $\partial E$, enters into play in the discrete formulation. Therefore, we may think that the degrees of freedom for $U_{h}(E)$ are linear functionals which uniquely determine $\mathbf{v}_{h}\in U_{h}(E)$ on $\partial E$. For instance, one may take the point values of $\mathbf{v}_{h}$ at the vertices of $E$.

3.1.1 Computation of $\Sigma^{f}_{h}(E)$

We now show how the space $\Sigma^{f}_{h}(E)$ can be described, i.e. how a suitable set of degrees of freedom can be selected. We first number the edges of $\partial E$ as $e_{1},e_{2}\ldots,e_{n_{E}}$, once and for all. We recall that we have supposed $\mathbf{f}_{|E}\in{\mathcal{P}}_{0}(E)$. The following result holds.

Lemma 3.1**.**

Let $\mathbf{f}_{|E}\in{\mathcal{P}}_{0}(E)$. Then, the space $\Sigma^{f}_{h}(E)$ is characterised by:

[TABLE]

where we have set $\mathbf{c}_{i}=\mathbf{c}_{e_{i}}$ and $d_{i}=d_{e_{i}}$.

Proof.

We notice, see [4], that $\int_{e_{n_{E}}}s\,\mathbf{n}\cdot(\mathbf{x}-\mathbf{x}_{C})^{\perp}\neq 0$. Hence, the right-hand side of (23) is well-defined. Now, the proof easily follows by a direct computation from (19) and (20). ∎

Remark 2**.**

Obviously, the space $\Sigma^{0}_{h}(E)$ is defined by (23) by choosing $\mathbf{f}_{|E}={\bf 0}$. In addition, using the results in [4], we infer that $\Sigma^{f}_{h}(E)$ and $\Sigma^{0}_{h}(E)$ are completely characterised once the quantities $\mathbf{c}_{i}$ and $d_{i}$ are given, for $i=1,\ldots,n_{E}-1$. Therefore, the dimension of both the spaces is $3n_{E}-3$.

3.2 The local bilinear forms

Given $E\in{\mathcal{T}}_{h}$, we first notice that, for every {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E) and $\mathbf{v}_{h}\in U_{h}(E)$, the term

[TABLE]

is computable from the knowledge of the degrees of freedom. Therefore, there is no need to introduce any approximation in the structure of the terms b({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}^{0},\mathbf{u}) and b({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0},\mathbf{v}) in problem (7). Instead, the term

[TABLE]

is not computable for a general couple ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h})\in\Sigma_{h}(E)\times\Sigma_{h}(E). As usual in the VEM approach (see [10], for instance), we then need to introduce a suitable approximation $a_{E}^{h}(\cdot,\cdot)$. We first define the projection operator

[TABLE]

where ${\mathcal{P}}_{\ast}(E)^{2\times 2}_{s}$ is a suitable space of polynomial symmetric tensors. We then set

[TABLE]

where $s_{E}(\cdot,\cdot)$ is a suitable stabilization term. We propose the following:

[TABLE]

where $\kappa_{E}$ is a positive constant to be chosen (for instance, any norm of $\mathbb{D}_{|E}$). A variant of (28) is provided by

[TABLE]

Morevover, we will make two different choices for ${\mathcal{P}}_{\ast}(E)^{2\times 2}_{s}$, namely:

[TABLE]

and

[TABLE]

Remark 3**.**

We remark that choice (31) leads to a VEM projection onto a reacher space than choice (30). Although we do not have any improvement in the convergence rate, we nonetheless expect more accurate discrete solutions when selecting (31) instead of (30). This behaviour is generally confirmed by the numerical experiments presented in Section 5.

3.3 The local term a_{E}(\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt})

In order to discretise F({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}), see (7), we need to consider the term, see (8):

[TABLE]

To this aim, we will proceed in two steps:

we first choose \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{h,f}, which may be considered as a suitable approximation of \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}; 2. 2.

we then consider a $\widetilde{a}_{E}^{h}(\cdot,\cdot)$, which may be considered as a suitable approximation of $a_{E}(\cdot,\cdot)$, and we finally form \widetilde{a}_{E}^{h}(\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{h,f},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}) to discretise (32).

In particular, recalling that $\mathbf{f}_{|E}\in{\mathcal{P}}_{0}(E)$, we define \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{h,f} as:

[TABLE]

Above, for a given vector $\mathbf{v}$, $(\mathbf{v})_{i}$ denotes its $i$-th component. Furthermore, we set

[TABLE]

Remark 4**.**

A different option could be to select \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{h,f}\in\Sigma^{f}_{h}(E), exploiting (23) and Remark 2. For instance, we may set:

[TABLE]

and compute $\mathbf{c}_{n_{E}}$ and $d_{n_{E}}$ according to (23). However, other choices can be made (more balanced among the polygon edges). To approximate (32), we could consider:

[TABLE]

where $a_{E}^{h}(\cdot,\cdot)$ is defined as in Section 3.2.

3.4 The discrete scheme

We are now ready to introduce the discrete scheme. We introduce global approximation spaces, see (3), (21) and (22):

[TABLE]

Furthermore, given a local approximation of $a_{E}(\cdot,\cdot)$, see (27), we set

[TABLE]

The method we consider is then defined by

[TABLE]

Above, $F_{h}(\cdot)$ is given by, see (8), (33) and (34):

[TABLE]

while $G_{h}(\mathbf{v}_{h})$ reads, see (8):

[TABLE]

Introducing the bilinear form ${\mathcal{A}}_{h}:(\Sigma_{h}^{0}\times U_{h}^{0})\times(\Sigma_{h}^{0}\times U_{h}^{0})\to\mathbb{R}$ defined by

[TABLE]

problem (39) can be written as

[TABLE]

4 Theoretical analysis using a mesh-dependent norm

A significant part of the theoretical analysis follows the guidelines developed in [4]. Hence, in many points we limit to state the results. However, we develop here an error analysis using a mesh-dependent norm for the stresses, which turns out to be a flexible tool in our case. In addition, we remark that the inf-sup condition of Proposition 4.9 is new, and therefore its proof is given in full details.

According to the assumption $\mathbf{f}_{|E}\in{\mathcal{P}}_{0}(E)$, we select the locally self-equilibrated solution \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f} as

[TABLE]

Therefore, we have \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}=\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{h,f}, see (33). As a consequence, we get (cf. (8), (40) and (41))

[TABLE]

4.1 Stability conditions

We introduce the following mesh-dependent quantity in $\Sigma_{h}:=\prod_{E\in{\mathcal{T}}_{h}}\Sigma_{h}(E)$:

[TABLE]

It is easily seen that (46) defines a norm on $\Sigma_{h}$, see (16) along with (19) and (20). Moreover, we have the following lemma.

Lemma 4.1**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, it holds:

[TABLE]

Proof.

Fix $E\in{\mathcal{T}}_{h}$. Using Lemma 5.1 of [4] we get

[TABLE]

Recalling that \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in RM(E), an integration by parts, the Agmon inequality (see for instance [2]) and an inverse estimate on polygons (see Lemma 6.3 of [8]) give:

[TABLE]

Hence it holds

[TABLE]

Combining (48) with (50) we obtain

[TABLE]

from which the first estimate of (47) follows.

We now prove the second estimate in (47). For every edge $e$ in $E$, let us denote with $b_{e}$ the function defined on $\partial E$ such that: on $e$ it is the quadratic bubble with $\sup\,b_{e}=1$, and $b_{e}=0$ elsewhere. One has:

[TABLE]

where we have also used a 1D inverse estimate. Hence, after exploiting a scaled trace inequality (cf. [8]), we get

[TABLE]

Now, using the technique of Lemma 6.3 in [8], we obtain

[TABLE]

Therefore, we have:

[TABLE]

from which the second inequality in (47) follows.

∎

As it is well-known (see for instance [12] or [13]), stability for problems with the format as in (42), is implied by the satisfaction of two conditions: the ellipticity-on-the-kernel condition, and the inf-sup condition. As far as the first one is concerned, we first prove that it holds:

[TABLE]

Above $K_{h}\subseteq\Sigma^{0}_{h}$ is the discrete kernel, defined by:

[TABLE]

Estimate (55) is a consequence of the following stronger result.

Lemma 4.2**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. Then it holds:

[TABLE]

Proof.

By the norm definition (5) and the space definition (21), we immediately infer

[TABLE]

Using (47) and standard VEM arguments as in [10], we get

[TABLE]

Estimate (55) follows from (58) and (59).

∎

From Lemma 4.2 and (47) we infer the coercivity property:

[TABLE]

The following lemma will be useful to prove the inf-sup condition.

Lemma 4.3**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled, and fix $E\in{\mathcal{T}}_{h}$. Take any $\mathbf{w}\in E$, node of ${\mathcal{T}}_{h}$. Then it holds:

[TABLE]

where ${\mathcal{M}}_{h}(E)$ is the set of edge mid-points of $E$.

Proof.

Denote with $\{\mathbf{w}=\mathbf{w}_{1},\mathbf{w}_{2}\ldots,\mathbf{w}_{n_{E}}\}$ the set of nodes for ${\mathcal{T}}_{h}$ on $\partial E$, ordered counter-clockwise. Furthermore, let $e_{i}=[\mathbf{w}_{i},\mathbf{w}_{i+1}]$ ($i=1,\ldots,n_{E}$) be the edges of $E$ (here we have set $\mathbf{w}_{n_{E}+1}=\mathbf{w}_{1}$), and let $\mathbf{m}_{i}$ be the midpoint of $e_{i}$. Fix $\mathbf{v}_{h}\in U_{h}(E)$; by using the Cavalieri-Simpson rule, we get

[TABLE]

Due to assumption $\mathbf{(A2)}$, we get

[TABLE]

We now notice that, since each component of $\mathbf{v}_{h}$ is a piecewise linear and continuous function on $\partial E$, it follows that for $i=2,\ldots,n_{E}$, the quantity $\mathbf{v}_{h}(\mathbf{w}_{i})$ is uniquely determined by $\mathbf{v}_{h}(\mathbf{w}_{1})$ and $\{\mathbf{v}_{h}(\mathbf{m}_{1}),\ldots,\mathbf{v}_{h}(\mathbf{m}_{i-1})\}$. Indeed, a direct computation shows that

[TABLE]

Hence, assumption $\mathbf{(A1)}$ and $\mathbf{(A2)}$ implies that

[TABLE]

Estimate (61) now follows from a combination of (63) and (64).

∎

Remark 5**.**

It is easy to see that, if $n_{E}$ is odd, then the values $\mathbf{v}_{h}(\mathbf{w}_{i})$ ($1\leq i\leq n_{E}$) can be determined without using $\mathbf{v}_{h}(\mathbf{w})$, but only the midpoint values $\{\mathbf{v}_{h}(\mathbf{m}_{1}),\ldots,\mathbf{v}_{h}(\mathbf{m}_{n_{E}})\}$.

We now prove the following local inf-sup condition.

Lemma 4.4**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. Then, there exists $\beta>0$ such that

[TABLE]

where the bilinear form $b_{E}(\cdot,\cdot)$ is defined by (cf. (8)):

[TABLE]

Proof.

Fix $\mathbf{v}_{h}\in U_{h}(E)$ and choose {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E) such that:

[TABLE]

Due to (16), the above choice is admissible. Thus, using also the mid-point rule, we get:

[TABLE]

Applying to $\mathbf{v}_{h}\cdot\mathbf{n}$ the same argument as in (62)-(63), from (68) we infer

[TABLE]

Above, $\mathbf{n}_{i}^{-}$ and $\mathbf{n}_{i}^{+}$ are the two normals to the edges which share $\mathbf{w}_{i}$ as a common vertex. We now notice that, due to assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, there exists a node $\mathbf{w}_{m}$ ($1\leq m\leq n_{E}$), for which $|(\mathbf{v}_{h}\cdot\mathbf{n}_{m}^{-})(\mathbf{w}_{m})|^{2}+|(\mathbf{v}_{h}\cdot\mathbf{n}_{m}^{+})(\mathbf{w}_{m})|^{2}\approx|\mathbf{v}_{h}(\mathbf{w}_{m})|^{2}$. Therefore, from (69) we get

[TABLE]

Applying Lemma 4.3 we thus infer

[TABLE]

Furthermore, from (67) we immediately get

[TABLE]

Estimate (65) is now a consequence of (71) and (72).

∎

Introducing the local discrete kernel $K_{h}(\partial E)$ defined by:

[TABLE]

Lemma 4.4 implies the following result (cf. [12]).

Lemma 4.5**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. Then, there exists $\beta>0$ such that

[TABLE]

We are now ready to prove the Proposition:

Proposition 4.6**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, there exists a linear operator $\pi_{E}\,:\,\Sigma(E)\to\Sigma_{h}(E)$ such that:

[TABLE]

Proof.

Fix {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma(E). Due to Lemma 4.4, the linear system in the first line of (75) is solvable, and two solutions differ up to an element of $K_{h}(\partial E)$, cf. (73). To prove the continuity estimate in (75), let us take \pi_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma_{h}(E) as the solution which minimizes ||(\pi_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt})\mathbf{n}||_{0,\partial E}. From (74) and the first of (75), we thus get:

[TABLE]

By recalling (66), a (scaled) duality estimate and a 1D inverse estimate for piecewise linear polynomials, give:

[TABLE]

Therefore, from (76) and (77) we obtain

[TABLE]

The continuity estimate in (75) now follows from a trace estimate and definition (46).

∎

We notice that {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma^{0}(E) implies \pi_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma^{0}_{h}(E), cf. (4) and (21). Indeed, by definition (16), \mathop{\bf div}\nolimits(\pi_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt})\in RM(E). Since $RM(E)\subseteq U_{h}(E)$, we can take \mathbf{v}_{h}=\mathop{\bf div}\nolimits(\pi_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}) in (75), to obtain (using also the integration by parts):

[TABLE]

This observation, together with Proposition 4.6, immediately leads to the following result.

Corollary 4.7**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, there exists a linear operator $\pi_{E}\,:\,\Sigma^{0}(E)\to\Sigma_{h}^{0}(E)$ such that:

[TABLE]

Recalling (3), we define the linear operator $\pi_{h}\,:\,\Sigma^{0}\to\Sigma_{h}^{0}$ by adding the local contributions $\pi_{E}$, i.e.:

[TABLE]

Obviously, Corollary 4.7 and (8) give:

Proposition 4.8**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, there exists a linear operator $\pi_{h}\,:\,\Sigma^{0}\to\Sigma_{h}^{0}$ such that:

[TABLE]

With Proposition 4.8 at hand, the following inf-sup condition is easily proved (it is nothing but Fortin’s trick, see [12]).

Proposition 4.9**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, we have:

[TABLE]

Proof.

Fix $\mathbf{v}_{h}\in U^{0}_{h}$. Using (82) and the inf-sup condition for the continuous problem (see (10)), we infer

[TABLE]

∎

4.2 An interpolation operator for the stresses

We now recall the interpolation operator introduced in [4]. We first set, given $r>2$:

[TABLE]

To continue, we locally define the operator ${\mathcal{I}}_{E}:W^{r}(E)\to\Sigma_{h}(E)$ as follows. Given {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in W^{r}(E), {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma_{h}(E) is determined by:

[TABLE]

where

[TABLE]

If \tau$$\tau$$\tau is not sufficiently regular, the integral in the right-hsnd side of (86) is intended as a duality between $W^{-\frac{1}{r},r}(\partial E)^{2}$ and $W^{\frac{1}{r},r^{\prime}}(\partial E)^{2}$. If \tau$$\tau$$\tau is a regular function, the above condition is equivalent to require:

[TABLE]

The global interpolation operator ${\mathcal{I}}_{h}:W^{r}(\Omega)\to\Sigma_{h}$ is then defined by simply glueing the local contributions provided by ${\mathcal{I}}_{E}$. More precisely, we set ({\mathcal{I}}_{h}\tau)_{|E}:={\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{|E} for every $E\in{\mathcal{T}}_{h}$ and {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in W^{r}(\Omega). It can be proved, see [4], that the commuting diagram property:

[TABLE]

holds true, $\Pi_{RM}$ being the $L^{2}$-projection operator onto the space of local rigid body motions. Furthermore, the following estimates have been proved in [4].

Proposition 4.10**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, for the interpolation operator ${\mathcal{I}}_{E}$ defined in (88), the following estimates hold:

[TABLE]

[TABLE]

4.3 Error estimates

We denote with ${\mathcal{P}}_{0}({\mathcal{T}}_{h})$ the space of piecewise constant functions with respect to the given mesh ${\mathcal{T}}_{h}$.

Once one has the stability conditions of estimate (60) and Proposition 4.9, along with the interpolation estimates of Proposition 4.10, an error analysis can be derived using the techniques of [15] or [4]. Indeed, one can prove:

Proposition 4.11**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. The following error estimate holds:

[TABLE]

*where ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u})\in\Sigma^{f}\times U^{0}/H is the solution to the continuous Problem 13, and ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h},\mathbf{u}_{h}) is such that {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h}={\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0}_{h}+\widehat{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h,f}, ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h},\mathbf{u}_{h})\in\Sigma^{f}_{h}\times U^{0}_{h}/H being the solution to the discrete problem 39. Furthermore, C(\Omega,{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u}) is independent of $h$ but depends on the domain $\Omega$ and on the Sobolev regularity of \sigma$$\sigma$$\sigma and $\mathbf{u}$. *

Proof.

We take {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I}\in\Sigma_{h}^{f} defined as {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I}={\mathcal{I}}_{h}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}, see section 4.2. Due to the splitting {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}={\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0}+\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}, with {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0}\in\Sigma^{0}, we have

[TABLE]

From (89), we get {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0}_{I}\in\Sigma^{0}_{h} and \widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{I,f}\in\Sigma^{f}_{h}. We also take $\mathbf{u}_{I}\in U^{0}_{h}$ as the usual VEM interpolant of $\mathbf{u}$, see for example [10] or [8].

We now form ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h}-{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I},\mathbf{u}_{h}-\mathbf{u}_{I})\in\Sigma_{h}^{0}\times U_{h}^{0}. Then, using the ellipticity-on-the-kernel condition of estimate (60) and the inf-sup condition of Proposition 4.9, there exists ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h},\mathbf{v}_{h})\in\Sigma_{h}^{0}\times U_{h}^{0} such that (see [12] and [13], for instance):

[TABLE]

and

[TABLE]

We have, considering the splitting {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}={\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}^{0}+\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f} and (93), and using both (7) and (43), together with (45):

[TABLE]

For both the choices (30) and (31), the term $T_{1}$ can be treated using the techniques of [4], to obtain:

[TABLE]

where {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{\pi} is the $L^{2}$-projection of \sigma$$\sigma$$\sigma onto ${\mathcal{P}}_{0}({\mathcal{T}}_{h})^{2\times 2}_{s}$.

Regarding $T_{2}$, using the Agmon’s trace inequality (see for example [2]), one has:

[TABLE]

Term $T_{3}$ can be treated using standard trace inequalities, to obtain:

[TABLE]

To estimate the term $T_{4}$ we first recall that, see (45):

[TABLE]

If $\Pi_{E}$ is selected according with (31), then \Pi_{E}\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}=\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}, and $T_{4}$ vanishes. If $\Pi_{E}$ is selected according with (30), we have \Pi_{E}\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}=0 (cf. (44)). Then we get, also using (47):

[TABLE]

A direct computation taking into account (44) gives ||\widehat{{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}}_{f}||_{0,E}\lesssim h_{E}||\mathbf{f}||_{0,E}. Therefore, we obtain

[TABLE]

From (95), (96), (97), (98), (99) and (102), we get:

[TABLE]

Using (94), standard approximation results and estimates (90)-(91), we infer:

[TABLE]

where C(\Omega,{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u}) is independent of $h$ but depends on the domain $\Omega$ and on the Sobolev regularity of \sigma$$\sigma$$\sigma and $\mathbf{u}$. We now use the triangle inequality and the estimate

[TABLE]

Exploiting (104), again standard approximation results and (90), we thus obtain:

[TABLE]

∎

5 Numerical results

The present section is devoted to validation of the proposed dual hybrid methods. First, in section 5.1 convergence and accuracy are numerically assessed on a couple of benchmarks having a closed-form solution. Subsequently, in section 5.2, an elastic problem stemming from a simple electromechanical application is considered, proving applicability of the method to the analysis and simulation of real structures. In all presented tests reference is made for comparison to the standard displacement-based linear Virtual Element Method detailed in [3].

5.1 Convergence and accuracy assessment

A set of two boundary value problems on the unit square domain $\Omega=[0,1]^{2}$ is considered, for which an analytical solution is available [20]. Material parameters are assigned in terms of Lamé constants $\lambda=1$, $\mu=1$, assuming plane strain and homogeneous isotropy conditions. The tests are defined by choosing a required solution and deriving the corresponding body force $\mathbf{f}$, as reported in the following:

• $\bullet$

Test $a$

[TABLE]

• $\bullet$

Test $b$

[TABLE]

Test $a$ has Dirichlet non-homogeneous boundary conditions, zero body force and a polynomial solution; Test $b$ has Dirichlet homogeneous boundary conditions, trigonometric body force and a trigonometric solution.

The tests are run on two sets of structured [resp. unstructured] square, hexagon, and concave quadrilateral [resp. triangle, quadrilateral, and Voronoi] simple polygonal meshes, each type being represented and labeled in Fig. 1, for a sequence of five uniform mesh refinements.

Numerical solutions for the above tests are sought with the proposed dual hybrid method, with the two versions of projection operator Eq. (30) [resp. Eq. (31)], which will be labeled DH P[math] [resp. DH P$1$]. For comparison purposes, a further numerical solution with the linear displacement based VEM presented in [3] is computed and labeled $\rm DISP$.

Convergence rate and accuracy level are investigated computing the following error quantities:

• $\bullet$

Discrete relative error quantity for the stress field:

[TABLE]

• $\bullet$

Discrete relative weighted error quantity for the inter-element traction field:

[TABLE]

where $\mathbf{n}$ is one outward unit normal to the edge $e$ chosen once and for all. Quantities $\mathbf{t}_{\mathbf{n},h}$ are the average of the two contributions stemming from the two elements adjacent to edge $e$.

• $\bullet$

Discrete $H^{1}$-type error quantity for the inter-element displacement field:

[TABLE]

where $\mathbf{e}$ is the unit tangent to the edge $e$ chosen once and for all.

Inspection of Figs. 2-5 confirms expected convergence rates for the three compared methods. In terms of accuracy, for both tests and for all mesh types, dual hybrid virtual element methods DH P0 and DH P1 outperform displacement based virtual element method DISP, for stress and inter-element traction field, respectively. Comparatively, DH P1 shows the highest edge on Test b (cf. also Remark 3). It is noted that the three stress [resp. traction] fields coincide for Test a adopting triangles. Figs. 6-7 show that DH P1 and DISP methods are fairly comparable in terms of displacement field accuracy, for both tests and for all mesh types, with some selected cases in which DISP presents the lowest error levels while DH P0 the highest ones. It is noted that the three displacement fields coincide for Test a adopting structured quadrilaterals, and for both tests adopting triangles.

5.2 Structural application: folded-beam suspension

A representative section of a typical microelectromechanical system (MEMS) [26], consisting of two bulky portions connected with four slender beams (see Fig. 8), is considered, as a structural application on which we test the VEM methods described in the previous section. As an electromechanical plate-like device, a two-dimensional linear elastic analysis under plane stress assumption is carried out. Geometry, dimensions, boundary conditions and loading are represented in Fig. 8. Material parameters are $E=60$ GPa, $\nu=0.22$; edge traction $q=10^{-2}$ N/m. The relevant half domain is meshed with triangles, quadrilaterals, and Voronoi polygons as previously done. The latter spatial discretization makes use of non-convex polygons in the zones surrounding re-entering corners, which, given the particular geometry under consideration, further enlights the broader mesh capability offered by VEM methods in respect with standard FEM discretizations.

Progressively finer meshes are considered for DISP, DH P0, DH P1 method, respectively, while a reference solution is computed with CPE4H hybrid element implemented in COMSOL [19] on a very fine mesh. Results in terms of relative error on the vertical displacement of target node $A$ (cf. Fig. 8) against total number of degrees of freedom are shown in Fig. 9, which confirms the efficiency of the proposed dual hybrid method as a tool for structural analysis.

6 Conclusions

We have presented a Virtual Element Method for 2D linear elastic problems, based on a dual hybrid variational formulation. The scheme offers two variants, which differ from each other according to the choice of the VEM stress projection. A stability and convergence analysis has been developed, and several numerical tests have been performed, confirming the theoretical predicitions. Our study shows that dual hybrid VEMs represents a valid alternative to standard displacement-based VEM schemes, especially if one is interested in an accurate description of the stress field. We finally remark that hybrid elasticity methods have been used, in the FEM framework, to tackle Structural Mechanics problems (e.g. plate problems, see [21]): our VEM scheme might be fruitfully employed within that context as well.

Aknowledgements

E.A. gratefully acknowledges the partial financial support of the University of Rome Tor Vergata Mission Sustainability Programme through project SPY-E81I18000540005.