3D Adaptive VEM with Stabilization-Free a Posteriori Error Bounds

Abstract

The present paper extends the theory of Adaptive Virtual Element Methods (AVEMs) started in (Beirão da Veiga et al. in SIAM J Numer Anal 61(2):457–494, https://doi.org/10.1137/21M1458740, 2023) to the three-dimensional meshes showing the possibility to bound the stabilization term by the residual-type error estimator. This new bound establishes the equivalence between a stabilization-free residual-type a posteriori error estimator and the energy error, enabling the formulation of a 3D AVEM algorithm and providing the necessary results to prove its convergence. Following the recent studies for the bi-dimensional case, we investigate the case of tetrahedral elements with aligned edges and faces. We believe that the AVEMs can be an efficient strategy to address the mesh conforming requirements of standard three-dimensional Adaptive Finite Element Methods (AFEMs), which typically extend the refinement procedure to non-marked mesh cells. Indeed, numerical tests show that this method can reduce the number of three-dimensional cells generated in the refinement process up to $30\%$ with respect to standard AFEMs, for a given error threshold.

1 Introduction

Adaptive methods for approximating the solutions of partial differential equations are fundamental tools in scientific computing and engineering, whose mathematical background has been studied for decades [11, 15, 22]. The standard adaptive methods are based on the iteration of the classical paradigm

$$\begin{aligned} \textsf{SOLVE}\rightarrow \textsf{ESTIMATE}\rightarrow \textsf{MARK}\rightarrow \textsf{REFINE}, \end{aligned}$$

to achieve at each loop a better approximation of the solution. The main idea behind the adaptivity resides in the definition of an a posteriori error estimator which indicates which elements bring the major error contribution, and then need to be refined. In standard Adaptive Finite Element Methods (AFEMs), mesh conformity requirements often necessitate the refinement of adjacent elements for purely geometric reasons when refining the marked elements. On the other hand, since their introduction [3, 4], Virtual Element Methods (VEMs) offer a valuable approach to refine only the marked elements without requiring additional refinement of neighboring elements, thanks to the presence of aligned edges and faces on the tessellation elements. However, the extension of the adaptive FEM theory to the VEM is not trivial. The first contraction results for an adaptive VEM algorithm has been recently proposed in [6] and the study of its optimality in [5] in the lowest order case for triangular elements admitting the presence of aligned edges. The extension to the higher-order case has been introduced in [14]. In VEM framework, the definition of an a-posteriori error estimator has to take into account the presence of a stabilization term. Indeed, the original VEM a posteriori error analysis showed in [13] presented bounds of the type

(1.1)

where is the error in the energy norm between the solution u of the continuous problem and the discrete solution $u_\mathcal {T}$ computed, $\eta _\mathcal {T}(u_\mathcal {T})$ is residual type estimator and $S_\mathcal {T}(u_\mathcal {T})$ is a stability term that needs to be added to the discrete bilinear form to get the coercivity. In particular, as remarked in [18], $S_\mathcal {T}(u_\mathcal {T})$ presents two main complications: (i) the stabilization term is not robust with respect to the polynomial degree of the method, and (ii) one should be really careful in the choice of $S_\mathcal {T}(u_\mathcal {T})$ since the refinement procedure may produce elements with very small edges and faces. One of the first studies of an a posteriori error bound that does not involve the stability term, under specific assumptions, can be traced back to [7]. Two of the possible strategies to overcome the presence of the stability term in (1.1) are the analysis of the recently introduced stabilization-free VEM schemes [8,9,10] or the definition of stabilization-free a posteriori error bounds for the classical VEM. In this work we investigate the latter scenario. The bi-dimensional first-order VEM bound of the type

$$\begin{aligned} S_\mathcal {T}(u_\mathcal {T})\,\lesssim \eta _\mathcal {T}(u_\mathcal {T}), \end{aligned}$$

is due to [6]. The purpose of this paper is to extend this result to the three-dimensional case and, consequently, to get

This stabilization-free a posteriori error bound brings to the definition of an adaptive algorithm, whose convergence result retraces the one of [6].

The main novelty that arises in the higher-dimensional case is a new proof of the scaled Poincaré inequality, which results essential to prove the stabilization-free a posteriori error bound. In this paper we have decided to report the main propositions and theorems introduced in [6] for the sake of readability, highlighting the differences which occur. Furthermore, we focus on the case of tetrahedral elements with aligned faces and aligned edges, since a refinement strategy for a general three-dimensional polyhedron preserving quality is an open problem. Some practical strategies of polyhedral refinement has been proposed in the last years, [2].

The paper is structured as follows. In Sect. 2 the continuous and the discrete problems are presented. We give here the definition of the global index of a node for the three dimensional case and the two preliminary properties on which the rest of the paper is based, i.e. the scaled Poincaré inequality and the Clément interpolation estimate. Section 3 presents the stabilization-free a posteriori error bound and Sect. 4 discusses the convergence of the method. Finally, Sect. 5 is devoted to the numerical experiments.

2 The Problem and the Discretization

We focus on the second-order Dirichlet boundary-value problem on a polyhedral domain $\varOmega \subset \mathbb {R}^3$,

$$\begin{aligned} {\left\{ \begin{array}{ll} -\nabla \cdot (\mathcal {K}\nabla u) + c u = f & \text { in }\varOmega ,\\ u = 0 & \text { on }\partial \varOmega , \end{array}\right. } \end{aligned}$$

(2.1)

where the data, $\mathcal {D}= (\mathcal {K}, c, f)$, are $\mathcal {K}\in (L^\infty (\varOmega ))^{3\times 3}$ symmetric and positive-definite in $\varOmega$, $c\in L^\infty (\varOmega )$ non negative in $\varOmega$, and $f\in L^2(\varOmega )$. For simplicity, we consider homogeneous Dirichlet boundary conditions. The extension to more general boundary conditions is more technical and follows similar arguments. From problem (2.1), we can derive its variational formulation, which reads as

$$\begin{aligned} {\left\{ \begin{array}{ll} \text {find }u \in \mathbb {V} {:}{=}H^1_{0}(\varOmega ) & \text { such that}\\ \mathcal {B}(u,v) = (f,v), & \forall \;v \in \mathbb {V}, \end{array}\right. } \end{aligned}$$

(2.2)

where $(\cdot ,\cdot )$ is the scalar product in $L^2(\varOmega )$ and $\mathcal {B}(u,v) {:}{=}a(u,v) + m(u,v)$ is the bilinear form associated with problem, i.e.

$$\begin{aligned} a(u,v) {:}{=}(\mathcal {K}\; \nabla u , \nabla v), \quad \quad \quad m(u,v) {:}{=}( c\; u,v). \end{aligned}$$

From the assumption on the data, $\mathcal {B}$ is continuous and coercive, and problem (2.2) admits the existence and uniqueness of the solution u. We indicate the energy norm as , which satisfies

(2.3)

where is the $H^1$-seminorm and $c_\mathcal {B}$ and $c^\mathcal {B}$ are two suitable strictly positive constants.

In the perspective of presenting a discrete formulation of the variational problem, we first introduce an initial conforming partition $\mathcal {T}_0$ of $\overline{\varOmega }$ made of tetrahedra which satisfies the initial conditions given in [23]. Let $\mathcal {T}$ be any successive partition obtained from $\mathcal {T}_0$ by the refinement of some elements via the three-dimensional newest-vertex bisection (NVB) algorithm [22].

Given a tetrahedron $E\in \mathcal {T}_0$, this refinement consists in selecting an edge $\tilde{e}$ of E and connecting its midpoint with the two remaining vertices “opposite” to $\tilde{e}$. The two new edges so built are part of the new triangular face formed in the refinement procedure that splits the tetrahedron into two elements (see Fig. 2 for a graphical representation).

We remark that we have decided to use the term tetrahedron to denote the three-dimensional simplex, and we denote by tetrahedral-shape elements all the elements obtained during the refinement process that might contain aligned edges and faces.

Given a tetrahedral-shape element $E \in \mathcal {T}$, we define a straight side of E each one of the six edges of the smallest simplex that contains E. Moreover, let $\mathcal {N}_E$ denote the collection of all the nodes of E. We refer to each of the four endpoints of the straight side of E as proper node of E, or simply vertex. These four vertices are collected in the set $\mathcal {P}_E$, namely the set of proper nodes of E. Finally, we indicate by $\mathcal {H}_E$ the set of all the nodes sitting on the boundary of E that are not propers for E, i.e., $\mathcal {H}_E:= \mathcal {N}_E \setminus \mathcal {P}_E$. Note that, in the classic case where E is a tetrahedron, the vertices and the nodes of E coincide. Specifically, $\mathcal {N}_E = \mathcal {P}_E$ and $\mathcal {H}_E = \emptyset$. See Fig. 1 for an example.

At a global level, we denote by $\mathcal {N}$ the set of all the nodes of the partition, we define the proper nodes set as $\mathcal {P}:= \left\{ \varvec{x}\in \mathcal {N}: \varvec{x}\in \mathcal {P}_E, \ \forall E \text { containing }\varvec{x}\right\}$, and we call hanging the nodes that are not proper nodes, i.e. $\mathcal {H}= \mathcal {N}\setminus \mathcal {P}$. We recall that the presence of hanging nodes at the boundary of the elements E, can be handled by the VEM theory, which treats tetrahedra with hanging nodes as polyhedra with aligned faces and edges. We denote by $\mathcal {F}_E$ the set of faces of each element $E \in \mathcal {T}$, by $\partial F$ the boundary of the face F, and by $\mathcal {E}_F$ the set of edges of any face F. In the paper, we set $h_E = |E|^{1/3}$, where |E| denotes the volume of E, and $h=\max _{E\in \mathcal {T}} h_E$.

A crucial concept, firstly introduced in [6] for the bi-dimensional case and extended in [14], is the global index of a node. To define it, we recall that in the NVB refinement strategy, a hanging node $\varvec{x}\in \mathcal {H}$ is always obtained by bisecting an edge of the tetrahedron. We denote by $\{\varvec{x'}, \varvec{x''}\}$ the set of the endpoints of the refined edge for which $\varvec{x}$ is the midpoint.

Definition 1

(global index of a node) Given a node $\varvec{x}\in \mathcal {N}$, we define its global index $\lambda$ as:

• $\lambda (\varvec{x})=0,$ if $\varvec{x}\in \mathcal {P}$;

• $\lambda (\varvec{x})= \max \{\lambda (\varvec{x'}), \lambda (\varvec{x''})\}+1$, if $\varvec{x}\in \mathcal {H}$.

An example of the evolution of the global index of a node, after three refinements is shown in Fig. 2.

Fig. 1

The light blue element on the left has five nodes. A, B, C, D are the proper nodes of the element, defining the six straight-sides, i.e. $\overline{\text {AB}}, \overline{\text {BC}}$, $\overline{\text {CD}}$, $\overline{\text {DA}}$, $\overline{\text {BD}}$, $\overline{\text {AC}}$. The node E lies on the straight-side $\overline{\text {DA}}$ and it is a hanging node for the element. This element, for instance, appears in the partition shown on the right side of the figure, which can be obtained through the refinement sequence depicted in the Fig. 2

Fig. 2

Three successive refinements of a tetrahedron. Blue lines represent the edges added at each refinement. Black points and red points represent the proper and the hanging nodes, respectively. The numbers are the global indices

In the analysis presented, we suppose that the largest global index of the partition does not increase so much. More specifically, we have the following assumption.

Assumption 1

We assume the existence of a constant $\varLambda \ge 1$ such that

$$\begin{aligned} \varLambda _\mathcal {T}{:}{=}\max _{\varvec{x}\in \mathcal {N}}\lambda (\varvec{x})\le \varLambda . \end{aligned}$$

We call this type of partition $\varLambda$-admissible.

Remark 1

The case $\varLambda =0$ is not included in Assumption 1 and it is not considered in the following analysis. Indeed, if $\varLambda =0$ the partition does not admit hanging nodes, then the VEM scheme reduces to the FEM one.

Remark 2

If a partition $\mathcal {T}$ is $\varLambda$-admissible, the number of hanging nodes on each side of an element can be upper bounded by $2^\varLambda -1$ hanging nodes (see also [6, Remark 2.3]).

2.1 VEM Spaces and Projectors

The construction of the three-dimensional VEM spaces is an extension of the bi-dimensional ones [1, 3, 4]. Let E be a tetrahedral-shape element in the partition $\mathcal {T}$, and F a face in $\mathcal {F}_E$. On F, we define the enhanced bi-dimensional space of order one, starting from the boundary of F

$$\begin{aligned} \mathbb {V}_{\partial F}{:}{=}\{ v \in C^0(\partial F): v|_e \in \mathbb {P}_1(e),\; \forall e\in \mathcal {E}_F\}, \end{aligned}$$

and extending it in the interior of F,

$$\begin{aligned} \begin{aligned} \mathbb {V}_F {:}{=}\{v \in H^1(F): v|_{\partial F} \in \mathbb {V}_{\partial F}, \;\varDelta v|_F \in \mathbb {P}_1(F),\; \\ \int _F \left( v- \varPi ^\nabla _F v\right) {{p_1}}=0 ,\; \forall {{p_1}} \in \mathbb {P}_1(F)\}, \end{aligned} \end{aligned}$$

where the projector $\varPi ^\nabla _F: H^1(F)\rightarrow \mathbb {P}_1(F)$ is such that $\forall v \in H^1(F)$

$$\begin{aligned} \int _F \nabla \left( \varPi ^\nabla _F v - v \right) \cdot \nabla {{p_1}} = 0\quad \forall {{p_1}} \in \mathbb {P}_1(F),\quad \quad \int _{\partial F}\left( \varPi ^\nabla _F v -v\right) =0. \end{aligned}$$

In the element E, we define the space $\mathbb {V}_E$ such that it reduces to $\mathbb {V}_F$ on each face of $\mathcal {F}_E$. In particular, we define on the boundary of E,

$$\begin{aligned} \mathbb {V}_{\partial E}{:}{=}\{v \in C^0(\partial E): \; v|_F \in \mathbb {V}_F,\; \forall F \in \mathcal {F}_E\}, \end{aligned}$$

and the enhanced local space

$$\begin{aligned} \begin{aligned} \mathbb {V}_E {:}{=}\{v \in H^1(E): \; v|_{\partial E}\in \mathbb {V}_{\partial E}, \;\varDelta v|_E \in \mathbb {P}_1(E), \; \\ \int _E \left( v - \varPi ^\nabla _ E\right) p_1 =0, \; \forall p_1 \in \mathbb {P}_1(E)\},\end{aligned} \end{aligned}$$

(2.4)

in which $\varPi ^\nabla _E: H^1(E)\rightarrow \mathbb {P}_1 (E)$ is defined as $\forall v \in H^1(E)$

$$\begin{aligned} \int _E \nabla \left( \varPi ^\nabla _E v - v \right) \cdot \nabla w = 0\quad \forall w \in \mathbb {P}_1(E), \quad \quad \quad \int _{\partial E}\left( \varPi ^\nabla _E v -v\right) =0. \end{aligned}$$

(2.5)

We denote the global space as

$$\begin{aligned} \mathbb {V}_\mathcal {T}{:}{=}\{v \in \mathbb {V}: v|_E \in \mathbb {V}_E,\; \forall E \in \mathcal {T}\}, \end{aligned}$$

the space of piecewise polynomials of order one on $\mathcal {T}$ as

$$\begin{aligned} \mathbb {W}_\mathcal {T}{:}{=}\{w \in L^2(\varOmega ): w|_E\in \mathbb {P}_1 (E),\; \forall E\in \mathcal {T}\}, \end{aligned}$$

and the subspace $\mathbb {V}^0_\mathcal {T}$, i.e.

$$\begin{aligned} \mathbb {V}^0_\mathcal {T}{:}{=}\mathbb {V}_\mathcal {T}\cap \mathbb {W}_\mathcal {T}, \end{aligned}$$

(2.6)

which plays a crucial role in the proofs.

Remark 3

Notice that in general $\mathbb {W}_\mathcal {T}\not \subset \mathbb {V}_\mathcal {T}$, indeed $\mathbb {V}_\mathcal {T}\subset C^0(\varOmega )$, while if $v\in \mathbb {W}_\mathcal {T}$ in general $v\not \in C^0(\varOmega )$, but only piece-wise polynomial on each element E of the partition $\mathcal {T}$.

Let us define the projection operator $\varPi ^\nabla _\mathcal {T}: \mathbb {V}_\mathcal {T}\rightarrow \mathbb {W}_\mathcal {T}$ such that when we restrict it to an element $E\in \mathcal {T}$ we get $\varPi ^\nabla _E$. Finally, let $\mathcal {I}_E:\mathbb {V}_E\rightarrow \mathbb {P}_1(E)$ be the Lagrange interpolation operator at the nodes of E and the respective global Lagrange interpolation operator $\mathcal {I}_{\mathcal {T}}:\mathbb {V}_\mathcal {T}\rightarrow \mathbb {W}_\mathcal {T}$ such that $\mathcal {I}_{\mathcal {T}}|_E = \mathcal {I}_E$.

We briefly recall some properties of the projection operator $\varPi ^\nabla _E$, which will be used in the paper.

Lemma 1

For the projection operator $\varPi ^\nabla _E$ it holds the following

(2.7)

where the hidden constant does not depend on the diameter $h_E$. Furthermore,

(2.8)

Proof

The proof of (2.7) can be found in [12, Lemma 3.4]. On the other hand, (2.8) is a direct consequence of the $H^1$ orthogonality of operator $\varPi ^\nabla _E$ (2.5), i.e.

$\square$

2.2 The Discrete Problem

We assume the data to be piecewise constant on each element E of the partition. In the following, we denote the restriction of the data $(\mathcal {K}, f,c)$ over the element E as $(\mathcal {K}_E, f_E, c_E)$. We now define the discrete bilinear forms $a_\mathcal {T}$ and $m_\mathcal {T}:$ $\mathbb {V}_\mathcal {T}\times \mathbb {V}_\mathcal {T}\rightarrow \mathbb {R}$, as

$$\begin{aligned}&a_E(v,w){:}{=}\int _E \mathcal {K}_E \left( \nabla \varPi ^\nabla _E v\right) \cdot \left( \nabla \varPi ^\nabla _E w\right) ,&a_\mathcal {T}(v,w) {:}{=}\sum _{E\in \mathcal {T}}a_E(v,w); \\&m_E(v,w){:}{=}\int _E c_E \ \varPi ^\nabla _E v \;\varPi ^\nabla _E w,&m_\mathcal {T}(v,w) {:}{=}\sum _{E\in \mathcal {T}}m_E(v,w). \end{aligned}$$

We introduce the local stabilization form $s_E: \mathbb {V}_E\times \mathbb {V}_E \rightarrow \mathbb {R}$ as

$$\begin{aligned} s_E (v,w) = h_E\sum ^{\mathcal {N}_E}_{i=1} v(\varvec{x}_i) w(\varvec{x}_i), \end{aligned}$$

(2.9)

where $\{\varvec{x}_i\}_{i=1}^{\mathcal {N}_E}$ is the set of the coordinates of the vertexes of element E, and $\{v(\varvec{x}_i)\}_i^{\mathcal {N}_E}$ is the set of the degrees of freedom for the case $k = 1$. We remark that this classical choice of the stabilization term differs from the bi-dimensional classical one for the presence of the $h_E$ factor. Furthermore, we recall that the choice of the form of the stabilization term is not restrictive and the results presented can be extended to other types of the stabilization term. In [12] it is shown that, for the choice (2.9), it holds

(2.10)

with $C_s$ and $c_s$ positive constants. We define the stabilization form on the partition $\mathcal {T}$ as

$$\begin{aligned} S_\mathcal {T}(v,w){:}{=}\sum _{E\in \mathcal {T}}s_E(v-\mathcal {I}_E v, w - \mathcal {I}_E w)\qquad \forall v,w\in \mathbb {V}_\mathcal {T}, \end{aligned}$$

which yields to

(2.11)

where is the broken $H^1$-seminorm over the mesh $\mathcal {T}$. We have now all the tools to define the complete bilinear form $\mathcal {B}_\mathcal {T}:\mathbb {V}_\mathcal {T}\times \mathbb {V}_\mathcal {T}\rightarrow \mathbb {R}$, associated to problem (2.1),

$$\begin{aligned} \mathcal {B}_\mathcal {T}(v,w){:}{=}a_\mathcal {T}(v,w) + m_\mathcal {T}(v,w) + \gamma S_\mathcal {T}(v,w), \end{aligned}$$

(2.12)

where $\gamma \in \mathbb {R}$, $\gamma \ge \gamma _0 > 0$. Finally, to approximate the forcing term we define $\mathcal {F}_\mathcal {T}: \mathbb {V}_\mathcal {T}\rightarrow \mathbb {R}$ such that

$$\begin{aligned} \mathcal {F}_\mathcal {T}(v){:}{=}\sum _{E\in \mathcal {T}}\int _E f_E \varPi ^\nabla _E v, \qquad \forall \,v\, \in \mathbb {V}_\mathcal {T}. \end{aligned}$$

(2.13)

We define the discrete formulation of the problem as

$$\begin{aligned} {\left\{ \begin{array}{ll} \text {find }u_\mathcal {T}\in \mathbb {V}_\mathcal {T}& \text {such that}\\ \mathcal {B}_\mathcal {T}(u_\mathcal {T}, v)=\mathcal {F}_\mathcal {T}(v), & \forall \;v \in \mathbb {V}_\mathcal {T}, \end{array}\right. } \end{aligned}$$

(2.14)

which admits a unique solution.

The importance of $\mathbb {V}^0_\mathcal {T}$ can be already viewed in the following Lemma, whose proof can also be found in [6, Lemma 2.6].

Lemma 2

(Galerkin quasi-orthogonality) For any $v\in \mathbb {V}^0_\mathcal {T}$, it holds

$$\begin{aligned} \mathcal {B}(u - u_\mathcal {T}, v)=0, \end{aligned}$$

where u is the solution of (2.2) and $u_\mathcal {T}$ the solution of the discrete problem (2.14).

Proof

From the continuous and the discrete formulation of the variational problem we get, $\forall v \in \mathbb {V}_\mathcal {T}^0$,

$$\begin{aligned} \mathcal {B}(u- u_\mathcal {T},v)= \mathcal {B}(u, v)- \mathcal {B}(u_\mathcal {T}, v)= (f,v)- \mathcal {F}_\mathcal {T}(v)+\mathcal {B}_\mathcal {T}(u_\mathcal {T}, v)-\mathcal {B}(u_\mathcal {T},v). \end{aligned}$$

Since $v\in \mathbb {V}_\mathcal {T}^0$, $\varPi ^\nabla _\mathcal {T}v = v$, we have from (2.13) $\mathcal {F}_\mathcal {T}(v)= (f,v)$. Furthermore, we notice that from (2.5) we get

$$\begin{aligned} a_\mathcal {T}(u_\mathcal {T}, v)&= \sum _{E\in \mathcal {T}}\int _E \mathcal {K}_E\left( \nabla \varPi ^\nabla _E u_\mathcal {T}\right) \cdot \left( \nabla \varPi ^\nabla _E v \right) =\sum _{E\in \mathcal {T}}\int _E \mathcal {K}_E\left( \nabla \varPi ^\nabla _E u_\mathcal {T}\right) \cdot \left( \nabla v \right) \\&= \sum _{E\in \mathcal {T}}\int _E \mathcal {K}_E\left( \nabla u_\mathcal {T}\right) \cdot \left( \nabla v \right) = a(u_\mathcal {T},v). \end{aligned}$$

Finally, from the definition of the enhanced Virtual Element space (2.4), $\forall v \in \mathbb {V}^0_\mathcal {T}$

$$\begin{aligned} m_\mathcal {T}(u_\mathcal {T}, v) = \sum _{E\in \mathcal {T}}\int _E c_E \varPi ^\nabla _E u_\mathcal {T}\;\varPi ^\nabla _E v&=\sum _{E\in \mathcal {T}}\int _E c_E \varPi ^\nabla _E u_\mathcal {T}\; v \\ &= \sum _{E\in \mathcal {T}}\int _E c_E u_\mathcal {T}\;v = m(u_\mathcal {T},v), \end{aligned}$$

which concludes the proof. $\square$

2.3 The Preliminary Properties

We present two key properties that will be essential in the rest of the paper. The first of these properties is the scaled Poincaré inequality, which retraces the result presented in the bi-dimensional case [6, Proposition 3.1], but with crucial differences in the proof. We present previously the following Lemmas.

Lemma 3

Let E be an arbitrary element in $\mathcal {T}$, E cannot have all its vertices with the same global index $\tilde{\lambda }>0$.

Proof

We prove it by contradiction. By hypothesis, E has only hanging nodes as vertices ($\tilde{\lambda }>0$), then it is generated after a partition. We reference to Fig. 3, one of the vertices of E, say $\varvec{x}_{0}$, is the midpoint of the edge with endpoints $\varvec{x}'_{0}$ and $\varvec{x}''_{0}$. Since one of these endpoints, say $\varvec{x}''_0$, is another vertex of E, the global index of $\varvec{x}''_0$ should be $\tilde{\lambda }$. However, by Definition 1, the global index of $\varvec{x}_{0}$ results

$$\begin{aligned} \tilde{\lambda }=\lambda (\varvec{x}_{0}) = \max \{\lambda (\varvec{x}'_{0}), \lambda (\varvec{x}''_{0})\}+ 1\ge \tilde{\lambda } +1, \end{aligned}$$

which yields to a contradiction. $\square$

In particular, we notice that if we set $\varLambda =1$ in Assumption 1 all the elements of the partition contain at least a proper node.

Fig. 3

Element $E\in \mathcal {T}$ obtained after the refinement of element T via the NVB, where the edge bisected is the one with endpoints $\varvec{x}'_0$ and $\varvec{x}''_0$. $\lambda _j,\lambda _k,\lambda _\ell ,\lambda _m$ are the global indices of the vertices of the element T, whereas $\lambda _i,\lambda _j,\lambda _k,\lambda _m$ are the global indices of the vertices of E

Lemma 4

Let E be an element with only hanging nodes as vertices, that is generated by the refinement of an element T. Passing from E to its parent T, the sum of the global indices of the vertices of E strictly increases by at least one unit compared to the sum of the global indices of the vertices of T.

Proof

In this proof, we refer again to Fig. 3. We denote by $\lambda _i, \lambda _j, \lambda _k, \lambda _m$ the global indices of the vertices i, j, k, m of E respectively that are assumed to be strictly positive. Let T be the parent of E and $\lambda _\ell \ge 0$ the global index of the vertex $\ell$ of T not belonging to E. We claim that

$$\begin{aligned} \lambda _\ell < \lambda _i. \end{aligned}$$

To prove this, we observe that i is the midpoint of the edge whose global indices are $\lambda _\ell$ and $\lambda _k$ (see Fig. 3). By Definition 1, $\lambda _i= \max \{\lambda _\ell , \lambda _k\} +1$, whence if $\lambda _k\le \lambda _\ell$, then $\lambda _i =\lambda _\ell +1 > \lambda _\ell$, whereas if $\lambda _\ell <\lambda _k$, then $\lambda _i = \lambda _k+1>\lambda _\ell +1 >\lambda _\ell$. $\square$

Lemma 5

Let E be an element with only hanging nodes as vertices and let us denote the chain of ancestors of E as:

$$\begin{aligned} \mathcal {A}_N(E) =\{T_0 =E, T_1, \dots , T_N\}, \end{aligned}$$

(2.15)

where $T_N$ is the first element containing a proper node. The length of this chain is bounded by $3 (\varLambda -1)$, and it holds

$$\begin{aligned} h_{T_N} \lesssim 2^{\varLambda -1}h_{E}. \end{aligned}$$

(2.16)

Proof

Let E be an element fixed with only hanging nodes as vertices and $\mathcal {A}_N(E)$ the chain of its ancestors (2.15). The number of elements N of $\mathcal {A}_N(E)$ can be bounded as

$$\begin{aligned} N\le 3(\varLambda -1). \end{aligned}$$

(2.17)

With the intent to prove this, let us firstly consider the longest not-admissible chain (according to Lemma 3) starting from an element with vertices having the same global index $\varLambda$ to the one with the same global index 1. From Lemma 4 this chain has maximum length $4 (\varLambda -1)$. To ensure the admissiblity of the chain, according to Lemma 3, we need to exclude all the $\varLambda$ elements of the chain with the same global index. In this way, we have obtained a chain that ends with an element with no proper nodes. Finally, to close the definition of $\mathcal {A}_N(E)$ we add the element with at least one proper node, $T_N$, proving the inequality, i.e.

$$\begin{aligned} N\le 4(\varLambda -1) -\varLambda +1 = 3(\varLambda -1). \end{aligned}$$

See for instance Fig. 4, representing $\mathcal {A}_N(E)$ of the light-blue tetrahedral-shape element E in a mesh in which $\varLambda =2$ and, from (2.17), $N\le 3$. In the picture, E needs three steps to reach the element $T_3$, i.e. the one having at least one proper node as a vertex. This chain is the longest admissible one, since a further refinement of E would bring a global index equal to 3. Furthermore, we remark that the presence of hanging nodes in the elements of the chain does not influence on the length of the chain: as shown in Fig. 4, the length of the chain from $T_2$ to $T_3$ does not consider the red point.

Each time an element is refined, its volume halves, i.e. if T is the parent of E, $h^3_{T} \simeq 2 h_E^3$. Then, employing the bound (2.17), for the smallest element of the chain (2.15) and the largest one it holds

$$\begin{aligned} h^3_{T_N} \lesssim 2^{3(\varLambda -1)}h^3_{E}. \end{aligned}$$

(2.18)

Taking the cube root, we get (2.16). $\square$

Fig. 4

The four elements of the chain (2.15) from the light blue tetrahedral-shape element E to the element with a proper node $T_3$. This chain can also be seen as three successive refinements of element $T_3$. Blue lines represent the edges added at each refinement. Black points represent the vertices of the considered elements, red points denote the hanging nodes. The numbers next to the vertices are their global indices. The fixed maximum global index is $\varLambda =2$

Fig. 5

Black and red points represent the proper and the hanging nodes respectively. The light blue element has only hanging nodes as its vertices. Notice that this is obtained by a further refinement of the last element of Fig. 2

Proposition 1

(scaled Poincaré inequality in $\mathbb {V}_\mathcal {T}$) There exists a constant $C_P>0$ depending on $\varLambda$ and independent of $h_E$, such that

(2.19)

Proof

We assume $v\in \mathbb {V}_\mathcal {T}$ such that $v(\varvec{x})=0$, $\forall \varvec{x}\in \mathcal {P}$, and we divide the proof into the following steps. Let $E\in \mathcal {T}$ be an arbitrary element. If at least one of its vertices is a proper node, we can use the classical Poincaré inequality [17, 22]

(2.20)

On the other hand, if E has all hanging nodes as vertices, as shown for instance in Fig. 5, we build for E the chain $\mathcal {A}_{N}(E)$ (2.15) in Lemma 5. Let $T_N$ be the last element of the chain, from (2.16) we get

$$\begin{aligned} h^{-2}_E{\Vert }v{\Vert }^2_{0,E}\le h_E^{-2}{\Vert }v{\Vert }^2_{0, T_N}= \left( \frac{h_{T_N}}{h_E}\right) ^2 h^{-2}_{T_N} {\Vert }v{\Vert }^2_{0,T_N}\; \lesssim 2^{2(\varLambda -1)}{|}v{|}^2_{1,T_N}, \end{aligned}$$

in which the last inequality exploits the presence of a proper node in $T_N$.

We introduce now $\textbf{T}_{T_N}$ the union set of all the elements E generated by the refinements of $T_N$ (all these elements have a chain $\mathcal {A}_{N}(E)$ ending with $T_N$). Recalling the inequality (2.18), the number of the elements of $\textbf{T}_{T_N}$ is bounded by $2^{3(\varLambda -1)}$. Thus, we have

$$\begin{aligned} \sum _{E \in \textbf{T}_{T_N}} h^{-2}_{E}{\Vert }v{\Vert }^2_{0,E}\lesssim 2^{5(\varLambda -1)}{|}v{|}^2_{1,T_N}. \end{aligned}$$

Finally, summing up all the elements in $\mathcal {T}$, and denoting by $\mathcal {T}^\mathcal {H}$ the set of elements with only hanging nodes, we have

which concludes the proof. $\square$

The second preliminary property involves the interpolation error $v-\mathcal {I}_E v$ on the boundary of an element $E\in \mathcal {T}$, for any function $v \in \mathbb {V}_\mathcal {T}$. Let first L be one of the six straight sides of the element E. If L has not been refined, by definition $v|_L$ is a polynomial of degree one, and it is fully determined by the values of v at the endpoints of L, which coincides with two vertices of E. It results that $(v-\mathcal {I}_E v)|_L=0$. On the other hand, if L contains one hanging node $\tilde{\varvec{x}}$, $(\mathcal {I}_E v) (\tilde{\varvec{x}}) = \frac{1}{2} v(\varvec{x}')+ \frac{1}{2} v(\varvec{x}'')$, with $\varvec{x}'$ and $\varvec{x}''$ the endpoints of L. This suggests the definition, firstly given in [6], of the function of the hierarchical details d of v, as

$$\begin{aligned} d(v,\varvec{z})= {\left\{ \begin{array}{ll} v(\varvec{z})& \text {if }\varvec{z}\in \mathcal {P}_E,\\ v(\varvec{z}) -\frac{1}{2}\left( v(\varvec{z}')+ v(\varvec{z}'')\right) & \text {if }\varvec{z}\in \mathcal {H}_E, \end{array}\right. } \end{aligned}$$

(2.21)

where $\{\varvec{z}',\varvec{z}''\}$ is the set of endpoints of the edge having $\varvec{z}$ as midpoint. The links of this function with the global interpolation error has been clarified in [6, Lemma 7.1 and Corollary 7.2] and we here present them in the three-dimensional case, providing the necessary changes in the proofs.

Lemma 6

(global interpolation error vs hierarchical details) If $\mathcal {T}$ is a $\varLambda$-admissible partition, it holds

where the hidden constants depend on $\varLambda$.

Proof

Let $E\in \mathcal {T}$ be an element. From (2.10), and the definition of $\mathcal {I}_E$, i.e. $\mathcal {I}_E v(\varvec{x}) = v(\varvec{x})$, $\forall \varvec{x}\in \mathcal {P}_E$, we have the following equivalence of norms, $\forall v \in \mathbb {V}_E$

(2.22)

We now define the subspace of $\mathbb {V}_E$, $\tilde{\mathbb {V}}_E{:}{=}\left\{ w \in \mathbb {V}_E: w(\varvec{x})=0,\;\forall \varvec{x}\in \mathcal {P}_E\right\}$, observing that, in particular, $v-\mathcal {I}_E v \in \tilde{\mathbb {V}}_E$. On this finite-dimensional space, we can consider these two equivalent norms

$$\begin{aligned} h_E\sum _{\varvec{x}\in \mathcal {H}_E} d^2(w,\varvec{x}) \simeq h_E\sum _{\varvec{x}\in \mathcal {H}_E} {|}w(\varvec{x}){|}^2, \quad w \in \tilde{\mathbb {V}}_E. \end{aligned}$$

(2.23)

We remark that the hidden constants depend on $\varLambda$, since, from Assumption 1, both the dimension of the space $\tilde{\mathbb {V}}_E$ and the number of the possible patterns of hanging nodes on $\partial E$ (from which the function $d(w,\varvec{x})$ depends) can be bounded by $\varLambda$. Taking now $(v - \mathcal {I}_E v) \in \tilde{\mathbb {V}}_E$, and using (2.22) and (2.23) we have that

Finally, summing over all the elements of the partition, we conclude the proof. $\square$

We now move to the space $\mathbb {V}^0_\mathcal {T}$, defined in (2.6), where we set the Lagrange interpolation operator

$$\begin{aligned} \mathcal {I}_{\mathcal {T}}^0:\mathbb {V}_\mathcal {T}\rightarrow \mathbb {V}^0_\mathcal {T}, \end{aligned}$$

such that $(\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x}) = v(\varvec{x})$, $\forall \varvec{x}\in \mathcal {P}$. We recall that a function $v\in \mathbb {V}^0_\mathcal {T}$ is a polynomial of degree one in three dimensions, thus it can be fully determined by the values at the four proper nodes $\mathcal {P}_E$ at each element E. On the other hand, the value of v at a hanging node $\varvec{x}\in \mathcal {H}_E$ is determined as $(\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x}) = \frac{1}{2} (\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x}') + \frac{1}{2} (\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x}'')$, where $\{\varvec{x'}, \varvec{x''}\}$ are the endpoints of the edge for which $\varvec{x}$ is the midpoint, in a recursive way.

Lemma 7

Let $\mathcal {T}$ be a $\varLambda$-admissible partition and $\delta (v, \varvec{x}): = v(\varvec{x})- (\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x})$, it holds

Proof

We fix an element $E\in \mathcal {T}$. If $v\in \mathbb {V}_\mathcal {T}$, $\mathcal {I}_E v$ and $(\mathcal {I}_{\mathcal {T}}^0v)|_E$ differs only at the vertices of E. Remarking that $\mathcal {I}_E v(\varvec{x})=v(\varvec{x})$ $\forall \varvec{x}\in \mathcal {P}_E$ we get

We conclude the proof by summing up all the elements of the partition and recalling that $(\mathcal {I}_{\mathcal {T}}^0v)(\varvec{x}) = v(\varvec{x})$, $\forall \varvec{x}\in \mathcal {P}$, and that $\bigcup _{E\in \mathcal {T}}\mathcal {P}_E = \mathcal {N}$, since any node is a vertex for at least one element,

Proposition 2

(comparison between interpolation operators) Given a $\varLambda$-admissible partition $\mathcal {T}$ it holds

where the hidden constant is independent of $\mathcal {T}$.

Proof

We first notice that using the triangle inequality it holds $\forall v \in \mathbb {V}_\mathcal {T}$

then it is enough to bound the last term, i.e.

Thanks to Lemmas 6 and 7, the proof reduces to show that

$$\begin{aligned} h \sum _{\varvec{x}\in \mathcal {H}} \delta ^2(v, \varvec{x}) \lesssim h \sum _{\varvec{x}\in \mathcal {H}} d^2(v, \varvec{x}), \quad \quad \quad \forall v \in \mathbb {V}_\mathcal {T}, \end{aligned}$$

or, in an equivalent form, simplifying the h,

$$\begin{aligned} {\Vert }\varvec{\delta }{\Vert }_{l^2(\mathcal {H})}\lesssim {\Vert }\varvec{d}{\Vert }_{l^2(\mathcal {H})}, \end{aligned}$$

(2.24)

where $\varvec{\delta }=(\delta (\varvec{x}))_{\varvec{x}\in \mathcal {H}}{:}{=}(\delta (v, \varvec{x}))_{\varvec{x}\in \mathcal {H}}$ and $\varvec{d} = (d(\varvec{x})){:}{=}(d(v,\varvec{x}))_{\varvec{x}\in \mathcal {H}}$. Given $\varvec{x}\in \mathcal {H}$, $\{\varvec{x}',\varvec{x}'\}$ are the endpoints of the edge with $\varvec{x}$ as midpoint, and exploiting 2.21, we have that

$$\begin{aligned} \delta (\varvec{x})&= v(\varvec{x}) - \mathcal {I}_{\mathcal {T}}^0v(\varvec{x}) = v(\varvec{x})- \frac{1}{2}\mathcal {I}_{\mathcal {T}}^0(\varvec{x}') -\frac{1}{2}\mathcal {I}_{\mathcal {T}}^0(\varvec{x}'') \\ &= v(\varvec{x}) -\frac{1}{2}\left( v(\varvec{x}') + v(\varvec{x}'')\right) + \frac{1}{2}\left( v(\varvec{x}') - \mathcal {I}_{\mathcal {T}}^0v(\varvec{x}')\right) + \frac{1}{2}\left( v(\varvec{x}'') - \mathcal {I}_{\mathcal {T}}^0v(\varvec{x}'')\right) \\&=d(\varvec{x}) +\frac{1}{2}\delta (\varvec{x}') + \frac{1}{2}\delta (\varvec{x}''). \end{aligned}$$

We can define a matrix $\varvec{W}:l^2(\mathcal {H})\rightarrow l^2(\mathcal {H})$ such that $\varvec{\delta } = \varvec{W}\varvec{d}$. Then, (2.24) verifies if

$$\begin{aligned} ||\varvec{W}||_2\lesssim 1. \end{aligned}$$

We organize the global indices $\lambda \in [1,\varLambda _\mathcal {T}]$ and we define the set $\mathcal {H}_\lambda =\{\varvec{x}\in \mathcal {H}:\lambda (\varvec{x})=\lambda \}$ and $\mathcal {H}=\bigcup _{1\le \lambda \le \varLambda _\mathcal {T}} \mathcal {H}_\lambda$ so that the matrix $\varvec{W}$ can be factorized in a block-wise manner as

$$\begin{aligned} \varvec{W} = \varvec{W}_{\varLambda _\mathcal {T}} \varvec{W}_{\varLambda _\mathcal {T}-1} \cdots \varvec{W}_{1}, \end{aligned}$$

where $\varvec{W}_1$ is the identity matrix since $\delta (\varvec{x}') = \delta (\varvec{x}'')=0$, and each matrix $\varvec{W}_\lambda$ is obtained changing from the identity matrix the rows of block $\lambda$, i.e.,

We can now apply the Hölder inequality as $\Vert \varvec{W}_\lambda \Vert _2\le \left( \Vert \varvec{W}_\lambda \Vert _1 \Vert \varvec{W}_\lambda \Vert _\infty \right) ^{1/2}$, $\Vert \varvec{W}_\lambda \Vert _\infty ~\le ~ \frac{1}{2}~+\frac{1}{2}~+~1~ =2$, and $\Vert \varvec{W}_\lambda \Vert _1\le \frac{1}{2} c +1$. In particular, for the second term we exploit the fact that the number of hanging nodes (and so the constant c) depends on the quality of the mesh, i.e. on the minimum three-dimensional angle required between each face of the partition. Thus,

$$\begin{aligned} \Vert \varvec{W}\Vert _2\le \prod _{2\le \lambda \le \varLambda _\mathcal {T}}\Vert \varvec{W}_\lambda \Vert _2 \le (c+2)^{\frac{\varLambda -1}{2}}, \end{aligned}$$

which ends the proof. $\square$

We finally introduce $\tilde{\mathcal {I}}^0_\mathcal {T}: \mathbb {V}_\mathcal {T}\rightarrow \mathbb {V}^0_\mathcal {T}$ the classical Clément interpolation operator on $\mathbb {V}^0_\mathcal {T}$, which takes at each internal node the average target function on the support of the associated basis function. The following Lemma has been introduced in [6, Lemma 3.3], it is based on Proposition 1, and it does not depend on the geometric dimension.

Lemma 8

(Clément interpolation estimate) It holds

Proof

The observed difference is not due to Let $\tilde{\mathcal {I}}_\mathcal {T}: \mathbb {V}\rightarrow \mathbb {V}_\mathcal {T}$ be the Clément operator on the virtual-element space, see [19] for the details. Given $v\in \mathbb {V}$, we denote by $\tilde{v} =\tilde{\mathcal {I}}_\mathcal {T}v \in \mathbb {V}_\mathcal {T}$. We write $v - \tilde{\mathcal {I}}^0_\mathcal {T}v = (v - \tilde{v}) + (\tilde{v}- \tilde{\mathcal {I}}^0_\mathcal {T}\tilde{v}) +(\tilde{\mathcal {I}}^0_\mathcal {T}\tilde{v}- \tilde{\mathcal {I}}^0_\mathcal {T}v)$ and from the local stability of $\tilde{\mathcal {I}}^0_\mathcal {T}$ in $L^2$, we get

Furthermore, from the invariance of $\tilde{\mathcal {I}}^0_\mathcal {T}$ in $\mathbb {V}^0_\mathcal {T}$, i.e. $\tilde{\mathcal {I}}^0_\mathcal {T}(\mathcal {I}_{\mathcal {T}}^0\tilde{v}) = \mathcal {I}_{\mathcal {T}}^0\tilde{v}$, and we have $\tilde{v}- \tilde{\mathcal {I}}^0_\mathcal {T}\tilde{v}= (\tilde{v}- \mathcal {I}_{\mathcal {T}}^0\tilde{v}) + \tilde{\mathcal {I}}^0_\mathcal {T}(\mathcal {I}_{\mathcal {T}}^0\tilde{v}- \tilde{v})$. We can then conclude the proof using Proposition 1 and the stability of $\tilde{\mathcal {I}}^0_\mathcal {T}$ in $L^2$. $\square$

3 A Posteriori Error Analysis

In this section we prove the bound of the stabilization term by the residual. We essentially rely on the validity of the scaled Poincaré inequality Proposition 1, Proposition 2 and Lemma 8. We highlight that this analysis follows the same steps of the one presented in [6], as it is independent on the geometric dimension of the problem.

Following the a posteriori error analysis for the VEM in [13], we introduce the internal residual associated with element $E\in \mathcal {T}$ as

$$\begin{aligned} r_\mathcal {T}(E; v, \mathcal {D}){:}{=}f_E - c_E \varPi ^\nabla _E v, \quad \forall \, v \in \mathbb {V}_E, \end{aligned}$$

where $\mathcal {D}$ is the set of the problem data, and the jump residual over a face $F{:}{=}E_1 \cap E_2 \subset \varOmega$ of two elements $E_1$, $E_2$

$$\begin{aligned} j_\mathcal {T}(F; v, \mathcal {D}){:}{=}[[ \mathcal {K}\nabla \varPi ^\nabla _\mathcal {T}v]]_F = \left( \mathcal {K}_{E_1} \nabla \varPi ^\nabla _{E_1} v|_{E_1}\right) \cdot \varvec{n_1} + \left( \mathcal {K}_{E_2} \nabla \varPi ^\nabla _{E_2} v|_{E_2}\right) \cdot \varvec{n_2}, \end{aligned}$$

where $\varvec{n_i}$ denotes the unit normal vector to F that points outward to the element $E_i$, $i=1,2$. If the face $F\subset \partial \varOmega$, we set $j_\mathcal {T}(F; v, \mathcal {D}){:}{=}0$. We define the local estimator as

(3.1)

and the global estimator as the sum over all the elements of the partition

$$\begin{aligned} \eta ^2_\mathcal {T}(v, \mathcal {D}){:}{=}\sum _{E\in \mathcal {T}}\eta _\mathcal {T}^2(E;v,\mathcal {D}), \quad \forall v \in \mathbb {V}_\mathcal {T}. \end{aligned}$$

With this error estimator, proceeding as in [13], and [6, Proposition 4.1] the following upper and lower bounds are stated. We remark that in the following propositions, the quality of $\mathcal {T}$ affects the constants. However, we have chosen to omit this dependence in the statements, as it is not the primary focus of this paper.

Proposition 3

(upper bound) There exists a $C_{apost}>0$ depending only on $\varLambda$ and $\mathcal {D}$, such that

Proposition 4

(lower bound) There exists a constant $c_{apost}>0$ depending only on $\varLambda$, such that

The preliminary properties obtained in the previous section, allow us to obtain the bound of the stabilization term, operating exactly as [6, Proposition 4.4]. For the sake of completeness, we report the full proof, highlighting in particular the use of the scaled Poincaré inequality (Proposition 1).

Proposition 5

(bound of the stabilization term) There exists a constant $C_B$ depending only on $\varLambda$, such that

$$\begin{aligned} \gamma ^2 S_\mathcal {T}(u_\mathcal {T},u_\mathcal {T})\le C_B \eta ^2_\mathcal {T}(u_\mathcal {T}, \mathcal {D}). \end{aligned}$$

(3.2)

Proof

Let $w\in \mathbb {V}^0_\mathcal {T}$ and $e_\mathcal {T}: = u_\mathcal {T}-w$. From the definition of the bilinear form $\mathcal {B}_\mathcal {T}(\cdot , \cdot )$ (2.12), we get

$$\begin{aligned}&\gamma S_\mathcal {T}(u_\mathcal {T},u_\mathcal {T}) \\ &\quad = \gamma S_\mathcal {T}(u_\mathcal {T},e_\mathcal {T}) \\ &\quad = \mathcal {B}_\mathcal {T}(u_\mathcal {T},e_\mathcal {T}) -a_\mathcal {T}\left( u_\mathcal {T}, e_\mathcal {T}\right) - m_\mathcal {T}\left( u_\mathcal {T}, e_\mathcal {T}\right) \\&\quad = \mathcal {F}_\mathcal {T}(e_\mathcal {T}) -a_\mathcal {T}\left( u_\mathcal {T}, e_\mathcal {T}\right) - m_\mathcal {T}\left( u_\mathcal {T}, e_\mathcal {T}\right) \\&\quad = \sum _{E\in \mathcal {T}} \int _E(f_E- c_E\varPi ^\nabla _E u_\mathcal {T})\varPi ^\nabla _E e_\mathcal {T}-\sum _{E\in \mathcal {T}}\int _E \mathcal {K}_E \left( \nabla \varPi ^\nabla _E u_\mathcal {T}\right) \cdot \left( \nabla \varPi ^\nabla _E e_\mathcal {T}\right) . \end{aligned}$$

We apply (2.5) and the integration by part for the second term, noting that $\mathcal {K}_E$ is constant on E and $\varDelta \varPi ^\nabla _E u_\mathcal {T}= 0$, which yields to

$$\begin{aligned} \gamma S_\mathcal {T}(u_\mathcal {T},u_\mathcal {T})&= \sum _{E\in \mathcal {T}} \int _E(f_E- c_E\varPi ^\nabla _E u_\mathcal {T})\varPi ^\nabla _E e_\mathcal {T}-\sum _{E\in \mathcal {T}}\int _{\partial E} \varvec{n} \cdot \left( \mathcal {K}_E \nabla \varPi ^\nabla _E u_\mathcal {T}\right) e_\mathcal {T}. \end{aligned}$$

Employing now the continuity of $\varPi ^\nabla _E$ as stated in (2.7), we get

For any $\delta >0$, the latter inequality brings to

$$\begin{aligned} \gamma S_\mathcal {T}(u_\mathcal {T}, u_\mathcal {T}) \le \frac{1}{2 \delta } \eta ^2_\mathcal {T}(u_\mathcal {T}, \mathcal {D}) + \frac{\delta }{2} \varPhi _\mathcal {T}(e_\mathcal {T}) \end{aligned}$$

where

Choosing now $w = \mathcal {I}_{\mathcal {T}}^0u_\mathcal {T}$, we can apply the scaled Poincaré (1), obtaining

$$\begin{aligned} \varPhi _\mathcal {T}(u_\mathcal {T}- \mathcal {I}_{\mathcal {T}}^0u_\mathcal {T})\lesssim {|}u_\mathcal {T}- \mathcal {I}_{\mathcal {T}}^0u_\mathcal {T}{|}^2_{1,\varOmega }. \end{aligned}$$

Employing Proposition 2 and (2.11), we get $\varPhi _\mathcal {T}(u_\mathcal {T}- \mathcal {I}_{\mathcal {T}}^0u_\mathcal {T})\lesssim S_\mathcal {T}(u_\mathcal {T}, u_\mathcal {T})$, which concludes the proof for a suitable choice of parameter $\delta$. $\square$

The propositions 3, 4 and 5, bring to the stabilization-free upper and lower bounds.

Corollary 1

(stabilization-free a posteriori error estimates) There exist two constants $C_L>~0$ and $C_U>0$ such that

(3.3)

4 The Module GALERKIN

In analogy of [6], we present the extended adaptive algorithm, AVEM, consisting of an outer loop made of two modules DATA and GALERKIN. The module DATA approximates at each step the original smooth $\mathcal {D}$ of the variational problem (2.2) with piecewise constants on the elements of the partition with prescribed accuracy. The GALERKIN module takes as input a $\varLambda$-admissible partition, the approximated data and a given tolerance $\epsilon$, it iterates on the loop

$$\begin{aligned} \textsf{SOLVE} \rightarrow \textsf{ESTIMATE} \rightarrow \textsf{MARK} \rightarrow \textsf{REFINE}, \end{aligned}$$

(4.1)

producing a sequence of $\varLambda$-admissible partitions and the associated Galerkin approximate solution $u_k$, and stopping when the condition $\eta (u_k, \mathcal {D})<\epsilon$ is satisfied.

Remark 4

We here discuss only the GALERKIN module, leaving the complexity and the quasi-optimality of the AVEM algorithm in a forthcoming paper, as done in [5] for dimension two. In particular, the study of optimality would bring to get the optimal decay for the energy error with respect to the number of degrees of freedom $\texttt {Ndof}_k$ at the k-th iteration

where s is the worse decay rate between all the approximations $u_k$ and u. In dimension three, one expects to obtain $s= 1/3$, as we measure in the numerical results.

In particular, we highlight that in our analysis the modules in loop (4.1) involve:

• the Dörfler criterion [16] in the MARK step, i.e., given an input a parameter $\theta \in (0,1)$, it produces an almost minimal set $\mathcal {M}\subset \mathcal {T}$ such that

  $$\begin{aligned} \theta \; \eta ^2_\mathcal {T}(u_\mathcal {T}, \mathcal {D})\le \sum _{E\in \mathcal {M}}\eta ^2_\mathcal {T}(E;u_\mathcal {T}, \mathcal {D}); \end{aligned}$$

• the newest-vertex bisection (NVB) refinement strategy for tetrahedral-shape elements in the REFINE step [22, 23].

We also remark that in the REFINE step a further module, MAKE_ADMISSIBLE, is needed to reestablish the $\varLambda$-admissibility.

Fig. 6

Two elements E (peach) and $E'$ (light blue) sharing an edge and the node $\varvec{\hat{x}}$ (red) that violates Assumption 1 (a). If $\varvec{\hat{x}}$ does belong to the newest-edge, one refinement of E restores the $\varLambda$-admissibility (b). If $\varvec{\hat{x}}$ does not belong to the newest-edge, two and three refinements of E are needed (c) and (d), respectively. Purple, green and yellow are the new faces separating the tetrahedral-shape elements. In b–d $E'$ has not been represented for graphic purposes

Indeed, let us suppose that after a refinement a node $\varvec{\hat{x}}$ of an element $E'$ has global $\lambda (\varvec{\hat{x}}) = \varLambda +1$, violating the Assumption 1. By definition of global index of a node, $\varvec{\hat{x}}$ is a hanging node for at least an other tetrahedral-shape element E (see Fig. 6a). Two possible strategies can be adopted to restore the $\varLambda$-admissibility, which depend on the position of the node $\varvec{\hat{x}}$. If $\varvec{\hat{x}}$ belongs to the newest-edge of the element E (Fig. 6b), one refinement of element E brings to the condition $\lambda (\varvec{\hat{x}})\le \varLambda$. On the other hand, if the edge in which $\varvec{\hat{x}}$ lays cannot be immediately refined (as in Fig. 6c, d), then the NVB produces, respectively, three and four new elements.

4.1 The Convergence of GALERKIN

We discuss the reduction of the a posteriori error estimator under the refinement, to prove the convergence of GALERKIN from a mesh $\mathcal {T}$ to the refined mesh $\mathcal {T}_*$. We consider an element $E \in \mathcal {T}$ that splits into two elements $E_i \in \mathcal {T}_*$, $i=1,2$, with the introduction of a new face $F=E_1\cap E_2$. The value of a function $v\in \mathbb {V}_\mathcal {T}$ is known at all the nodes and edges of E, and then at the newest-vertex produced by the refinement. The new edges introduced do not contain any internal nodes, and then also the value of v at $\partial F$ is known. We can then deduce that v is known at all the nodes on $\partial E_1$ and $\partial E_2$. We then associate to any $v \in \mathbb {V}_\mathcal {T}$ a function $v_*\in \mathbb {V}_{\mathcal {T}_*}$ such that $v|_F = v_*|_F$.

We denote by $\eta _{\mathcal {T}_*}(E; v_*, \mathcal {D})$ the sum of the local estimators on the two new elements, i.e.

$$\begin{aligned} \begin{aligned} \eta ^2_{\mathcal {T}_*}(E; v_*, \mathcal {D}){:}{=}\sum ^2_{i=1}\eta ^2_{\mathcal {T}_*}(E_i;v_*, \mathcal {D})\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \\ =\sum ^2_{i=1}h^2_{E_i}{\Vert }r_{\mathcal {T}_*}(E; v_*, \mathcal {D}){\Vert }^2_{0,E_i} + \sum ^2_{i=1}\frac{h_{E_i}}{2}\sum _{F\in \partial \mathcal {F}_{E_i}} {\Vert }j_{\mathcal {T}_*}(F; v_*, \mathcal {D}){\Vert }^2_{0, F}, \end{aligned} \end{aligned}$$

where $h_{E_i}=\frac{1}{\root 3 \of {2}}h_E$, since the volume of E halves after the refinement. Furthermore, we remark that $\mathcal {K}_E|_{E_i} = \mathcal {K}_{E_i}$, $c_E|_{E_i}= c_{E_i}$, and $f_{E}|_{E_i}= f_{E_i}$. As a first step, we analyze the behavior of the estimator under the refinement in analogy to [6, Lemma 5.2].

Lemma 9

(local residual under refinement) Let E be an element of the partition $\mathcal {T}$ which is split into two elements $E_1$ and $E_2$ in $\mathcal {T}_*$. There exists a strictly positive constant $\mu <1$ such that $\forall v \in \mathbb {V}_\mathcal {T}$

$$\begin{aligned} \eta _{\mathcal {T}_*}(E;v_*,\mathcal {D})\le \mu \, \eta _\mathcal {T}(E;v, \mathcal {D})+ c_{er,1} S^{1/2}_{\mathcal {T}(E)}(v,v), \end{aligned}$$

where $c_{er,1}>0$ and

$$\begin{aligned} \mathcal {T}(E){:}{=}\{E'\in \mathcal {T}: \mathcal {F}_E \cap \mathcal {F}_{E'}\ne \emptyset \}, \ S_{\mathcal {T}(E)}(v,v){:}{=}\sum _{E'\in \mathcal {T}(E)}S_{E'}(v,v). \end{aligned}$$

Proof

As a first step, we write the residual on E and on its children $E_i$, $i=1,2$,

$$\begin{aligned} r_E = f-c\,\varPi ^\nabla _E v, \quad \quad \quad r_{E_i} = f-c\,\varPi ^\nabla _{E_i} v_*, \end{aligned}$$

where we have dropped the E from the piecewise constant data f and c. Note that $r_{E_i}= r_E +c \left( \varPi ^\nabla _E v -\varPi ^\nabla _{E_i} v_*\right)$ and, given $\epsilon >0$

$$\begin{aligned} \sum ^2_{i=1}h^2_{E_i}{\Vert }r_{E_i}{\Vert }_{0,E_i}^2\le \sum _{i=1}^2 h^2_{E_i}\left( (1+ \epsilon ){\Vert }r_E{\Vert }^2_{0,E_i} +\left( 1 + \frac{1}{\epsilon }\right) {\Vert }c(\varPi ^\nabla _E v - \varPi ^\nabla _{E_i} v_*){\Vert }^2_{0,E_i}\right) . \end{aligned}$$

Recalling that , we get

$$\begin{aligned} \sum ^2_{i=1}h^2_{E_i}{\Vert }r_{E_i}{\Vert }_{0,E_i}^2\le \frac{1+\epsilon }{\root 3 \of {4}}h^2_E {\Vert }r_E{\Vert }^2_{0,E}+ h^2_E \left( 1 + \frac{1}{\epsilon }\right) \frac{c^2}{\root 3 \of {4}}\sum ^2_{i=1}{\Vert }\varPi ^\nabla _E v - \varPi ^\nabla _{E_i} v_*{\Vert }^2_{0,E_i} \end{aligned}$$

(4.2)

From the triangle inequality, one has

$$\begin{aligned} {\Vert }\varPi ^\nabla _E v - \varPi ^\nabla _{E_i} v_*{\Vert }^2_{0,E_i}\le {\Vert }v-\varPi ^\nabla _E v{\Vert }^2_{0, E} + \sum _{i=1}^2 {\Vert }v_*- \varPi ^\nabla _{E_i}v_*{\Vert }^2_{0,E_i}, \end{aligned}$$

and applying the minimality property of $\varPi ^\nabla$ (2.8), we get

If we set, for instance, $\epsilon = \frac{1}{2}$, and $\mu {:}{=}\frac{1+\epsilon }{\root 3 \of {4}}=\frac{3}{2\root 3 \of {4}}<1$, the inequality (4.2) becomes

with a suitable constant $C>0$. Employing (2.11) we get

For the jump term, we have that for any $F\in \partial E$,

$$\begin{aligned} j_{\mathcal {T}_*}(F;v_*) = j_\mathcal {T}(F; v) +\left( j_{\mathcal {T}_*}\left( F; v_* \right) - j_{\mathcal {T}} (F; v)\right) . \end{aligned}$$

Fix $\epsilon >0$, we operate as we have done for the residual term, i.e.,

$$\begin{aligned} \sum ^2_{i=1} \sum _{F\in \mathcal {F}_{E_i} }h_{E_i}{\Vert }j_{\mathcal {T}_*}(F; v_*){\Vert }^2_{0, F}\le (1+\epsilon )\, I + \left( 1 + \frac{1}{\epsilon }\right) \,II \end{aligned}$$

where

$$\begin{aligned} I = \sum ^2_{i=1} \sum _{F\in \mathcal {F}_{E_i} }h_{E_i} {\Vert }j_\mathcal {T}(F;v){\Vert }^2_{0,F},\ II = \sum ^2_{i=1} \sum _{F\in \mathcal {F}_E}h_{E_i}{\Vert }j_{\mathcal {T}_*}(F;v_*)-j_\mathcal {T}(F;v){\Vert }^2_{0,F}. \end{aligned}$$

On the new face that is generated in the refinement the jump term is zero, therefore

$$\begin{aligned} I\le \frac{1}{\root 3 \of {2}}\sum _{F\in \mathcal {F}_E}h_E {\Vert }j_\mathcal {T}(F;v){\Vert }^2_{0,F}. \end{aligned}$$

For the second term, we start by denoting by $\mathcal {T}_*(E_i){:}{=}\{E' \in \mathcal {T}_*:\mathcal {F}_{E_i}\cap \mathcal {F}_{E'}\ne \emptyset \}$ and the element $E_{i,F}\in \mathcal {T}_*(E_i)$ such that, for any face $F\in \mathcal {F}_{E_i}$, $F = \partial E_i \cap \partial E_{i,F}$. Indicating by $\hat{E}_{i,F}$ the parent of $E_{i,F}$, we have

$$\begin{aligned} \begin{aligned} {\Vert }j_{\mathcal {T}_*}(F;v_*) -j_\mathcal {T}(F;v){\Vert }_{0,F} ={\Vert }[[ \mathcal {K}\nabla \left( \varPi ^\nabla _{\mathcal {T}_*} - \varPi ^\nabla _{\mathcal {T}}\right) v]] {\Vert }_{0,F}\quad \quad \quad \quad \quad \quad \quad \quad \quad \\ \le {\Vert }\mathcal {K}_E \nabla \left( \varPi ^\nabla _{E_i} -\varPi ^\nabla _E\right) v{\Vert }_{0,F} + {\Vert }\mathcal {K}_{\hat{E}_{i,F}} \nabla (\varPi ^\nabla _{E_i, F} - \varPi ^\nabla _{\hat{E}_i,F})v {\Vert }_{0,F}, \end{aligned} \end{aligned}$$

Employing the trace inequality, we get

$$\begin{aligned} II&\lesssim \sum ^2_{i=1}\sum _{E'\in \mathcal {T}_*(E_i)} {\Vert }\nabla (\varPi ^\nabla _{E'} -\varPi ^\nabla _{\hat{E}'} ) v{\Vert }^2_{0,E'}\\&\lesssim \sum ^2_{i=1}\sum _{E'\in \mathcal {T}_*(E_i)}\left( {\Vert }\nabla v - \nabla \varPi _E^\nabla v{\Vert }^2_{0,E'} + {\Vert }\nabla v - \nabla \varPi ^\nabla _{\hat{E}'}v{\Vert }^2_{0,E'}\right) \end{aligned}$$

From the minimality of the operators $\varPi ^\nabla _{E'}$ and $\varPi ^\nabla _{\hat{E}'}$, we get

$$\begin{aligned} II \lesssim \sum _{E'\in \mathcal {T}(E)}{\Vert }\nabla (v-\mathcal {I}_{E'}v){\Vert }^2_{0,E'}\lesssim \sum _{E'\in \mathcal {T}(E)} S_{E'} (v,v), \end{aligned}$$

which brings to the conclusion of the proof, with a sufficiently small $\epsilon$. $\square$

This result, in addition to the Lipschitz continuity of the estimator with respect to the function v (whose proof can be found in [6, Lemma 5.3] and does not depend on the dimension), leads to the following Proposition.

Proposition 6

(estimator reduction) There exist three positive constants independent of $\mathcal {T}$, namely $\rho <1$, $C_{er,1}$, and $C_{er,2}$, such that it holds $\forall w \in \mathbb {V}_{\mathcal {T}_*}$

$$\begin{aligned} \eta ^2_{\mathcal {T}_*}(\mathcal {T}_*;w,\mathcal {D}) \le \rho \eta ^2_\mathcal {T}(\mathcal {T}; u_\mathcal {T},\mathcal {D}) + C_{er,1} S_\mathcal {T}(u_\mathcal {T},u_\mathcal {T}) + C_{er,2} {|}u_\mathcal {T}-w{|}^2_{1,\varOmega }, \end{aligned}$$

(4.3)

where $u_\mathcal {T}\in \mathbb {V}_\mathcal {T}$ is the Galerkin solution and $\mathcal {T}_*$ the refined partition obtained by applying the step REFINE.

The second fundamental Proposition that involves the bound of the stabilization term is the quasi-orthogonality property [6, Corollary 5.8].

Proposition 7

(quasi-orthogonality of energy errors) Given $\delta \in \left( 0, \frac{1}{4}\right)$, it holds the following

(4.4)

where $u_{\mathcal {T}_*}\in \mathbb {V}_{\mathcal {T}_*}$ is the solution of problem (2.14) on the refine mesh $\mathcal {T}_*$.

These two properties guarantee the convergence of the GALERKIN module with the same steps of [6, Theorem 5.1].

Theorem 1

(contraction property of GALERKIN) Let $u_{\mathcal {T}}\in \mathbb {V}_\mathcal {T}$ the Galerkin solution, and $\mathcal {M}\subset \mathcal {T}$ the set of marked elements obtained by the procedure MARK. There exist a constant $\alpha \in (0,1)$ and $\beta >0$ such that

where $\mathcal {T}_*$ is a refinement of the partition $\mathcal {T}$ obtained by applying the step REFINE.

5 Numerical Experiments

In this section, we present several numerical tests. In all of them, we set the Dörfler parameter $\theta$ for the MARK step equal to 0.5. Furthermore, we carry on a study on the value $\varLambda$. In particular, we analyze the two extreme cases in which $\varLambda =0$ and $\varLambda =10$. We recall from Remark 1 that setting $\varLambda = 0$ we force our mesh not to admit hanging nodes, and the Adaptive Virtual Element Method algorithm coincides to the classical Adaptive Finite Element Method built on tetrahedral elements. On the other hand, $\varLambda =10$ does not represent a limit in the number of hanging nodes in the mesh for the tests reported. We label AFEM and AVEM the cases $\varLambda =0$ and $\varLambda =10$, respectively.

Fig. 7

The Fichera corner domain test solution and the refined meshes obtained

Fig. 8

Error analysis of the Fichera corner domain test

5.1 The Fichera Corner Domain

We consider the Fichera test [20, 21] consisting in a three-dimensional domain $\varOmega = (-1,1)^3\setminus (0,1)^3$, which is a natural extension of the bi-dimensional L-shape domain. We solve the Poisson problem (2.1) with $\mathcal {K}= I$, where I is the identity matrix, $c=1$, and Dirichlet boundary conditions such that the exact solution results

$$\begin{aligned} u_{ex}(\rho , \theta , \psi )= \rho ^{\alpha }, \end{aligned}$$

with $\rho$, $\theta$ and $\psi$ the polar coordinates, and $\alpha \in (0,1)$.

It is possible to take $u_{ex} \in W^k_p(\varOmega )$, $\varOmega \subset \mathbb {R}^d$, by exploiting the bound $k- \frac{d}{p} +\frac{d}{q}>0$ [11, Theorem 6.29], with $q=2$ and $p\ge 1$. In particular, there exists $p\in \left( \frac{d}{k+d/2},\frac{d}{k-\alpha }\right)$ such that the adaptive algorithm can achieve the optimal order of convergence, i.e. ${\texttt {Ndofs}}^{k/d}$, where $\texttt {Ndofs}$ indicates the number of degrees of freedom. In our case, we take $\alpha =1/2$, $d=3$, and $k=1$, and we get $p\in \left( \frac{6}{5},6\right)$. We recall that a VEM function is not known explicitly inside the elements E. Then, to compute the error, we project the solution onto $\mathbb {P}_1(E)$ via $\varPi ^\nabla _E$, and we define the error as

We set a maximum number of $\texttt {Ndofs}$ equal to 38,000.

The solution is shown in Fig. 7b, with a focus on the resulting refined meshes in the origin of the domain reported in Figs. 7c, d. As shown by the distribution of the red dots in Fig. 7d the number of hanging nodes grows in the proximity of the origin, where the solution is less regular. This is the region where the AFEM propagates more the refinements (see Fig. 7c) and the number of cells of the AFEM is greater with respect to the AVEM. In the simulation the value of $\varLambda$ does not grow up, Fig. 7a shows that the maximum global index of the partition floats between the values of one and two.

Figure 8a reports the convergence plots of $H^1_{err}$ and $\eta _\mathcal {T}$ with the optimal decay rate for both the AFEM and AVEM. Figure 8b represents a numerical confirmation of the bound found in Proposition 5, in which the stability term $S_\mathcal {T}(u_\mathcal {T},u_\mathcal {T})$ results bounded by $\eta _\mathcal {T}(u_\mathcal {T},\mathcal {D})$. The quantities $S_\mathcal {T}$ and $\eta _\mathcal {T}$ represented in all the plots are scaled with ${\Vert }\nabla u_{\text {ex}}{\Vert }_{0, \varOmega }$.

We observe that the AVEM generates a partition with the 27.33% less of three-dimensional cells $\#\mathcal {T}$ (161937 vs 222866) with respect to AFEM, reaching the same $H^1_{err}$ order of magnitude, see Fig. 8c. Finally, Fig. 8d shows that, despite the AFEM and the AVEM have approximately the same number of marked three-dimensional cells, having fixed the $\texttt {Ndofs}$, the number of refined elements in the AVEM is around 30% less than the AFEM. We remark that the observed difference stems from AFEM’s enforcement of mesh conformity, which requires extending the refinement to neighboring elements. Furthermore, reducing the number of cells results in a proportional and significant decrease in the computational cost of the global refinement process.

Fig. 9

Diffusion jumps tests: the diffusion tensor $\mathcal {K}(x,y)$ in the different portions of the domain

Fig. 10

Diffusion jumps tests: the solution and cells type distribution

5.2 Diffusion Jumps

We test our method in the presence of diffusion jumps. Let $\varOmega =[0,1]^3$ and the diffusion tensor $\mathcal {K}(x,y) = k(x,y) I$, where I is the identity matrix and

$$\begin{aligned} k(x,y) {:}{=}{\left\{ \begin{array}{ll}5 \quad \text {in } [0, 0.5] \times [0,0.5]\times [0, 0.5],\quad 1 \quad \text {in } [0, 0.5] \times [0,0.5]\times (0.5,1 ],\\ 1 \quad \text {in } (0.5, 1] \times [0,0.5]\times [0, 0.5],\quad 5 \quad \text {in } (0.5, 1] \times [0,0.5]\times (0.5, 1],\\ 1\quad \text {in } [ 0,0.5] \times (0.5,1]\times [0, 0.5],\quad 5\quad \text {in } [ 0,0.5] \times (0.5,1]\times (0.5,1],\\ 5 \quad \text {in } (0.5, 1] \times (0.5,1]\times [0, 0.5],\quad 1 \quad \text {in } (0.5, 1] \times (0.5,1]\times (0.5,1] \end{array}\right. }. \end{aligned}$$

Figure 9 illustrates the distribution of $\mathcal {K}(x, y)$ across the different regions of the domain. We assume $c=0$, $f=-1$, and impose homogeneous Dirichlet boundary conditions in Problem (2.1). After 25 refinement iterations, the solution is depicted on different sections in Fig. 10a. To better emphasize the gradient discontinuity at the diffusion jump, the solution is scaled by a factor of 10 in the sectional views at $z \in \{0.2, 0.5, 0.8, 0.98\}$.

After 25 refinement iterations, the distribution of cells obtained for the AVEM consists of $84\%$ classic tetrahedra, $12\%$ polyhedra with exactly 5 vertices, and $4\%$ polyhedra with more than 5 vertices. Figure 10b illustrates the arrangement of tetrahedra from the latter group. As expected, these tetrahedra are concentrated at the interfaces, where the diffusion term exhibits jumps.

Similar to the Fichera corner domain test, the maximum value of the global indices, $\varLambda _{\max }$, observed during the experiment is 2, which remains well below the threshold value $\varLambda = 10$, as required by Assumption 1.

Finally, we compare the results obtained using AVEM with those generated by AFEM, with the maximum number of degrees of freedom set to $10^5$. As shown in Fig. 11a and b, both methods achieve the same level of accuracy for a fixed number of degrees of freedom ${\texttt {Ndofs}}$. However, AVEM demonstrates a reduction of approximately $15\%-20\%$ in the number of three-dimensional cells.

Fig. 11

Error analysis of the diffusion jumps tests

Fig. 12

CAD applications: error analysis of the tests

5.3 CAD Applications

We finally test the AVEM in different geometries characterized by an increasing complexity. In particular, we consider three different domains coming from the real applications https://www.traceparts.com/it: Slot key (03260-10), Clearance spacer (r4781-02), and Crank Handle (gdcr10-80). We solve problem (2.1) with $\mathcal {K}= I$, where I is the identity matrix, $c=0$, and we seek the exact solution:

$$\begin{aligned} \begin{aligned}&u_{ex}(x,y,z) = -\frac{8}{x^2_l y^2_w z^2_h}\\&\times (x-x_0)(y-y_0)(z-z_0)(x-(x_0 +x_l))(y -(y_0 +y_w))(z - (z_0+z_h)), \end{aligned} \end{aligned}$$

where the two points $(x_0,y_0,z_0)$ and $(x_l,y_w,z_l)$ identifies the minimum bounding box of the considered domain. The results are summarised in Fig. 12.

Figure 12a, d, and g show the distribution of tetrahedra with aligned edges and faces obtained in the final iteration. Since the imposed solution is regular, the distribution of degrees of freedom is consequently uniform across all three geometries. Moreover, as in the previous tests, the value of $\varLambda _{\max }$ remains well below the fixed threshold of 10.

Figure 12b, e, and h show the convergence of the estimator. Refinement iterations are stopped once the threshold of $10^5$ ${\texttt {Ndofs}}$ is exceeded. All three geometries exhibit the expected theoretical slope for both AFEM and AVEM, with comparable accuracy.

Finally, the reduction of three-dimensional cells during the refinement iterations is reported in Fig. 12c, f, and i. We observe that the cell reduction strongly depends on the geometry used. Specifically, the Slot key case shows a reduction of approximately 10–15%, the Clearance spacer test shows a reduction of 20–25%, and the Crank handle demonstrates a reduction of around 25–30%. Notably, we observe a greater reduction in the number of cells as the domain geometry becomes more complex. As a noteworthy results, this leads to a significant reduction in the computational time required for the refinement process in complex geometric domains.

6 Conclusions

In this paper, we investigate the three-dimensional Adaptive Virtual Element Method (AVEM) algorithm, supported by computational experiments. We focus on tetrahedra with hanging nodes, i.e., aligned faces and edges, extending the results obtained for the two-dimensional case [6].

By introducing a new proof of the scaled Poincaré inequality, we provide an extension to the three-dimensional case of the stabilization-free a posteriori bound for the energy error and compare the proposed algorithm with the Adaptive Finite Element Method (AFEM). The AVEM shows a reduced number of mesh elements in the refinement procedure, for a fixed number of degrees of freedom and error, leading to a proportional decrease in the computational cost associated with the refinement process.