Abstract
We present here a novel approach to handling curved meshes in polytopal methods within the framework of hybrid high-order methods. The hybrid high-order method is a modern numerical scheme for the approximation of elliptic PDEs. An extension to curved meshes allows for the strong enforcement of boundary conditions on curved domains and for the capture of curved geometries that appear internally in the domain e.g. discontinuities in a diffusion coefficient. The method makes use of non-polynomial functions on the curved faces and does not require any mappings between reference elements/faces. Such an approach does not require the faces to be polynomial and has a strict upper bound on the number of degrees of freedom on a curved face for a given polynomial degree. Moreover, this approach of enriching the space of unknowns on the curved faces with non-polynomial functions should extend naturally to other polytopal methods. We show the method to be stable and consistent on curved meshes and derive optimal error estimates in \(L^2\) and energy norms. We present numerical examples of the method on a domain with curved boundary and for a diffusion problem such that the diffusion tensor is discontinuous along a curved arc.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
In recent years, there has been a trend in the computational literature toward arbitrary order polytopal methods for the approximation of partial differential equations. Such methods have a greater flexibility in the mesh requirements and can capture more intricate geometric and physical details in the domain. Being of arbitrary order, they also benefit from better convergence rates with respect to the global degrees of freedom. A short list of such methods includes discontinuous Galerkin and hybridisable discontinuous Galerkin methods [13, 18, 24], virtual element methods [1, 6, 9, 15], weak Galerkin methods [29], and polytopal finite elements [36]. However, it is well known that any approximation method on a polytopal mesh of a smooth domain (i.e. with a first-order representation of the boundary) will yield at best an order two convergence rate [35, 37]. Thus, any high-order method on curved domains requires a high-order (or exact) representation of the boundary for optimal convergence.
Developed in [23, 25], hybrid high-order (HHO) schemes are modern polytopal methods for the approximation of elliptic PDEs. A key aspect of HHO is its applicability to generic meshes with arbitrarily shaped polytopal elements. This article focuses on the extension of HHO methods to allow for curved meshes, with unknowns that capture the geometry exactly, yet still achieve optimal convergence. While the approach is presented within an HHO framework for a diffusion problem, the key ideas are more general and can be extended to related polytopal methods and to other models such as linear elasticity, or the Stokes and Navier–Stokes equations.
There has been much work on the development of discontinuous Galerkin (DG) methods on curved meshes [12, 14, 28]. We also make note of the article [31] which analyses several approaches to high-order finite element methods on curved meshes. However, for the aforementioned methods the problem is much simpler than for hybrid high-order methods due to the lack of unknowns on the mesh faces. The addition of unknowns on the mesh faces is one of the key benefits hybrid methods have over DG and other non-hybrid methods due to the strong enforcement of boundary conditions and the reduction of the global degrees of freedom via static condensation [22, Appendix B.3.2].
The article [5] proposes a virtual element method (VEM) in two dimensions for meshes possessing curved edges. For each curved edge, the authors consider the space of polynomials on a linear reference segment in \(\mathbb {R}\) and map this space onto the curved edge via a sufficiently smooth parameterisation. A similar approach is taken in the articles [4, 21]. A typical approach for hybridisable discontinuous Galerkin (HDG) methods on curved domains is to map the boundary data onto a polytopal sub-domain [cf. 19, 20]. We make note of the articles [32, 33] which also use this approach to curved boundaries.
While there has been some work on the development of hybrid high-order methods on curved meshes [8, 10, 11], the approach we take in this paper is quite different. Indeed, the article [8] approaches the issue of defining unknowns on curved faces by considering a polynomial mapping from a planar reference face onto the curved face. While this naturally requires the mesh faces to be polynomial, it also reduces the approximation order [7]. Indeed, if the mapping onto the face has effective mapping order m (see [8, Equation (5) & Remark 1]), then defining face unknowns of degree l in the reference frame will yield approximation properties of at best order \(\lfloor \frac{l}{m} \rfloor \) [8, Equation (8)]. To recover the optimal approximation order observed for straight meshes, the degree of the face polynomials in the reference frame is increased by a factor of m, yielding a very large global stencil for high-order mappings. Moreover, approximation properties in the curved faces are unknown, and the authors assume them to be true [cf. 30, botti.dispspietro:2018:assessment] in order to obtain optimal error estimates. We also make note of the conference proceedings [34] which follows the same approach using reference frame polynomials to define unknowns on curved faces for a HDG method.
An alternative approach, first considered in [10, 11], is to increase the polynomial degree of the element unknowns and weakly enforce the boundary or interface conditions without defining any unknowns on curved faces. This procedure has also been implemented for a fourth order bi-harmonic problem in a curved domain [26]. Such an approach ensures stability of the system and that optimal convergence rates are achieved. However, this method does not capture the geometry exactly and requires a finely tuned Nitsche parameter to achieve stability and consistency [cf. 30]. Moreover, without unknowns defined on curved faces, it is not clear how to design an enriched method such as that proposed in [38], whereas the method devised in this paper works seamlessly with enrichment.
In this paper, we take inspiration from the article [38] and consider unknowns on the faces to include the Neumann traces of higher-order polynomials. We note that this approach does not consider reference elements or faces but rather directly defines non-polynomial spaces on curved faces. Such an approach is therefore more closely analogous to an enriched or extended method than it is to any of the previously mentioned methods of defining unknowns on curved faces. Using this approach, we are not restricted to consider polynomial faces, but can rather take any \(C^1\) manifold. Moreover, the number of degrees of freedom on curved faces is strictly bounded above and does not grow arbitrarily large for high-order mappings. We are able to prove consistency of the scheme, and by including the space of constant functions on the faces, the method is shown to be stable. In Sect. 3, we prove optimal error estimates in energy and \(L^2\)-norm, and in Sect. 5 we present a method for the design of quadrature rules on curved elements. The paper is concluded with some numerical tests in two dimensions in Sect. 6.
1.1 Model and Assumptions on the Mesh
We take a domain \(\Omega \subset \mathbb {R}^d\), \(d\ge 2\), and consider the Dirichlet–diffusion problem: find \(u\in H^{1}_0( \Omega ) \) such that
where \({\textrm{a}}(u, v) {:}{=}({\textbf {K}}\nabla u, \nabla v)_{\Omega }\) and \(\ell (v) {:}{=}(f, v)_{\Omega }\) for some source term \(f\in L^{2}(\Omega )\) and diffusion tensor \({\textbf {K}}\) assumed to be a symmetric, piecewise constant matrix-valued function satisfying, for all \({\varvec{x}}\in \mathbb {R}^d\),
for two fixed real numbers \(0 < \underline{K}\le \overline{K}\). Here and in the following, \((\cdot , \cdot )_{X}\) is the \(L^{2}\)-inner product of scalar- or vector-valued functions on a set \(X\) for its natural measure. We shall also denote by \(\Vert \cdot \Vert _{X}\) the \(L^2\)-norm.
Let \(\mathcal {H}\subset (0, \infty )\) be a countable set of mesh sizes with a unique cluster point at \(0\). For each \(h\in \mathcal {H}\), we partition the domain \(\Omega \) into a mesh \(\mathcal {M}_{h}=(\mathcal {T}_{h}, \mathcal {F}_{h})\), where \(\mathcal {T}_{h}\) denotes the mesh elements and \(\mathcal {F}_{h}\) the mesh faces.
We suppose that the mesh elements \(\mathcal {T}_{h}\) are a disjoint set of bounded simply connected domains in \(\mathbb {R}^d\) with piece-wise \(C^1\) boundary \(\partial T\). We further suppose that \(\overline{\Omega } = \bigcup _{T\in \mathcal {T}_{h}} \overline{T}\).
We suppose that the mesh faces \(\mathcal {F}_{h}\) are a disjoint set of non-intersecting, finite, \((d-1)\)-dimensional \(C^1\) manifolds which partition the mesh skeleton: \(\bigcup _{T\in \mathcal {T}_{h}} \partial T = \bigcup _{F\in \mathcal {F}_{h}} \overline{F}\) and that for each \(F\in \mathcal {F}_{h}\) there either exists two distinct elements \(T_1,T_2\in \mathcal {T}_{h}\) such that \(F \subset \partial T_1 \cap \partial T_2\) and F is called an internal face, or there exists one element \(T\in \mathcal {T}_{h}\) such that \(F \subset \partial T \cap \partial \Omega \) and F is called a boundary face. Interior faces are collected in the set \(\mathcal {F}_{h}^i\) and boundary faces in the set \(\mathcal {F}_{h}^b\).
The parameter \(h\) is given by \(h{:}{=}\max _{T\in \mathcal {T}_{h}}h_{T}\) where, for \(X=T\in \mathcal {T}_{h}\) or \(X=F\in \mathcal {F}_{h}\), \(h_X\) denotes the diameter of \(X\). We shall also collect the faces attached to an element \(T\in \mathcal {T}_{h}\) in the set \({\mathcal {F}_{T}}:=\{F\in \mathcal {F}_{h}:F\subset T\}\). The unit normal to \(F\in {\mathcal {F}_{T}}\) pointing outside \(T\) is denoted by \({\varvec{n}}_{{T}{F}}\), and \({\varvec{n}}_{{T}}:{\partial {T}}\rightarrow \mathbb {R}^d\) is the unit normal defined by \(({\varvec{n}}_{{T}})|_F={\varvec{n}}_{{T}{F}}\) for all \(F\in {\mathcal {F}_{T}}\). We note that as each F is \(C^1\), the normal \({\varvec{n}}_{{T}{F}}\) is well defined. It is also worth noting that the normal vector \({\varvec{n}}_{{T}{F}}\) will not be constant on curved faces.
We consider the following regularity assumption on the mesh elements.
Assumption 1
(Regular mesh sequence) There exists a constant \(\varrho >0\) such that, for each \(h\in \mathcal {H}\), every \(T\in \mathcal {T}_{h}\) is connected by star-shaped sets with parameter \(\varrho \) (see [22, Definition 1.41]).
Remark 1
(Assumptions on the mesh) Assumption 1 is taken from [27, Assumption 1] however we have removed the assumption that the faces are connected by star shaped sets. We shall also note that there is no requirement that the mesh elements be polytopal or for the mesh faces to be planar.
We further require that the elements of the mesh align with the discontinuities of the diffusion tensor i.e. for each \(T\in \mathcal {T}_{h}\), \({\textbf {K}}|_T{:}{=}{\textbf {K}}_{T}\) is a constant matrix. In an analogous manner to (1.2) we define quantities \(0 < \underline{K}_{T}\le \overline{K}_{T}\) to satisfy
and we define the local diffusion anisotropy ratio \(\alpha _{T}{:}{=}\frac{\overline{K}_{T}}{\underline{K}_{T}}\).
From hereon we shall write \(f \lesssim g\) if there exists some constant C which is independent of the quantities f and g, the mesh diameter h, and of the diffusion tensor \({\textbf {K}}\), such that \(f \le C g\).
Under Assumption 1, the following continuous trace inequality holds: for all \(v\in H^{1}(T)\),
A proof of (1.4) is provided in [27]. We note that no assumption on T being polytopal is required. We also note the following inverse Sobolev inequality, a proof of which is provided for highly generic and potentially curved elements in [12, Lemma 4.23]:
where we denote by \(\mathbb {P}^{\ell }(T)\) the space of polynomials on T of degree \(\le \ell \), \(\ell \in \mathbb {N}\). Combining (1.4) and (1.5) yields the following discrete trace inequality:
2 Discrete Model
A standard hybrid high-order method on polytopal meshes defines the local discrete space as
This makes sense on polytopal meshes where F is a \((d-1)\)-dimensional hyperplane as there is no ambiguity in what is meant by \(\mathbb {P}^{k}(F)\). Indeed, on such meshes it holds that \(\mathbb {P}^{k}(F) = \mathbb {P}^{k}(\Omega )|_F = \mathbb {P}^{k}(\Omega )^d \cdot {\varvec{n}}_{{F}}\). On curved meshes, it is not so obvious what the discrete space should be.
We find that the appropriate local discrete space is that of
where we define
and \({\varvec{n}}_{{F}}\) is an arbitrary unit normal to the face F. The choice of unit normal \({\varvec{n}}_{{F}}\) does not affect the definition of \(\mathcal {P}^k(F)\). We note that, even for a curved face, there is no ambiguity in the term \(\mathbb {P}^{0}(F)\) as it represents the space of functions which are constant on the face F. We emphasise that as the unit normal \({\varvec{n}}_{{F}}\) is not constant, the space \(\mathcal {P}^k(F)\) will be non-polynomial on curved faces.
Remark 2
If F is planar (that is, a \((d-1)\)-dimensional hyperplane) then it holds that \(\mathcal {P}^k(F) = \mathbb {P}^{k}(F)\) and thus the discrete space in (2.1) coincides with the usual HHO space.
Remark 3
It suffices to take the space of unknowns on the faces as \(\mathbb {P}^{0}(F) + ({\textbf {K}}_{T_1}\nabla \mathbb {P}^{k+1}(\Omega )) \cdot {\varvec{n}}_{F} + ({\textbf {K}}_{T_2}\nabla \mathbb {P}^{k+1}(\Omega )) \cdot {\varvec{n}}_{F}\) with \(\{T_1,T_2\}=\mathcal {T}_{F}\) for stability and consistency to hold. However, we define the space as \(\mathcal {P}^k(F) = \mathbb {P}^{0}(F) + \mathbb {P}^{k}(\Omega )^d \cdot {\varvec{n}}_{{F}}\) for simpler implementation and robustness of more general models.
We shall denote by \(\mathcal {P}^{k}(\mathcal {F}_{T})\) the space
and for a given \(\underline{v}_T\in \underline{U}_{{T}}^{k}\) we shall write \(\underline{v}_T=(v_T,v_{{\mathcal {F}_{T}}})\) with \(v_T\in \mathbb {P}^{k}(T)\) and \(v_{{\mathcal {F}_{T}}}\in \mathcal {P}^k(\mathcal {F}_{T})\). The potential reconstruction \({\textrm{p}}_{{\textbf {K}},{T}}^{k+1}:\underline{U}_{{T}}^{k}\rightarrow \mathbb {P}^{k+1}(T)\) is defined as the unique solution to
We denote by \(\pi _{T}^{0, k}:L^1(T) \rightarrow \mathbb {P}^{k}(T)\) and \(\pi _{F}^{0, k}:L^1(F) \rightarrow \mathcal {P}^{k}(F)\) the \(L^2\)-orthogonal projectors onto the spaces \(\mathbb {P}^{k}(T)\) and \(\mathcal {P}^{k}(F)\), respectively. We denote by \(\pi _{{\textbf {K}}, T}^{1, k + 1}:L^1(T) \rightarrow \mathbb {P}^{k+1}(T)\) the oblique elliptic projector onto the space \(\mathbb {P}^{k+1}(T)\) satisfying
The following weighted inner-products and norms are taken from [27]. The weighted inner-product \((\cdot , \cdot )_{{\textbf {K}}, {\partial {T}}}:L^{2}({\partial {T}})\times L^{2}({\partial {T}})\rightarrow \mathbb {R}\) is defined for all \(v,w\in L^{2}({\partial {T}})\) via
For all \(r\ge 1\) and \(v\in H^{r}(T)\) the weighted \(H^{r}\)-seminorm \(|{\cdot }|_{{\textbf {K}}, H^{r}(T)}\) is defined as
Lemma 1
(Approximation properties of \(\pi _{{\textbf {K}}, T}^{1, k+1}\)) For all \(s= 1, \dots , k+1\) and \(v\in H^{k + 2}(T)\),
Proof
A proof is provided by [27, Lemma 9]. While that particular proof assumes the elements are polytopal, the proof only relies on [22, Theorem 1.50] which is provided for generic elements connected by star-shaped sets. \(\square \)
The interpolant \(\underline{I}_{{T}}^{k}:H^{1}(T)\rightarrow \underline{U}_{{T}}^{k}\) is defined by
where \(\pi _{\mathcal {F}_{T}}^{0, k}:L^{1}(T)\rightarrow \mathcal {P}^{k}(\mathcal {F}_{T})\) is defined such that \(\pi _{\mathcal {F}_{T}}^{0, k}|_F=\pi _{F}^{0, k}\) for all \(F\in \mathcal {F}_{T}\).
Lemma 2
The following commutation property holds:
Proof
It follows from the definitions of \({\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\) and \(\underline{I}_{{T}}^{k}\) that
However, as \(\nabla \cdot ({\textbf {K}}_{T}\nabla w) \in \mathbb {P}^{k}(T)\) and \(({\textbf {K}}_{T}\nabla w) \cdot {\varvec{n}}_{{T}}\in ({\textbf {K}}_{T}\nabla \mathbb {P}^{k+1}(T)) \cdot {\varvec{n}}_{{T}}\subset \mathcal {P}^{k}(\mathcal {F}_{T})\) each of the projectors \(\pi _{T}^{0, k}\) and \(\pi _{\mathcal {F}_{T}}^{0, k}\) can be removed to yield
where in the last two equalities we have integrated by parts and introduced the oblique elliptic projector using Eq. (2.4a). Taking \(w = {\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{I}_{{T}}^{k}v - \pi _{{\textbf {K}}, T}^{1, k+1}v\) we observe that
Combining with \(\int _T ({\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{I}_{{T}}^{k}v - \pi _{{\textbf {K}}, T}^{1, k+1}v) = 0\) (due to Eqs. (2.3b) and (2.4b)) we conclude that
\(\square \)
Remark 4
We note that the commutation property (2.9) is the key result required to prove consistency of the scheme and relies on the fact that \(({\textbf {K}}_{T}\nabla \mathbb {P}^{k+1}(T)) \cdot {\varvec{n}}_{{T}{F}}\subset \mathcal {P}^{k}(F)\) for each F. The additional condition that \(\mathbb {P}^{0}(F) \subset \mathcal {P}^{k}(F)\) is required for coercivity to hold.
We endow the discrete space \(\underline{U}_{{T}}^{k}\) with the seminorm
The local bilinear form \({\textrm{a}}_{{\textbf {K}},T}:\underline{U}_{{T}}^{k}\times \underline{U}_{{T}}^{k}\rightarrow \mathbb {R}\) is defined as
where \({\textrm{s}}_{{\textbf {K}},T}:\underline{U}_{{T}}^{k}\times \underline{U}_{{T}}^{k}\rightarrow \mathbb {R}\) is a local stabilisation term such that the following assumptions hold.
Assumption 2
(Local stabilisation term) The stabilisation term \({\textrm{s}}_{{\textbf {K}},T}\) is a symmetric, positive semi-definite bilinear form that satisfies:
-
1.
Stability and boundedness. For all \(\underline{v}_T\in \underline{U}_{{T}}^{k}\),
$$\begin{aligned} \alpha _{T}^{-1}\Vert \underline{v}_T\Vert _{1,{\textbf {K}},T}^2 \lesssim {\textrm{a}}_{{\textbf {K}},T}(\underline{v}_T,\underline{v}_T) \lesssim \alpha _{T}\Vert \underline{v}_T\Vert _{1,{\textbf {K}},T}^2. \end{aligned}$$(2.12) -
2.
Polynomial consistency. For all \(\underline{v}_T\in \underline{U}_{{T}}^{k}\) and \(w\in \mathbb {P}^{k+1}(T)\),
$$\begin{aligned} {\textrm{s}}_{{\textbf {K}},T}(\underline{v}_T, \underline{I}_{{T}}^{k}w) = 0. \end{aligned}$$(2.13)
An example of a stabilisation bilinear form satisfying Assumption 2 is provided in Sect. 4.
Lemma 3
(Consistency of \({\textrm{s}}_{{\textbf {K}},T}\)) Suppose \({\textrm{s}}_{{\textbf {K}},T}:\underline{U}_{{T}}^{k}\times \underline{U}_{{T}}^{k}\rightarrow \mathbb {R}\) satisfies Assumption 2. Then it holds for all \(w\in H^{k+2}(T)\) that
Proof
It follows from Eq. (2.13) that
Therefore, applying the upper bound in (2.12) and the definition (2.10) of \(\Vert \cdot \Vert _{1,{\textbf {K}},T}\) yields
Thus, we infer from Lemma 7 below that
The proof follows from the approximation properties (2.7) of \(\pi _{{\textbf {K}}, T}^{1, k+1}\). \(\square \)
2.1 Global Space and HHO Scheme
The global space of unknowns is defined as
To account for the homogeneous boundary conditions, the following subspace is also introduced,
For any \(\underline{v}_h\in \underline{U}_{h}^{k}\) we denote its restriction to an element \(T\) by \(\underline{v}_T=(v_T,v_{{\mathcal {F}_{T}}})\in \underline{U}_{{T}}^{k}\) (where, naturally, \(v_{{\mathcal {F}_{T}}}\) is defined form \((v_F)_{F\in {\mathcal {F}_{T}}}\)). We also denote by \(v_h\) the piecewise polynomial function satisfying \(v_h|_T=v_T\) for all \(T\in \mathcal {T}_{h}\). The global bilinear forms \({\textrm{a}}_{{\textbf {K}},h}:\underline{U}_{h}^{k}\times \underline{U}_{h}^{k}\rightarrow \mathbb {R}\) and \({\textrm{s}}_{{\textbf {K}},h}:\underline{U}_{h}^{k}\times \underline{U}_{h}^{k}\rightarrow \mathbb {R}\) are defined as
The HHO scheme reads: find \(\underline{u}_h\in \underline{U}_{h, 0}^{k}\) such that
where \(\ell _h:\underline{U}_{h, 0}^{k}\rightarrow \mathbb {R}\) is a linear form defined as
We define the discrete energy norm \(\Vert {\cdot }\Vert _{{\textrm{a}}, {\textbf {K}}, h}\) on \(\underline{U}_{h, 0}^{k}\) as
Lemma 4
The mapping \(\Vert {\cdot }\Vert _{{\textrm{a}}, {\textbf {K}}, h}:\underline{U}_{h, 0}^{k}\rightarrow \mathbb {R}\) defines a norm on \(\underline{U}_{h, 0}^{k}\).
Proof
As \(\Vert {\cdot }\Vert _{{\textrm{a}}, {\textbf {K}}, h}\) is clearly a seminorm we only need to prove that if \(\Vert {\underline{v}_h}\Vert _{{\textrm{a}}, {\textbf {K}}, h} = 0\) then \(\underline{v}_h= 0\). It follows from the boundedness (4.4) that
Thus, if \(\Vert {\underline{v}_h}\Vert _{{\textrm{a}}, {\textbf {K}}, h} = 0\) then it must hold that \(v_T= v_F= {\textrm{const}}\) for every \(T\in \mathcal {T}_{h}\), \(F\in \mathcal {T}_{h}\). However, we infer from the homogeneous boundary conditions that those constants must all be zero. \(\square \)
3 Error Estimates
Theorem 5
(Consistency error) The consistency error \(\mathcal {E}_h( w;\cdot ):\underline{U}_{h, 0}^{k}\rightarrow \mathbb {R}\) is the linear form defined for all \(\underline{v}_h\in \underline{U}_{h, 0}^{k}\) as
for any \(w\in H^{1}_0( \Omega ) \) such that \(\nabla \cdot ({\textbf {K}}_{T}\nabla w)\in L^{2}(\Omega )\). If such a \(w\) additionally satisfies \(w|_T\in H^{k+2}(T)\) for all \(T\in \mathcal {T}_{h}\), the consistency error satisfies
The global operators \({\textrm{p}}_{{\textbf {K}},{h}}^{k+1}:\underline{U}_{h}^{k}\rightarrow \mathbb {P}^{k+1}(\mathcal {T}_{h})\) and \(\pi _{{\textbf {K}}, h}^{1, k+1}:H^{1}(\mathcal {T}_{h})\rightarrow \mathbb {P}^{k+1}(\mathcal {T}_{h})\) are defined such that their actions restricted to an element \(T\in \mathcal {T}_{h}\) are that of \({\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\) and \(\pi _{{\textbf {K}}, T}^{1, k+1}\). The global interpolator \(\underline{I}_{h}^{k}:H^{1}(\Omega )\rightarrow \underline{U}_{h}^{k}\) is defined as \(\underline{I}_{h}^{k}v {:}{=}((\pi _{T}^{0, k}v)_{T\in \mathcal {T}_{h}},(\pi _{F}^{0, k}v)_{F\in \mathcal {F}_{h}})\).
Theorem 6
(Energy and \(L^2\) error estimates) Let \(u\in H^{1}_0(\Omega )\) be the exact solution to Eq. (1.1) and suppose the additional regularity \(u \in H^{k+2}(\mathcal {T}_{h})\). Let \(\underline{u}_h\) be the exact solution to the discrete problem (2.18). Then the following error estimates hold:
-
Energy estimate.
$$\begin{aligned} \Vert \underline{u}_h- \underline{I}_{h}^{k}u\Vert _{{\textrm{a}}, {\textbf {K}}, h} + |{\textrm{p}}_{{\textbf {K}},{h}}^{k+1}\underline{u}_h- u|_{H^{1}(\mathcal {T}_{h})} \lesssim \left( \sum _{T\in \mathcal {T}_{h}}\Big [\alpha _{T}h_{T}^{k+1}|u|_{{\textbf {K}},H^{k+2}(T)}\Big ]^2\right) _{}^\frac{1}{2}. \end{aligned}$$(3.2) -
\(L^2\) estimate. Suppose additionally that the domain \(\Omega \) is convex and \({\textbf {K}}={\textbf {I}}\) is the identity matrix, then optimal convergence in \(L^2\)-norm holds:
$$\begin{aligned} \Vert {\textrm{p}}_{{\textbf {K}},{h}}^{k+1}\underline{u}_h- u\Vert _{\Omega } \lesssim h^{k+2}|u|_{H^{k+2}(\mathcal {T}_{h})}, \end{aligned}$$(3.3)where the seminorm \(|\cdot |_{H^{s}(\mathcal {T}_{h})}\) is defined as the square-root of the sum of squares of \(|\cdot |_{H^{s}(T)}\) for any \(s\in \mathbb {N}\).
Remark 5
The \(L^2\)-error estimate is stated with identity diffusion, corresponding to a Poisson problem. However, the result follows trivially (with a hidden constant depending additionally on the anisotropy of \({\textbf {K}}\)) for any constant diffusion tensor \({\textbf {K}}\) (cf. [22], Remark3.21).
Proof of Theorems 5 and 6
The estimates (3.1) and (3.2) are provided in [27] and rely only on the design conditions stated in Assumption 2, the commutation property (2.9), the approximation properties of the elliptic projector (2.7), the consistency of \({\textrm{s}}_{{\textbf {K}},T}\) (2.14), Lemma 4, as well as standard trace and inverse estimates provided in Sect. 1.1.
To prove (3.3) we require a slightly different approach to that of [22, Theorem 2.32]. In particular, as \(\pi _{F}^{0, k}\) is not a polynomial projector [22, Equation (2.78)] does not hold in our case. However, the remainder of the proof is the same so we only have to show that
where \(z_g\) is the solution to the dual problem
As we have assumed \(\Omega \) to be convex, the following elliptic regularity holds:
Moreover, as \({\textbf {K}}={\textbf {I}}\), the following equality established in the proof of [22, Lemma 2.18] holds true:
The sum over the boundary term in (3.6) can be written as follows,
As \(\nabla \pi _{{\textbf {K}}, T}^{1, k+1}u \cdot {\varvec{n}}_{{T}{F}}\in \mathcal {P}^{k}(F)\) we may drop the projector \(\pi _{F}^{0, k}\) to write
As \(\nabla u \in \varvec{H}({{\,\textrm{div}\,}};\Omega )\), the fluxes of u are continuous across every internal face \(F\in \mathcal {F}_{h}^{{\textrm{i}}}\). Therefore, as \(\pi _{F}^{0, k} z_g = 0\) for all \(F\in \mathcal {F}_{h}^{{\textrm{b}}}\) (due to \(z_g = 0\) on \(\partial \Omega \)), it holds that
Substituting back into (3.6) yields
It follows from a Cauchy–Schwarz inequality and the consistency (2.14) that
It also follows from a Cauchy–Schwarz inequality, the continuous trace inequality (1.4) and the approximation properties (2.7) that
Thus, we need to prove that
and the proof follows from the elliptic regularity (3.5) and the bound \(\Vert g\Vert _{\Omega }\le 1\). By a continuous trace inequality and a Poincaré–Wirtinger inequality
The result holds due to the \(H^1\)-approximation properties of the \(L^2\)-projector [22, Lemma 1.43] which remain valid in curved domains.
4 Analysis of the Stabilisation
We consider here the stabilisation bilinear form defined by
however, the arguments we use to show robustness on curved meshes extend seamlessly to more general choices of stability such as those considered in [27, Section 4]. It is clear that \({\textrm{s}}_{{\textbf {K}},T}\) satisfies (2.13) so it remains to prove that (2.12) holds.
Lemma 7
It holds for all \(v\in H^{1}(T)\) that
Proof
We first note the bound
which follows from the ellipticity (1.3) of \({\textbf {K}}_{T}\). Consider, by a triangle inequality
First, we wish to bound the term \(\Vert v - \pi _{\mathcal {F}_{T}}^{0, k} v\Vert _{\partial T}\). As \(\pi _{\mathcal {F}_{T}}^{0, k}\) is the \(L^2\)-orthogonal projector on \(\mathcal {P}^{k}(\mathcal {F}_{T})\), it minimises its respective norm. Therefore, we may replace \(\pi _{\mathcal {F}_{T}}^{0, k} v\) with any element of \(\mathcal {P}^{k}(\mathcal {F}_{T})\). In particular, as \(\mathbb {P}^{0}(T)|_{\partial T} \subset \mathcal {P}^{k}(\mathcal {F}_{T})\) it holds that
It follows from the continuous trace inequality (1.4) and a Poincaré–Wirtinger inequality that
Similarly, we apply the continuous trace inequality and a Poincaré–Wirtinger inequality on the term \(h_{T}^{-1}\Vert v - \pi _{T}^{0, k} v\Vert _{\partial T}^2\) to yield
where we have applied a triangle inequality to reach the conclusion. It follows from [22, Equation (1.77)] (which invokes [22, Equation (1.74)] which does not rely on the elements being polytopal) that
Thus, we can conclude that
The proof follows by applying the ellipticity (1.3) of \({\textbf {K}}_{T}\) to yield
\(\square \)
Remark 6
We note that the inclusion \(\mathbb {P}^{0}(F)\subset \mathcal {P}^{k}(F)\) is crucial for the bound
to hold, and without this inclusion, coercivity cannot hold.
Lemma 8
(Coercivity) It holds for all \(\underline{v}_T\in \underline{U}_{{T}}^{k}\) that
Proof
It follows from the definition (2.10) of \(\Vert \cdot \Vert _{1,{\textbf {K}},T}\) that
where we have added and subtracted \(\pi _{T}^{0, k}{\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{v}_T\) to the volumetric term, and \(\pi _{T}^{0, k}{\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{v}_T\) and \(\pi _{\mathcal {F}_{T}}^{0, k}{\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{v}_T\) to the boundary term, and invoked triangle inequalities to reach the conclusion. Similar to the proof of Lemma 7, we apply the continuous trace inequality (1.4), a Poincaré–Wirtinger inequality (due to the zero mean value of \(\pi _{T}^{0, k}{\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{v}_T- v_T\)) and the ellipticity (1.3) of \({\textbf {K}}_{T}\) to yield
Therefore,
We can conclude from Lemma 7 that
which combined with the definition of \({\textrm{a}}_{{\textbf {K}},T}\) yields the result. \(\square \)
Lemma 9
(Boundedness) It holds for all \(\underline{v}_T\in \underline{U}_{{T}}^{k}\) that
Proof
Consider by a triangle inequality and Lemma 7
Similarly, by a triangle inequality and Lemma 7,
Thus, we need to prove that
It follows from the definition (2.3) of \({\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\) and an integration by parts that
where we have applied Cauchy–Schwarz inequalities on both inner-products and the discrete trace inequality (1.6). The proof follows by simplifying (4.5) by \(|{\textrm{p}}_{{\textbf {K}},{T}}^{k+1}\underline{v}_T|_{{\textbf {K}},H^{1}(T)}\) and squaring. \(\square \)
5 Integration on Curved Domains
The design of integration methods on curved domains is an active area of research. In the recent article [2], a quadrature rule for curved domains is developed by considering a decomposition into triangular or rectangular pyramids \(\mathcal {T}\) and a mapping \(\varvec{T}:[0,1]^d\rightarrow \mathcal {T}\) for each decomposition. With knowledge of the Jacobian of such a mapping, integration can be performed on the pre-image of each \(\mathcal {T}\). The article [17] develops an extension of the homogeneous integration rule developed in [16] by considering a curved triangulation of the domain and constructing a scaled boundary parameterisation on each curved triangle. Here, we also consider an extension of the homogeneous integration rule, but the approach we take is quite different. We avoid the need to split the curved domain into sub-regions and directly map the integral onto the boundary by constructing a Poincaré-type operator which inverts the divergence operator. Indeed, this operator was briefly mentioned in the appendix of [17]; however, we develop the ideas here without a sub-triangulation, and independent of dimension.
We begin with the formula developed in [16] to rewrite the integral onto the boundary of the element. This rule works by identifying a vector field
such that \(\nabla \cdot \varvec{F}= v\) for homogeneous functions v of degree q. Therefore,
We would like to extend this rule to non-homogeneous functions. We begin by searching for a vector field of the form
such that \(\nabla \cdot \varvec{F}= v\) where \(\hat{{\varvec{r}}}\) denotes the unit vector in the radial direction. We find that the unknown function g must satisfy
where we denote by \(r = |{\varvec{x}}|\). A solution is given by
Thus, we have found an inverse divergence
Therefore, an integral over the element T can be rewritten to its boundary as follows:
We note that if v is a homogeneous function of degree q (that is, \(v(t {\varvec{x}}) = t^q v({\varvec{x}})\)), then the inverse divergence formulae (5.1) and (5.3) coincide and thus so do the rules (5.2) and (5.4). In this sense, the method can be considered an extension of the homogeneous integration rule developed in [16].
If we instead consider a vector field of the form \(\varvec{F}= g\hat{{\varvec{r}}}_0\) where \(\hat{{\varvec{r}}}_0\) is the unit radial direction from a shifted origin \({\varvec{x}}_0\), we arrive at the more general formula
Therefore, we may write
This is very useful if the element contains one or more planar faces. For a vertex with coordinates \({\varvec{\nu }}\), we can set \({\varvec{x}}_0={\varvec{\nu }}\) and it holds that \(({\varvec{x}}- {\varvec{\nu }}) \cdot {\varvec{n}}_{{T}}= 0\) on any planar faces connected to the vertex \({\varvec{\nu }}\). We note that if T is not star-shaped with respect to \({\varvec{x}}_0\), then the integral \(\int _{0}^1 t^{d-1} v(t{\varvec{x}}+ (1-t){\varvec{x}}_0) \,{{\textrm{d}}} t\) will pass through points outside of T. Thus, one would require a sufficiently smooth extension of v outside of T. However, for polynomials or functions analytic over \(\Omega \) (such as an analytic source term), such an extension is trivial.
5.1 A Quadrature Rule for Curved Edges in Two Dimensions
For a given edge E, consider a parameterisation \(\gamma _E: [t_0, t_1] \rightarrow E\), \(t_0 < t_1\). Therefore, integration on curved edges is trivial:
The above integral can easily be approximated with a one-dimensional Gaussian quadrature rule. In particular, let \(w_i\), \(x_i\), \(i=1,\dots ,N\) be the weights and abscissae associated with a quadrature rule on [0, 1]. Then we can generate weights \(w_i^E\) and abscissae \({\varvec{x}}_i^E\) on the edge E as follows:
In practise, we generally store an arc length parameterisation for each edge and thus the term \(|\gamma _E'|\) is not required.
5.2 A Quadrature Rule for Elements in Two Dimensions
In two dimensions the faces are edges and thus the boundary integral in (5.6) can be evaluated on each edge \(F\in \mathcal {F}_{T}\) using the rule described in (5.7). We let \(w_i^F\) and \({\varvec{x}}_i^F\), \(i=1,\dots ,N\) be the quadrature weights and abscissae associated with an edge \(F\in \mathcal {F}_{T}\) and \(w_j\), \(x_j\), \(j=1,\dots ,M\) be the weights and abscissae associated with a quadrature rule on [0, 1]. We set \({\varvec{\nu }}\) to be the coordinate of a vertex of T connected to the highest number of straight edges in T. We then consider the quadrature rule
That is, we store weights
and abscissae
for each \(i=1,\dots ,N\), \(j=1,\dots ,M\) and on each edge \(F\in \mathcal {F}_{T}\) that is not a straight edge connected to the vertex \({\varvec{\nu }}\).
If T is polygonal, then there always exists two straight edges connected to a vertex \({\varvec{\nu }}\). Thus, the rule described by (5.8) consists of \((|\mathcal {F}_{T}| - 2) NM\) quadrature points. If we consider a Gauss-Legendre rule on each edge which is exact for polynomials of degree k, then we require to take \(N = \lceil \frac{k + 1}{2} \rceil \). However, for the inverse divergence formula (5.5) to reproduce polynomials of degree k exactly, we require to take \(M = \lceil \frac{k + 2}{2} \rceil \) due to the presence of the multiplier t. As \(({\varvec{x}}_i^F - {\varvec{\nu }})\cdot {\varvec{n}}_{{T}}({\varvec{x}}_i^F)\) is constant on polygonal T, Eq. (5.8) is exact for polynomials of degree k and consists of \((|\mathcal {F}_{T}| - 2) \lceil \frac{k + 1}{2} \rceil \lceil \frac{k + 2}{2} \rceil \) quadrature points. We note this is a slightly larger number of quadrature points than the usual \((|\mathcal {F}_{T}| - 2) \lceil \frac{k + 1}{2} \rceil ^2\) required by splitting the polygon T into \((|\mathcal {F}_{T}| - 2)\) sub-triangles (an optimal sub-triangulation) and considering a Gauss-Legendre rule on each sub-triangle. However, (5.8) avoids the complex process of generating such a sub-triangulation. To avoid these additional quadrature points, one would need to consider a Gauss-Legendre rule on each edge, but a weighted Gaussian rule with the weight function \(w(t)=t\) for the integral (5.5). This is not explored further here.
5.3 A Quadrature Rule for Elements in Three Dimensions
In three dimensions, a volumetric integral can be mapped onto the faces as follows,
Thus, given a quadrature rule for each face \(F\in \mathcal {F}_{T}\), a quadrature rule for the element T can be developed analogously to the two-dimensional case. On each face \(F\in \mathcal {F}_{T}\) let us define \(v_F({\varvec{x}}) = ({\varvec{x}}- {\varvec{x}}_0) \cdot {\varvec{n}}_{{T}{F}}\int _{0}^1 t^{2} v(t{\varvec{x}}+ (1-t){\varvec{x}}_0) \,{{\textrm{d}}} t\). Take the planar region \(\hat{F}\subset \mathbb {R}^2\) with (potentially curved) edges \(\hat{\mathcal {E}}_{\hat{F}}\) and a parameterisation \({\varvec{\gamma }}_F: \hat{F} \rightarrow F\). It holds that
where \(J(\hat{{\varvec{x}}}) = \sqrt{\det \varvec{J}^t(\hat{{\varvec{x}}}) \varvec{J}(\hat{{\varvec{x}}})}\) and \(\varvec{J}\) is the Jacobian matrix of the map \({\varvec{\gamma }}_F\). It then follows from (5.6) that
where \({\varvec{n}}_{\hat{F}\hat{E}}\) denotes the unit normal directed out of \(\hat{F}\) and toward \(\hat{E}\). Therefore, given the parameterisation \({\varvec{\gamma }}_F\) and a parameterisation of each mapped edge \(\hat{E}\in \hat{\mathcal {E}}_{\hat{F}}\), the integral (5.10) can be evaluated analogously to the 2D case (5.8).
5.3.1 A Note on Planar Faces
If the face F is planar, one can follow a procedure similar to that in [3] to rewrite the integrals on each face onto the edges \(E\in \mathcal {E}_F\). We take \({\varvec{\gamma }}_F(\hat{{\varvec{x}}}) = {\varvec{x}}_F + \varvec{E}\hat{{\varvec{x}}}\) where \({\varvec{x}}_F\) is a point in the face F and \(\varvec{E}\) is an orthonormal matrix. Then it holds that \(J(\hat{{\varvec{x}}})\equiv 1\) and
Thus, we can map the integral (5.10) back to the edges of the face F as follows,
However, as \(\varvec{E}\) is orthonormal it preserves distance and therefore it holds that \(({\varvec{\gamma }}_F^{-1}({\varvec{x}}) - \hat{{\varvec{x}}}_0) \cdot {\varvec{n}}_{\hat{F}\hat{E}} = ({\varvec{x}}- {\varvec{\gamma }}_F(\hat{{\varvec{x}}}_0)) \cdot {\varvec{n}}_{FE}\), where \({\varvec{n}}_{FE}\) denotes the unit normal directed out of F and toward E. Moreover, the mapping \({\varvec{\gamma }}_F\) is onto, so we can choose \(\hat{{\varvec{x}}}_0\) such that \({\varvec{\gamma }}_F(\hat{{\varvec{x}}}_0) = {\varvec{x}}_{F,0}\) for an arbitrary point \({\varvec{x}}_{F,0}\in F\). Therefore
Again, we may choose \({\varvec{x}}_{F,0}\) to be the vertex of the face F connected to the largest number of straight edges. The integral (5.11) is then evaluated in an identical manner as two-dimensional elements.
6 Implementation
The HHO method for curved edges is implemented using the open source C++ library PolyMesh [39]. We generate curved meshes by first considering uniform Cartesian meshes and ‘cutting’ along a curve. The integrals are computed using the quadrature rule described by (5.8) where we take the one-dimensional integration rules to be Gauss-Legendre rules of degree 30.
A basis is formed for the space \(\mathcal {P}^{k}(F)\) by first generating a spanning set by considering a canonical basis of \(\mathbb {P}^{k}(\Omega )^d\) and taking \(\mathbb {P}^{0}(F) + \mathbb {P}^{k}(\Omega )^d\cdot {\varvec{n}}_{{F}}\). The linearly dependent basis functions are removed algebraically using the FullPivLU class found in the Eigen library, with documentation available at https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html. This requires a threshold to be set which determines the point at which pivots are considered to be numerically zero. We set this value to \(10^{-15}\). We note that for sufficiently small h and large k this can result in certain linearly independent functions being removed from \(\mathcal {P}^{k}(F)\). However, as these functions are ‘close’ to being linearly dependent, the method seems unaffected by their removal. The bases of both \(\mathcal {P}^{k}(F)\) and \(\mathbb {P}^{k}(T)\) are orthonormalised via a Gram-Schmidt process.
6.1 Curved Doundary
We consider here the domain given by the rotated ellipse
where the level set \(L:\mathbb {R}^2\rightarrow \mathbb {R}\) is defined by
with \(\alpha = \frac{4}{5}\). We note the following parameterisation of \(\partial \Omega \): \(\gamma : [0, 2 \pi ) \rightarrow \partial \Omega \),
The exact solution to problem (1.1) is taken to be
with corresponding source term given by
The relative error of the scheme is measured through the following three quantities:
where the norm \(\Vert \cdot \Vert _{L^{2}(\mathcal {T}_{h})}\) is defined as the square-root of the sum of squares of \(\Vert \cdot \Vert _{L^{2}(T)}\). We note that if the mesh conforms to the domain \(\Omega \) then \(\Vert v\Vert _{L^{2}(\mathcal {T}_{h})} = \Vert v\Vert _{\Omega }\) for all \(v\in L^{2}(\Omega )\).
We consider here two sequence of meshes of the domain \(\Omega \). The curved meshes use an exact representation of the boundary, whereas the straight meshes take a piece-wise linear approximation of the boundary. The parameters of the mesh sequences are displayed in Table 1. Both sequences of meshes have the same parameters. Example curved meshes are plotted in Fig. 1 and straight meshes are plotted in Fig. 2.
In Fig. 3, we test both a curved HHO scheme and a classical HHO scheme (on straight meshes) with polynomial degrees given by \(k=1\) and \(k=3\). In both cases the curved HHO scheme on the fitted mesh observes significantly better convergence rates than the classical scheme on the straight mesh. While the scheme appears to converge optimally on curved meshes, it converges at most order 2 on straight meshes.
In Fig. 4, we test the performance of both methods as k increases on Mesh 2. While the scheme enjoys exponential convergence on the curved mesh, the classical method on the straight mesh does not converge. This is to be expected, as the straight mesh does not fit \(\Omega \) exactly and so by increasing k the scheme is converging to the solution of a different problem.
6.2 Heterogeneous Diffusion
We conclude the numerical section with a test of a diffusion problem with a piece-wise constant diffusion tensor. The HHO method requires the mesh to conform to any discontinuities in the diffusion. Thus, if the diffusion has a discontinuity along a curve, the mesh has to be curved to fit the discontinuity in the diffusion. Any polytopal mesh will require an approximation of the diffusion tensor.
We consider \(\Omega =\{(x,y):x^2+y^2<1\}\) to be the unit disc and \({\textbf {K}}\) a piece-wise constant diffusion tensor given by
We take \(R = 0.8\), \(\beta _1=10^{-6}\) and \(\beta _2=1\) which corresponds to anisotropic diffusion in the region \(r < R\), and a Poisson problem in \(r > R\). We take the source term to be \(f \equiv 1\).
Again, we consider two sequences of meshes of the domain \(\Omega \). We take both sequences to fit the domain \(\Omega \) exactly, however, the curved mesh we take to fit the discontinuity in \(\varvec{K}\) exactly and the straight mesh takes a piece-wise linear approximation of \(\varvec{K}\). The mesh data is presented in Table 2. We note that both sequences of meshes have the same parameters.
An example curved mesh and an example straight mesh is plotted in Fig. 5.
As we do not know the exact solution to this problem, we run the scheme on the finest curved mesh with \(k=7\). We denote by the discrete solution to this problem \({\textrm{p}}_{{\textbf {K}},{h}}^{k+1}\underline{u}_h= u_h^*\), which will play the role of the ‘exact’ solution. We measure the quantities
We would then like to test the performance of the scheme on coarser meshes (both curved and straight) with smaller k by investigating the behaviour of
We are less interested in the rate of convergence of these measures, but rather want to observe steady convergence, and investigate the difference between the two schemes. In Fig. 6, we plot the quantities \(E_1\) and \(E_2\) against increasing polynomial degree k where we fix the mesh to be Mesh 1. It is clear that the scheme on the straight mesh, where we consider an approximate diffusion tensor, stops converging for \(k > 1\) whereas the scheme on the curved mesh converges smoothly. In Fig. 7, we test convergence against decreasing mesh size h for polynomial degrees k. While for \(k=1\) the order of \(E_1\) and \(E_2\) are similar for both schemes (and at times smaller on the straight mesh), the convergence is much smoother on the curved mesh. For \(k=3\), the values are significantly smaller for the curved mesh, and the behaviour for the straight mesh does not differ much from the \(k=1\) case. This coincides with the previous observations that increasing k past 1 has little effect on the scheme when considering a piece-wise linear approximation of the discontinuity in the diffusion.
Finally, in Fig. 8 we show contour plots of the potential reconstructions of the discrete solutions on Mesh 1 with \(k=7\). We observe that the plot on the straight mesh seems to be distorted along the eigen vectors of \({\textbf {K}}\) (that is, \((1,1)^t\) and \((1,-1)^t\)) when compared to the plot on the curved mesh. We also plot the absolute value of the difference between the two schemes and observe that this value seems to be of greatest magnitude around the discontinuity in the diffusion tensor.
References
B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. “Equivalent projectors for virtual element methods”. In: Computers & Mathematics with Applications 66.3 (2013), pp. 376–391. https://doi.org/10.1016/j.camwa.2013.05.015.
P. Antolin, X. Wei, and A. Buffa. “Robust numerical integration on curved polyhedra based on folded decompositions”. In: Computer Methods in Applied Mechanics and Engineering 395 (2022), p. 114948. https://doi.org/10.1016/j.cma.2022.114948.
P. F. Antonietti, P. Houston, and G. Pennesi. “Fast numerical integration on polytopic meshes with applications to discontinuous Galerkin finite element methods”. In: Journal of Scientific Computing 77.3 (2018), pp. 1339–1370. https://doi.org/10.1007/s10915-018-0802-y.
E. Artioli, L. Beirão da Veiga, and M. Verani. “An adaptive curved virtual element method for the statistical homogenization of random fibre-reinforced composites”. In: Finite Elements in Analysis and Design 177 (2020), p. 103418. https://doi.org/10.1016/j.finel.2020.103418.
L. Beirão da Veiga, A. Russo, and G. Vacca. “The virtual element method with curved edges”. In: ESAIM: Mathematical Modelling and Numerical Analysis 53.2 (2019), pp. 375–404. https://doi.org/10.1051/m2an/2018052.
L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. “Basic principles of virtual element methods”. In: Mathematical Models and Methods in Applied Sciences 23.01 (2013), pp. 199– 214. https://doi.org/10.1142/S0218202512500492.
L. Botti. “Influence of reference-to-physical frame mappings on approximation properties of discontinuous piecewise polynomial spaces”. In: Journal of Scientific Computing 52.3 (2012), pp. 675–703. https://doi.org/10.1007/s10915-011-9566-3.
L. Botti and D. A. Di Pietro. “Assessment of hybrid high-order methods on curved meshes and comparison with discontinuous Galerkin methods”. In: Journal of Computational Physics 370 (2018), pp. 58–84. https://doi.org/10.1016/j.jcp.2018.05.017.
F. Brezzi, R. S. Falk, and L. D. Marini. “Basic principles of mixed virtual element methods”. In: ESAIM Math. Model. Numer. Anal. 48.4 (2014), pp. 1227–1240. https://doi.org/10.1051/m2an/2013138.
E. Burman, M. Cicuttin, G. Delay, and A. Ern. “An unfitted hybrid high-order method with cell agglomeration for elliptic interface problems”. In: SIAM Journal on Scientific Computing 43.2 (2021), A859–A882. https://doi.org/10.1137/19M1285901.
E. Burman and A. Ern. “An unfitted hybrid high-order method for elliptic interface problems”. In: SIAM Journal on Numerical Analysis 56.3 (2018), pp. 1525–1546. https://doi.org/10.1137/17M1154266.
A. Cangiani, Z. Dong, and E. H. Georgoulis. “hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements”. In: Mathematics of Computation 91.333 (2022), pp. 1–35. https://doi.org/10.1090/mcom/3667.
A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston. hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. SpringerBriefs in Mathematics. Springer, Cham, 2017, pp. viii+131. isbn: 978-3-319-67671-5; 978-3-319-67673-9.
A. Cangiani, E. H. Georgoulis, and Y. Sabawi. “Adaptive discontinuous Galerkin methods for elliptic interface problems”. In: Mathematics of Computation 87.314 (2018), pp. 2675–2707. https://doi.org/10.1090/mcom/3322.
A. Cangiani, G. Manzini, and O. J. Sutton. “Conforming and nonconforming virtual element methods for elliptic problems”. In: IMA J. Numer. Anal. 37.3 (2017), pp. 1317–1354. https://doi.org/10.1093/imanum/drw036.
E. B. Chin, J. B. Lasserre, and N Sukumar. “Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra”. In: Computational Mechanics 56.6 (2015), pp. 967–981. https://doi.org/10.1007/s00466-015-1213-7.
E. B. Chin and N Sukumar. “Scaled boundary cubature scheme for numerical integration over planar regions with affine and curved boundaries”. In: Computer Methods in Applied Mechanics and Engineering 380 (2021), p. 113796. https://doi.org/10.1016/j.cma.2021.113796.
B. Cockburn, B. Dong, J. Guzmán, M. Restelli, and R. Sacco. “A hybridizable discontinuous Galerkin method for steady-state convection-diffusion-reaction problems”. In: SIAM J. Sci. Comput. 31.5 (2009), pp. 3827–3846. https://doi.org/10.1137/080728810.https://doi.org/10.1137/080728810
B. Cockburn, W. Qiu, and M. Solano. “A priori error analysis for HDG methods using extensions from subdomains to achieve boundary conformity”. In: Mathematics of Computation 83.286 (2014), pp. 665– 699. issn: 00255718, 10886842. https://doi.org/10.1090/S0025-5718-2013-02747-0.http://www.jstor.org/stable/24488232 (visited on 12/04/2022).
B. Cockburn and M. Solano. “Solving Dirichlet Boundary-value Problems on Curved Domains by Extensions from Subdomains”. In: SIAM Journal on Scientific Computing 34.1 (2012), A497–A519. https://doi.org/10.1137/100805200.
F. Dassi, A. Fumagalli, D. Losapio, S. Scialò, A. Scotti, and G. Vacca. “The mixed virtual element method on curved edges in two dimensions”. In: Computer Methods in Applied Mechanics and Engineering 386 (2021), p. 114098. https://doi.org/10.1016/j.cma.2021.114098.
D. A. Di Pietro and J. Droniou. The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications. Vol. 19. Modeling, Simulation and Applications. https://hal.archives-ouvertes.fr/hal- 02151813: Springer International Publishing, Jan. 2020, pp. xxxi + 525. isbn: 978-3-030-37202-6. https://doi.org/10.1007/978-3-030-37203-3.
D. A. Di Pietro and A. Ern. “A hybrid high-order locking-free method for linear elasticity on general meshes”. In: Computer Methods in Applied Mechanics and Engineering 283 (2015), pp. 1–21. https://doi.org/10.1016/j.cma.2014.09.009.
D. A. Di Pietro and A. Ern. Mathematical aspects of discontinuous Galerkin methods. Vol. 69. Springer Science & Business Media, 2012. isbn: 978-3-642-22979-4. https://doi.org/10.1007/978-3-642-22980-0.
D. A. Di Pietro, A. Ern, and S. Lemaire. “An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators”. In: Computational Methods in Applied Mathematics 14.4 (2014), pp. 461–472. https://doi.org/10.1515/cmam-2014-0018.
Z. Dong and A. Ern. “Hybrid high-order method for singularly perturbed fourth-order problems on curved domains”. In: ESAIM: Mathematical Modelling and Numerical Analysis 55.6 (2021), pp. 3091–3114. https://doi.org/10.1051/m2an/2021081.
J. Droniou and L. Yemm. “Robust Hybrid High-Order Method on Polytopal Meshes with Small Faces”. In: Computational Methods in Applied Mathematics 22.1 (2022), pp. 47–71. https://doi.org/10.1515/cmam-2021-0018.
F. Hindenlang, T. Bolemann, and C.-D. Munz. “Mesh Curving Techniques for High Order Discontinuous Galerkin Simulations”. In: IDIHOM: Industrialization of High-Order Methods - A Top-Down Approach. Vol. 128. Springer, Jan. 2015, pp. 133–152. isbn: 978-3-319-12885-6. https://doi.org/10.1007/978-3-319-12886-3_8.
L. Mu, J. Wang, and X. Ye. “Weak Galerkin finite element methods on polytopal meshes”. In: International Journal of Numerical Analysis and Modeling 12.1 (2015), pp. 31–53. issn: 2617-8710. http://globalsci.org/intro/article_detail/ijnam/477.html.
F. de Prenter, C. Lehrenfeld, and A. Massing. “A note on the stability parameter in Nitsche’s method for unfitted boundary value problems”. In: Computers & Mathematics with Applications 75.12 (2018), pp. 4322– 4336. https://doi.org/10.1016/j.camwa.2018.03.032.
R. Sevilla, S. Fernández-Méndez, and A. Huerta. “Comparison of high-order curved finite elements”. In: International Journal for Numerical Methods in Engineering 87.8 (2011), pp. 719–734. doi:org/10.1002/nme.3129.
M. Solano, S. Terrana, N.-C. Nguyen, and J. Peraire. “An HDG method for dissimilar meshes”. In: IMA Journal of Numerical Analysis 42.2 (2022), pp. 1665–1699. https://doi.org/10.1093/imanum/drab059.
M. Solano and F. Vargas. “An unfitted HDG method for Oseen equations”. In: Journal of Computational and Applied Mathematics 399 (2022), p. 113721. https://doi.org/10.1016/j.cam.2021.113721.
I. de Souza Ledoino and A. F. D. Loula. “Primal HDG Methods for Elliptic Problems on Curved Meshes”. In: XLI Ibero-Latin-American Congress on Computational Methods in Engineering (CILAMCE-2020). Nov. 2020.
G. Strang and A. E. Berger. “The change in solution due to change in domain”. In: Partial differential equations. 1973, pp. 199–205. isbn: 978-0-8218-9309-8. https://doi.org/10.1090/pspum/023.
N. Sukumar and A. Tabarraei. “Conforming polygonal finite elements”. In: International Journal for Numerical Methods in Engineering 61.12 (2004), pp. 2045–2066. https://doi.org/10.1002/nme.1141.
V. Thomée. “Polygonal Domain Approximation in Dirichlet’s Problem.” In: IMA Journal of Applied Mathematics 11 (Feb. 1973), pp. 33–44. issn: 0272-4960. https://doi.org/10.1093/imamat/11.1.33.
L. Yemm. “Design and analysis of the Extended Hybrid High-Order method for the Poisson problem”. In: Advances in Computational Mathematics 48.4 (2022), pp. 1–25. https://doi.org/10.1007/s10444-022-09958-y.
L. Yemm. PolyMesh. 2022. https://github.com/liamyemm/polymesh.
Funding
Open Access funding enabled and organized by CAUL and its Member Institutions
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The author declares that they have no conflict of interest.
Additional information
Communicated by Rob Stevenson.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Yemm, L. A New Approach to Handle Curved Meshes in the Hybrid High-Order Method. Found Comput Math (2023). https://doi.org/10.1007/s10208-023-09615-w
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10208-023-09615-w