# Class Vectors: Embedding representation of Document Classes

Devendra Singh Sachan
Enixta Innovations
Hyderabad, India
<EMAIL_ADDRESS>
&Shailesh Kumar
Enixta Innovations
Hyderabad, India
<EMAIL_ADDRESS>

###### Abstract

Distributed representations of words and paragraphs as semantic embeddings in
high dimensional data are used across a number of Natural Language
Understanding tasks such as retrieval, translation, and classification. In
this work, we propose ”Class Vectors” - a framework for learning a vector per
class in the same embedding space as the word and paragraph embeddings.
Similarity between these class vectors and word vectors are used as features
to classify a document to a class.

In experiment on several sentiment analysis tasks such as Yelp reviews and
Amazon electronic product reviews, class vectors have shown better or
comparable results in classification while learning very meaningful class
embeddings.

##  1 Introduction

Text classification is one of the important tasks in natural language
processing. In text classification tasks, the objective is to categorize
documents into one or more predefined classes. This finds application in
opinion mining and sentiment analysis (e.g. detecting the polarity of reviews,
comments or tweets etc.) [2008], topic categorization ( e.g. aspect
classification of web-pages and news articles such as sports, technical etc.)
and legal document discovery etc.

In text analysis, supervised machine learning algorithms such as Naive Bayes
(NB) [1998], Logistic Regression (LR) and Support Vector Machine (SVM) [1998]
are used in text classification tasks. The bag of words [1954] approach is
commonly used for feature extraction and the features can be either binary
presence of terms or term frequency or weighted term frequency. It suffers
from data sparsity problem when the size of training data is small but it
works remarkably well when size of training data is not an issue and its
results are comparable with more complex algorithms [2012].

Using the co-occurring words information, we can learn distributed
representation of words and phrases [2005] in which each term is represented
by a dense vector in embedding space. In the skip-gram model [2013], the
objective is to maximize the prediction probability of adjacent surrounding
words given current word while global-vectors model [2014] minimizes the
difference between dot product of word vectors and the logarithm of words co-
occurrence probability.

One remarkable property of these vectors is that they learn the semantic
relationships between words i.e. in the embedding space, semantically similar
words will have higher cosine similarity. For example, the word ”gpu” will be
more similar to ”processor” than to ”camera”. To use these word vectors in
classification tasks, Le et al. [2014] proposed the Paragraph Vectors
approach, in which they learn the vectors representation for documents by
stochastic gradient descent and the gradient is computed by backpropagation of
the error from the word vectors. The document vectors and the word vectors are
learned jointly. Kim [2014] demonstrated the application of Convolutional
Neural Networks in sentence classification tasks using the pre-trained word
embeddings.

Taking inspiration from the paragraph vectors approach, we propose class
vectors method in which we learn a vector representation for each class. These
class vectors are semantically similar to vectors of those words which
characterizes the class and also give competitive results in document
classification tasks.

##  2 Model

We use skip-gram model [2013] to learn these vectors. In the skip-gram
approach, we learn the parameters of model to maximize the prediction
probability of the cooccurence of words. Let the words in the corpus be
represented as $w_{1}$, $w_{2}$, $w_{3}$, .., $w_{n}$. The objective function
is defined as,

$L=\sum_{i=1}^{N_{s}}\sum_{c\epsilon[-w,w],c\neq 0}\log p(w_{i+c}/w_{i})$ (1)

where $N_{s}$ is the number of words in the sentence(corpus) and $L$ denotes
the likelihood of the observed data. $w_{t}$ denotes the current word, while
$w_{t+c}$ is the context word within a window of size $w$. The prediction
probability $p(w_{i+c}/w_{i})$ is calculated using the softmax classifier as
below,

$p(w_{i+c}/w_{i})=\frac{\exp\left({v_{w_{i}}^{\intercal}v^{\prime}_{w_{i+c}}}\right)}{\sum_{w=1}^{T}\exp\left(v_{w_{i}}^{\intercal}{v^{\prime}_{w}}\right)}$
(2)

$T$ is number of unique words selected from corpus in the dictionary,
$v_{w_{i}}$ is the vectors representation of the current word from inner layer
of neural network while ${v^{\prime}_{w}}$ is the vector representation of the
context word from the outer layer of the neural network. In practice, since
the size of dictionary can be quite large, the cost of computing the
denominator in the above equation can be very expensive and thus gradient
update step becomes impractical.

Morin et al. [2005] proposed Hierarchical Softmax to speed up the training.
They construct a binary Huffman tree to compute the probability distribution
which gives logarithmic speedup $\log_{2}(T)$. Mikolov et al. [2013] proposed
negative sampling which approximates $\log{p(w_{i+c}/w_{i})}$ as,

$\log\sigma({v_{w_{i}}^{\intercal}v^{\prime}_{w_{i+c}}})+\sum_{j=1}^{k}\mathbb{E}_{w_{j}\sim
P_{n}(w)}\left(\log\sigma(-{v_{w_{i}}}^{\intercal}v^{\prime}_{w_{j}})\right)$
(3)

$\sigma(x)$ is the sigmoid function, the word ${w_{j}}$ is sampled from
probability distribution over words ${P_{n}(w)}$. The word vectors are updated
by maximizing the likelihood $L$ using stochastic gradient ascent.

Our model, shown in Figure 1, learns a vector representation for each of the
classes along with word vectors in the same embedding space. We represent each
class vector by its id (class_id). Each class id co-occurs with every sentence
and thus with every word in that class. Basically, each class id has a window
length of the number of words in that class. We call them as Class Vectors
(CV). Following eq1 new objective function becomes,

$\sum_{i=1}^{N_{s}}\sum_{c\epsilon[-w,w],c\neq 0}\log
p(w_{i+c}/w_{i})+\lambda\sum_{j=1}^{N_{c}}\sum_{i=1}^{N_{j}}\log
p(w_{i}/c_{j})$ (4)

$N_{c}$ is the number of classes, $N_{j}$ is the number of words in
$class_{j}$, $c_{j}$ is the class id of the $class_{j}$. We use skipgram
method to learn both the word vectors and class vectors.

class_idsen1sen2sen3 Figure 1: Class Vectors model. While training each class
vector is represented by an id. Every word in the sentence of that class co-
occurs with its class vector. Class vectors and words vectors are jointly
trained using skip-gram approach.

###  2.1 Class Vector based scoring

Converting class vector to word similarity to probabilistic score using
softmax function

$s(w_{j}/c_{i})=\frac{\exp\left({v_{c_{i}}^{\intercal}v_{w_{j}}}\right)}{\sum_{w=1}^{T}\exp\left(v_{c_{i}}^{\intercal}v_{w_{j}}\right)}$
(5)

$v_{c_{i}}$ and $v_{w_{j}}$ are the inner un-normalised $ith$ class vector and
$jth$ word vector respectively. To predict the class of test data, we use
different ways as described below

  * •

We do summation of probability score for all the words in sentence for each
class and predict the class with the maximum score. (CV Score)

$\underset{i=1,..,C}{\operatorname{arg\,max}}{\sum_{j=1}^{N_{s}}\log(s(w_{j}/c_{i}))}$
(6)

  * •

We take the difference of the probability score of the class vectors and use
them as features in the bag of words model followed by Logistic Regression
classifier. For example, in the case of sentiment analysis, the two class are
positive and negative. So, the expression becomes, (CV-LR)

$f(w)=\log(s(w/c_{pos}))-\log(s(w/c_{neg}))$ (7)

$w$ is the vector of the words in vocabulary.

  * •

We compute the similarity between class vectors and word vectors after
normalizing them by their l2-norm and using the difference between the
similarity score as features in bag of words model. (norm CV-LR)

$f(w)={v_{c_{pos}}^{\intercal}v_{w}}-{v_{c_{neg}}^{\intercal}v_{w}}$ (8)

###  2.2 Feature Selection

Important features in the corpus can be selected by information theoretic
criteria such as conditional entropy and mutual information. We assume the
entropy of the class to be maximum i.e. $H(C)=1$ irrespective of the number of
documents in each class. Realized information of class given a feature $w_{i}$
is defined as,

$I(C;w=w_{i})=H(C)-H(C/w=w_{i})$ (9)

where conditional entropy of class, $H(C/w_{i})$ is,

$H(C/w=w_{i})=-\sum_{c_{i}}^{N_{c}}p(c_{i}/w_{i})\log_{2}p(c_{i}/w_{i})$ (10)

$p(c/w_{i})=\frac{\exp\left({v_{c_{i}}^{\intercal}v_{w_{i}}}\right)}{\sum_{c_{i}}^{N_{c}}\exp\left(v_{c_{i}}^{\intercal}v_{w_{i}}\right)}$
(11)

We calculate expected information $I(C;w)$ also called mutual information for
each word as,

$I(C;w)=H(C)-\sum_{w}p(w)H(C/w)$ (12)

$p(w)$ is calculated from the document frequency of word. We plot expected
information vs realized information to see the important features in the
dataset.

##  3 Dataset description

We did experiments on Amazon Electronic Reviews corpus and Yelp Restaurant
Reviews. The task is to do sentiment classification among 2 classes ( i.e.
each review can belong to either positive class or negative class ) .

  * •

Amazon Electronic Product reviews 111http://riejohnson.com/cnn_data.html \-
This dataset is a part of large Amazon reviews dataset McAuley et
al.,[2013]222http://snap.stanford.edu/data/web-Amazon.html. This dataset
[2015] contains training set of 392K reviews split into various various sizes
and a test set of 25K reviews. We pre-process the data by converting the text
to lowercase and removing some punctuation characters.

  * •

Yelp Reviews corpus 333https://www.kaggle.com/c/yelp-recruiting/data \- This
reviews dataset was provided by Yelp as a part of Kaggle competition. Each
review contains star rating from 1 to 5. Following the generation of above
IMDB Movie Reviews and Amazon Electronic Product Reviews data we considered
ratings 1 and 2 as negative class and 4 and 5 as positive class. We separated
the files into ratings and do pre-processing of the corpus. 444We use the code
available at https://github.com/TaddyLab/deepir/blob/master/code/parseyelp.py
[2015] In this way, we obtain around 193K reviews for training and around 20K
reviews for testing.

Dataset | Pos Train | Neg Train | Test Set
---|---|---|---
Amazon | 196000 | 196000 | 25000
Yelp | 154506 | 38172 | 19931
Table 1:  Dataset summary. Pos Train: Number of training examples in positive
class. Neg Train: Number of training examples in negative class. Test Set:
Number of reviews in Test Set

##  4 Experiments

We do phrase identification in the data by two sequential iterations using the
approach as described in Kumar et al. [2014]. We select the top important
phrases according to their frequency and coherence and annotate the corpus
with phrases. To do experiments and train the models, we consider those words
whose frequency is greater than 5. We use this common setup for all the
experiments.

We did experiments with following methods. In the bag of words(bow) approach
in which we annotate the corpus with phrases as mentioned earlier. We report
the best results among the bag of words in table 2. In the bag of words
method, we extract the features by using

  1. 1.

presence/absence of words (binary)

  2. 2.

term frequency of the words (tf)

  3. 3.

inverse document frequency of words (idf)

  4. 4.

product of term frequency and inverse document frequency of words (tf-idf)

We also evaluate some of the recent state of the art methods for text
classification on the above datasets

  1. 1.

naive bayes features in bag of words followed by Logistic Regression (NB-LR)
[2012]

  2. 2.

inversion of distributed language representation (W2V inversion) [2015] 555We
use the code available at https://github.com/TaddyLab/deepir which builds on
top of gensim toolkit [2010]

  3. 3.

Convolutional Neural Networks for text categorization (CNN) [2015]

  4. 4.

Paragraph Vectors - Distributed Bag of Words Model (PV-DBOW) [2014]

Class Vector method based scoring and feature extraction. We extend the open-
source code https://code.google.com/p/word2vec/ to implement the class vectors
approach. We learn the class vectors and word embeddings using these
hyperparameter settings (window=10, negative=5, min_count=5, sample=1e-3,
hs=1, iterations=40, $\lambda$=1). For prediction, we experiment with the
three approaches as mentioned above. (2.1)

After the features are extracted we train Logistic Regression classifier in
scikit-learn [2011] to compute the results. 666http://scikit-
learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
Results of our model and other models are listed in table 2.

Model | Amazon | Yelp
---|---|---
bow binary | 91.29 | 92.48
bow tf | 90.49 | 91.45
bow idf | 92.00 | 93.98
bow tf-idf | 91.76 | 93.46
Naive Bayes | 86.25 | 89.77
NB-LR | 91.49 | 94.68
W2V inversion | – | 93.3
CNN | 92.86 | –
PV-DBOW | 90.07 | 92.86
CV Score | 84.06 | 87.85
norm CV-LR | 91.58 | 94.91
CV-LR | 91.70 | 94.83
Table 2:  Comparison of accuracy scores for different algorithms  Figure 2:
Expected information vs Realized information using normalized vectors for
$1500$ most frequent words in Yelp Reviews Corpus

##  5 Results and Discussion

  1. 1.

We found that annotating the corpus by phrases is important to give better
results. For example, the accuracy of PV-DBOW method on Yelp Reviews increased
from 89.67% (without phrases) to 92.86% (with phrases) which is more than 3%
increase in accuracy .

  2. 2.

Class vectors have high cosine similarity with words which discriminate
between classes. For example, when trained on Yelp reviews, positive class
vector was similar to words like ”very_very_good”, ”fantastic” while negative
class vector was similar to words like ”awful”, ”terrible” etc. More results
can be seen in Table 3 and Table 4.

  3. 3.

In Figure 2, we see that class informative words have greater values of both
expected information and realized information. One advantage of class vectors
based feature selection method over document frequency based method is that
low frequency words can have high mutual information value.

  4. 4.

On Yelp reviews dataset, we find that the class vectors based approach (CV-LR
and norm CV-LR) performs much better than normalized term frequency (tf), tf-
idf weighted bag of words, paragraph vectors and W2V inversion and it achieves
competitive results in sentiment classification.

  5. 5.

On Amazon reviews dataset, bow idf performs surprisingly well and outperforms
all other methods except CNN based approach.

  6. 6.

Shuffling the corpus is important to learn high quality class vectors. When
learning the class vectors using only the data of that class, we find that
class vectors lose their discriminating power. So, it is important to jointly
learn the model using full dataset.

Amazon Electronic Product Reviews
---
Top Similar Words to
Pos class vector | Neg class vector
very_pleased | unfortunately
product_works_great | very_disappointed
awesome | piece_of_crap
more_than_i_expected | piece_of_garbage
very_satisfied | hunk_of_junk
great_buy | awful_service
so_good | even_worse
great_product | sadly
very_happy | worthless
am_very_pleased | terrible
a_great_value | useless
it_works_great | never_worked
works_like_a_charm | horrible
great_purchase | terrible_product
fantastic | wasted_my_money
Table 3:  Top 15 similar words to the positive class vector and negative class
vector.  Yelp Restaurant Reviews
---
Top Similar Words to
Pos class vector | Neg class vector
very_very_good | awful
fantastic | terrible
awesome | horrible
amaz | fine_but
very_yummy | food_wa_cold
great_too | awful_service
excellent | horrib
real_good | not_very_good
spot_on | pathetic
great | tastele
food_wa_fantastic | mediocre_at_best
very_good_too | unacceptable
love_thi_place | disgust
food_wa_awesome | food_wa_bland
very_good | crappy_service
Table 4:  Top 15 similar words to the positive class vector and negative class
vector.

##  6 Conclusion and Future Work

We learned the class vectors and used its similarity with words in vocabulary
as features effectively in text categorization tasks.

There is a lot of scope for further work and research such as using pre-
trained word vectors to compute the class vectors. This will help in cases
when training data is small. In order to use more than 1-gram as features we
need approaches to compute the embeddings of n-grams from the composition of
its uni-grams. Recursive Neural Networks of Socher et al [2013] can be applied
in these case. We can also work on generative models of class based on word
embeddings and its application in text clustering and text classification.

##  7 Acknowledgements

We thank the anonymous reviewers for great discussion and feedback.

## References

  * [1954]  Zellig Harris.  1954\.  Distributional structure.  Word, 10(23):146–162.
  * [1998]  Thorsten Joachims.  1998\.  Text categorization with suport vector machines: Learning with many relevant features.  In Proceedings of the 10th European Conference on Machine Learning, ECML ’98, pages 137–142, London, UK, UK. Springer-Verlag.
  * [2015]  Rie Johnson and Tong Zhang.  2015\.  Effective use of word order for text categorization with convolutional neural networks.  In Proceedings of the 2015 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, pages 103–112, Denver, Colorado, May–June. Association for Computational Linguistics.
  * [2014]  Yoon Kim.  2014\.  Convolutional neural networks for sentence classification.  CoRR, abs/1408.5882.
  * [2014]  S. Kumar.  2014\.  Phrase identification in a sequence of words, November 18.  US Patent 8,892,422.
  * [2014]  Quoc V. Le and Tomas Mikolov.  2014\.  Distributed representations of sentences and documents.  In Proceedings of the 31 st International Conference on Machine Learning.
  * [2013]  J. J. McAuley and J. Leskovec.  2013\.  Hidden factors and hidden topics: understanding rating dimensions with review text.  In Recommender Systems.
  * [1998]  Andrew McCallum and Kamal Nigam.  1998\.  A comparison of event models for naive bayes text classification.
  * [2013]  Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S. Corrado, and Jeff Dean.  2013\.  Distributed representations of words and phrases and their compositionality.  In Advances in Neural Information Processing Systems, pages 3111–3119.
  * [2005]  Frederic Morin and Yoshua Bengio.  2005\.  Hierarchical probabilistic neural network language model.  In Proceedings of the International Workshop on Artificial Intelligence and Statistics, pages 246–252.
  * [2008]  Bo Pang and Lillian Lee.  2008\.  Opinion Mining and Sentiment Analysis.  Foundations and Trends in Information Retrieval, 1-2:1–135.
  * [2011]  F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830.
  * [2014]  Jeffrey Pennington, Richard Socher, and Christopher Manning.  2014\.  Glove: Global vectors for word representation.  In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 1532–1543, Doha, Qatar, October. Association for Computational Linguistics.
  * [2010]  Radim Řehůřek and Petr Sojka.  2010\.  Software Framework for Topic Modelling with Large Corpora.  In Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks, pages 45–50, Valletta, Malta, May. ELRA.  http://is.muni.cz/publication/884893/en.
  * [2013]  Richard Socher, Alex Perelygin, Jean Y. Wu, Jason Chuang, Christopher D. Manning, Andrew Y. Ng, and Christopher Potts.  2013\.  Recursive deep models for semantic compositionality over a sentiment treebank.  In Proceedings of the conference on empirical methods in natural language processing (EMNLP), volume 1631, page 1642.
  * [2015]  Matt Taddy.  2015\.  Document classification by inversion of distributed language representations.  In Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics.
  * [2012]  Sida I. Wang and Christopher D. Manning.  2012\.  Baselines and bigrams: Simple, good sentiment and topic classification.  In Proceedings of the ACL, pages 90–94.


Design and analysis of the Extended Hybrid High-Order method for the Poisson problem
Liam Yemm
School of Mathematics, Monash University, Melbourne, Australia<EMAIL_ADDRESS>
We propose an Extended Hybrid High-Order scheme for the Poisson problem with solution possessing weak singularities. Some general assumptions are stated on the nature of this singularity and the remaining part of the solution. The method is formulated by enriching the local polynomial spaces with appropriate singular functions. Via a detailed error analysis, the method is shown to converge optimally in both discrete and continuous energy norms. Some tests are conducted in two dimensions for singularities arising from irregular geometries in the domain. The numerical simulations illustrate the established error estimates, and show the method to be a significant improvement over a standard Hybrid High-Order method.

Key words: Hybrid High-Order methods, enriched scheme, error analysis, singular solution, polytopal meshes.

MSC2010: 65N12, 65N15, 65N30.

§ INTRODUCTION

This century has seen a growing interest in so called polytopal methods for the discretisation of elliptic problems. These methods extend the classical Finite Element method to generic polytopal grids. A short list of such methods includes Discontinuous Galerkin and Hybridizable Discontinuous Galerkin methods [20, 12, 15], Multi-Point Flux Approximation finite volume methods [1], Hybrid Mimetic Mixed methods [25], Virtual Element methods [6, 2, 11, 13], weak Galerkin methods [31], and polytopal Finite Elements [32]. The Hybrid High-Order (HHO) method is a recent addition to these techniques, and also benefits from being of arbitrary order and dimension independent. Additionally, HHO methods allow for static condensation of the system matrix, and depend on polynomial reconstructions which account for the physics of the problem. For a comprehensive discussion of the method and its applications we refer the interested reader to the monograph [18].

Hybrid High-Order methods have been used to model diffusion and diffusion-advection-reaction equations [22, 19], elasticity problems [21, 10], Leray Lions and p-Laplace equations [16, 18], and the Stokes and Navier-Stokes equations [23, 9, 18]. The method has also been shown to be robust on highly irregular grids consisting of elements with arbitrarily many small faces [26] and on skewed meshes such that each element possesses a linear map to an an isotropic element [24]. All of these schemes, however, rely on error estimates which assume certain regularity of the exact solution. This is typical of high-order approximations which require high-order regularity to obtain optimal approximation rates. On non-smooth domains (such as regions with non-convex corners or those possessing cracks) it is expected that the exact solution to elliptic problems will contain weak singularities [27]. This lack of regularity is well documented in the finite element literature and is typically overcome through enriched approximations based on the partition of unity method [29, 5]. The extended Finite Element method [7, 30] is one such method, originally designed to handle discontinuities in crack growth models. In particular, by enriching the local spaces with basis functions that are discontinuous across the crack, the method allows for optimal approximation without the need for mesh refinement near the discontinuity. A similar approach was taken for polygonal finite elements [33]. More recently, an enriched Virtual Element method for the Poission problem was designed in [8], however no estimates of the error were given. Following this work, the same authors have proposed in [8] an enriched VEM for a linear elasticity fracture problem. The Hybrid High-Order method is closely linked (c.f. <cit.>) to the nonconforming Virtual Element method (NCVEM) [4]. The recent article [3] designs an enriched NCVEM for harmonic singularities arising from irregular domains in two dimensions. Moreover, the enriched NCVEM is capable of handling highly irregular harmonic singularities, including those arising from cracked domains. While the assumptions we make in Assumption <ref> do not cover cracked domains, the method presented here is robust for all other boundary singularities and has the particular advantage of not requiring the singular functions to be harmonic. This can be particularly useful if the irregularity of the problem is due to singularities in the source term, and not due to non-smooth domains. We also note that the analysis carried out in this paper does not require that inverse inequalities hold for the enriched polynomial spaces.

In this paper we propose an Extended Hybrid High-Order (XHHO) method for the Poisson problem. The work is presented dimension independent, and capable of handling arbitrary singular functions satisfying Assumption <ref>. Specifically, we assume the exact solution consists of a `weakly singular' part lying in a finite dimensional singular space. The local polynomial spaces on the mesh elements and faces are enriched with the appropriate singular space. By then adjusting the local projectors and potential reconstruction accordingly, an optimal XHHO scheme is developed that mimics the standard method <cit.>. However, the analysis is far more involved as we can no longer rely on the Lebesgue/Sobolev embeddings and discrete trace inequalities that apply to polynomial spaces. In Section <ref> we provide a thorough error analysis of the scheme under minimal regularity assumptions. The paper is concluded with an analysis of the choice of stabilisation term in Section <ref> and a discussion on the implementation of the scheme, its numerical limitations and some benchmark tests in Section <ref>.

§.§ Model Problem

Let us consider the typical Dirichlet problem in a polytopal domain \(\Omega\subset\R^d\), \(d \ge 2\):

find \(u\in\HONEzr(\Omega)\) such that

\begin{equation}\label{eq:weak.form}
\a(u,v) = \ell(v),\qquad\forall\,v\in \HONEzr(\Omega) ,
\end{equation}

where \(\a(u,v) \defeq \brac[\Omega]{\nabla u, \nabla v}\), \(\ell(v)\defeq\brac[\Omega]{f,v}\) for some \(f\in \LTWO(\Omega)\). Here and in the following, \(\brac[X]{\cdot, \cdot} \) is the \(\LTWO\)-inner product of scalar- or vector-valued functions on a set \(X\) for its natural measure. To ease the analysis we consider only homogeneous Dirichlet boundary conditions in this paper. However, the scheme extends quite naturally to more general Dirichlet and Neumann problems. Such an extension is outlined in <cit.>.

Consider a partition of the domain \(\Omega\) into a mesh \(\Mh=( \Th,\Fh)\) where the set of mesh elements \(\Th\) are a set of disjoint polytopes such that \(\ol{\Omega} = \bigcup_{T\in\Th}\ol{T}\) and the set of mesh faces \(\Fh\) form the mesh skeleton \(\bigcup_{T\in\Th}{\bdryT} = \bigcup_{F\in\Fh}\ol{F}\). A detailed definition of this structure can be found in <cit.>. The parameter $h$ denotes the maximal element diameter $h\defeq\max_{T\in\Th}\hT$ where, for $X=T\in\Th$ or $X=F\in\Fh$, $h_X$ denotes the diameter of $X$. We shall also collect the set of faces attached to an element $T\in\Th$ in the set $\mathcal{F}_T:=\{F\in\Fh:F\subset T\}$. Similarly, the set containing the one or two elements attached to a face \(F\in\Fh\) is defined as \(\Th[F] \defeq \{T\in\Th:F\subset\bdryT\}\). For each \(T\in\Th\) we denote by \(\norT\) the unit normal directed out of \(T\), and its restriction to a face \(F\in\Fh[T]\) is given by \(\norTF = \norT|_F\). We further make the following assumption on the mesh inline with that stated in [26].

[Connected by star-shaped sets]
There exists a constant \(\varrho>0\) such that for every \(h\in\mathcal{H}\), each \(T\in\Th\) and \(F\in\Fh\) are connected by star-shaped sets with parameter \(\varrho\) (see <cit.>).

It is worth noting that Assumption <ref> is independent of the size of and number of faces in each mesh element. Thus, as in [26], all error estimates in this work remain robust with respect to small and numerous faces.

A typical Hybrid High-Order discretisation of problem (<ref>) relies on piecewise \(H^{k+2}\)-regularity of the solution where \(k\ge0\) is the polynomial degree of the face unknowns. We consider here an exact solution of the form  \( u = u_r + \tilde{u} \) where $u_r$ denotes the `regular part' and $\tilde{u}\in W(\Th)$ where $ W(\Th)$ is a finite dimensional, `weakly singular', enrichment space such that the following assumption holds.

[Assumptions on the enrichment space]
We assume that every \(\psi\in W(\Th)\) satisfies the following conditions:

*  \( \forall T\in \Th\), \(\psi \in H^1(T)\),

*  \( \forall T\in \Th\), \(\Delta \psi \in L^2(T)\),

*  \( \forall T\in \Th\), \(\forall F\in \Fh[T]\), \(\nabla \psi \cdot \norTF \in L^2(F)\).

As mentioned, we assume throughout this paper that the exact solution to (<ref>) can be written as the sum of a regular part and an element of \(W(\Th)\). More specifically, for some \(k\in\N\) we define the solution space to be

\begin{equation}\label{eq:solution.space}
\mathcal{V}^{k+2}(\Omega) \defeq \{w\in H^1_0(\Omega): \Delta w \in L^2(\Omega),\, w \in H^{k+2}(\Th) + W(\Th)\},
\end{equation}

where we denote by \(H^{k+2}(\Th)\) the broken Sobolev space

\begin{equation}
H^{k+2}(\Th) \defeq \{w\in L^{2}(\Omega):w|_T\in H^{k+2}(T)\ \forall T\in\Th\}, \nn
\end{equation}

and assume that the exact solution satisfies $u\in\mathcal{V}^{k+2}(\Omega)$. We note here that Assumption <ref> does not require the enrichment function to be harmonic. Such a case is considered in Section <ref>.

We note that the conditions on both the regular part and the singular part are purely local and the only global conditions enforced on the exact solution are those of \(u\in \HONEzr(\Omega)\) and \(\Delta u \in \LTWO(\Omega)\).

§ DISCRETE PROBLEM

On each $T\in\Th$ we define the discrete space

\begin{equation}\label{def:extended.space}
\POLYX{k+1}(T) \defeq \POLY{k+1}(T) +  W(T),
\end{equation}

where \(W(T)\) denotes the restriction to \(T\) of \(W(\Th)\), and \(\POLY{k+1}(T)\) denotes the space of polynomials on \(T\) of degree no more than \(k+1\), \(k\in\N\). The extended elliptic projector on the space \(\POLYX{k+1}(T)\) is defined as the unique \(\piXTe{k+1}:H^1(T)\to\POLYX{k+1}(T)\) such that for all \(v\in\HONE(T)\)

\begin{equation}\label{eq:elliptic.projector}
\brac[T]{\nabla(v-\piXTe{k+1} v),\nabla w} = 0,\qquad\forall\,w\in\POLYX{k+1}(T)
\end{equation}

\begin{equation}\label{eq:elliptic.projector.closure}
\brac[T]{v-\piXTe{k+1}  v,1} = 0.
\end{equation}

Analogous to <cit.>, we wish to define a discrete space \(\UXTk\), a reconstruction operator \(\pXT{k+1} :\UXTk\to\POLYX{k+1}(T)\), and an interpolator \(\IXTk:H^1(T)\to\UXTk\) such that \(\pXT{k+1} \circ\IXTk = \piXTe{k+1} \). Due to the regularity assumptions on the enrichment space the following integration by parts formula holds for all \(v\in H^1(T)\) and \(w\in\POLYX{k+1}(T)\),

\begin{equation}\label{eq:ibp}
\brac[T]{\nabla v, \nabla w} = -\brac[T]{v, \Delta w} + \brac[\bdryT]{v, \nabla w \cdot \norT}.
\end{equation}

To introduce projectors to (<ref>) we first define the local discrete spaces

\begin{equation}\label{eq:element.and.face.spaces}
\POLYD{k}(T) := \POLY{k}(T) + \Delta  W(T) \quad\textrm{and}\quad \POLYN{k}(F) := \POLY{k}(F) + \sum_{T\in\Th[F]}\nabla  W(T)\cdot\norTF.
\end{equation}

A common application of the method designed in this paper is to irregularities in the solution arising from corners in the domain. In such cases, the singular functions considered are harmonic (see Section <ref>) and the space $\POLYD{k}(T)$ coincides with the polynomial space $\POLY{k}(T)$.

Much of the time, the gradient of an enrichment function \(\psi \in W(\Th)\) is continuous across each face \(F\in\Fh\) and the definition of \(\POLYN{k}(F)\) is equivalent to

\begin{equation}
\POLYN{k}(F) \defeq \POLY{k}(F) + \nabla  W(\Th)\cdot\norF \nn
\end{equation}

for an arbitrary normal \(\norF\) to the face \(F\). However, definition (<ref>) is still well defined for discontinuous enrichment functions. As the space $\POLY{k}(F)$ is defined independent of any particular $T\in\Th$, it is essential in such cases to include the Neumann traces from both elements $T\in\Th[F]$ attached to the face $F$. This is particularly useful when considering locally enriched schemes. Further discussion on this matter is given in Section <ref>.

The discrete broken space on an element boundary is defined as

\[
\POLYN{k}(\Fh[T]) \defeq \{v\in L^1(\bdryT) : v|_F\in \POLYN{k}(F)\quad\forall F\in\Fh[T]\}.
\]

It follows from Assumptions <ref> and <ref> that

\begin{equation}
\POLY{k}(T) \subset \POLYD{k}(T) \subset L^2(T) \quad\textrm{and}\quad \POLY{k}(\Fh[T]) \subset \POLYN{k}(\Fh[T]) \subset L^2(\bdryT). \nn
\end{equation}

We denote by \(\piXTzr{k}\) and \(\piXFTzr{k}\) the \(\LTWO\)-orthogonal projectors on \(\POLYD{k}(T)\) and \(\POLYN{k}(\Fh[T])\) respectively.
Thus, it follows from (<ref>) and (<ref>), as well as the inclusions \(\Delta \POLYX{k+1}(T) \subset \POLYD{k}(T)\) and \(\nabla \POLYX{k+1}(T)\cdot\norT \subset \POLYN{k}(\Fh[T])\), that for all \(v\in H^1(T)\),

\begin{equation}\label{eq:elliptic.projector.formula}
\brac[T]{\nabla \piXTe{k+1} v, \nabla w} = -\brac[T]{\piXTzr{k} v, \Delta w} + \brac[\bdryT]{\piXFTzr{k} v, \nabla w \cdot \norT}\qquad\forall w\in\POLYX{k+1}(T).
\end{equation}

§.§ Local Space

The local space of unknowns is defined to be

\begin{equation}\label{def:local.space}
\UXTk \defeq \POLYD{k}(T) \times \POLYN{k}(\Fh[T]).
\end{equation}

For all \( \ulvXT=(\vXT,\vXFT)\in\UXTk \), the reconstruction operator \( \pXT{k+1} :\UXTk\to\POLYX{k+1}(T) \) is defined by

\begin{equation}\label{eq:reconstruction}
\brac[T]{\nabla \pXT{k+1} \ulvXT, \nabla w} = -\brac[T]{\vXT, \Delta w} + \brac[\bdryT]{\vXFT, \nabla w \cdot \norT} \qquad\forall\,w\in\POLYX{k+1}(T)
\end{equation}

\begin{equation}\label{eq:reconstruction.closure}
\brac[T]{\vXT-\pXT{k+1} \ulvXT,1} = 0.
\end{equation}

Naturally, the interpolator is defined as \( \IXTk v \defeq \brac{\piXTzr{k} v, \piXFTzr{k} v} \). By comparing equations (<ref>) and (<ref>) as well as the closure condition (<ref>) we observe the desired commutation property

\begin{equation}\label{eq:commutation}
\pXT{k+1} \circ\IXTk = \piXTe{k+1} .
\end{equation}

Due to the \(L^2\)-regularity of the unknown spaces, the local reconstruction and interpolator mimic those defined in [18], and the formulation of the discrete problem follows a standard procedure (albeit with non-standard analysis). The continuous form $\a(\cdot,\cdot)$ is approximated on each element by the discrete form $\aXT:\UXTk\times\UXTk\to\R$ defined by

\begin{equation}
\aXT(\uluXT,\ulvXT) \defeq \brac[T]{\nabla \pXT{k+1} \uluXT,\nabla \pXT{k+1} \ulvXT}+\sXT(\uluXT,\ulvXT), \nn
\end{equation}

where $\sXT:\UXTk\times\UXTk\to\R$ is a stabilisation term satisfying the following assumption. From hereon we shall write \(f\lesssim g\) to mean \(f \le Cg\) where \(C\) is a constant depending only on \(\Omega\), \(k\) and the mesh regularity parameter \(\varrho\).

[Local stabilisation term]
The stabilisation term $\sXT$ is a symmetric, positive semi-definite bilinear form that satisfies:

* Coercivity. For all \(\ulvXT\in\UXTk\) it holds that

\begin{equation}\label{eq:stab.vol.coercivity}
\hT^{-2}\norm[T]{\vXT-\pXT{k+1} \ulvXT}^2 \lesssim  \aXT(\ulvXT,\ulvXT)
\end{equation}

\begin{equation}\label{eq:stab.bdry.coercivity}
\hT^{-1}\norm[\bdryT]{\vXFT-\pXT{k+1} \ulvXT}^2 \lesssim  \aXT(\ulvXT,\ulvXT).
\end{equation}

* Consistency. For all \(w=w_r + \psi\) where \(w_r\in H^{k+2}(T)\) and \(\psi\in W(T)\), it holds that

\begin{equation}\label{eq:stab.consistency}
\sXT(\IXTk w, \IXTk w) \lesssim \locerrorRHS{w_r}.
\end{equation}

For the case of regular, polynomial unknowns, the two coercivity conditions (<ref>) and (<ref>) are equivalent to the single coercivity condition stated in <cit.>. However, we have to consider the two conditions here to account for the lack of regularity of the solution.

Some examples of stabilisation terms satisfying Assumption <ref> are given in Section <ref>.

The definition (<ref>) of the extended elliptic projector and the consistency (<ref>) infer for all \(\psi \in W(T)\) the identities

\begin{equation}\label{eq:psi.invariance}
\piXTe{k+1} \psi = \psi \qquad\textrm{and}\qquad \sXT(\IXTk\psi, \IXTk\psi) = 0,
\end{equation}

which, together with (<ref>) and the definition of \(\aXT\), yield

\begin{equation}\label{eq:consistency.on.psi}
\aXT(\IXTk\psi, \IXTk v) = \a(\psi, v)\qquad\forall v \in H^1(T).
\end{equation}

Equation (<ref>) shows the scheme to be exact on elements of \(W(T)\). More generally, \(\psi\) can be replaced by an arbitrary \(w\in\POLYX{k+1}(T)\) in equations (<ref>) and (<ref>) which, as in standard HHO, shows consistency on the discrete space.

§.§ Global Space

The global space of unknowns is defined as

\begin{align}\label{def:global.space}
\UXhkzr\defeq \Big\{\ulvXh=((\vXT)_{T\in\Th},(\vXF)_{F\in\Fh})\,:\,\vXT\in\POLYD{k}(T)\quad\forall T\in\Th, \nl
\vXF\in\POLYN{k}(F)\quad\forall F\in\Fh\,,\vXF=0\quad\forall F \subset \partial \Omega \Big\}.
\end{align}

For any \(\ulvXh\in\UXhkzr\)  we denote its restriction to an element \(T\) by \(\ulvXT = (\vXT , \vXFT ) \in \UXTk\) (where,
naturally, \(\vXFT\) is defined from \((\vXF)_{F\in\Fh[T]}\)). We also denote by \(\vXh\) the piecewise function satisfying \(\vXh|_T=\vXT\) for all \(T\in\Th\). The global problem reads: find \(\uluXh\in\UXhkzr\) such that

\begin{equation}
\label{eq:discrete.problem}
\aXh(\uluXh, \ulvXh) = \brac[\Omega]{f,\vXh}\qquad\forall \ulvXh\in \UXhkzr,
\end{equation}

\begin{equation}
\aXh(\uluXh, \ulvXh) \defeq \sum_{T\in\Th}\aXT(\uluXT, \ulvXT). \nn
\end{equation}

We endow the global space \(\UXhkzr\) with the norm

\begin{equation}
\energynorm{\ulvXh} \defeq \aXh(\ulvXh, \ulvXh)^\frac12. \nn
\end{equation}

We shall also denote by \(\IXhk\) and \(\pXh{k+1}\) the global operators whose restrictions to an element \(T\) are given by  \(\IXTk\) and \(\pXT{k+1}\) respectively. The global stabilisation term is defined as

\begin{equation}
\sXh(\uluXh, \ulvXh) \defeq \sum_{T\in\Th}\sXT(\uluXT, \ulvXT). \nn
\end{equation}

§ ERROR ANALYSIS

Providing a robust error analysis is often the pitfall of enriched schemes as many tools available for polynomial spaces (such as discrete trace and inverse inequalities) no longer apply. The lack of \(H^1\)-regularity of the element unknowns in \(\POLYD{k}\) means we cannot consider their gradients or traces. We are also restricted to considering approximation properties of the extended elliptic projector only in the \(H^1\)-seminorm. Despite these shortcomings, we are able to show consistency of the scheme and provide estimates for the discrete and continuous energy errors.

We first make note of the following continuous trace inequality	(c.f. [26]): for all \(T\in\Th\), \(v\in H^1(T)\),

\begin{equation}
\label{eq:trace.continuous}
\hT\norm[\bdryT]{v}^2 \lesssim \norm[T]{v}^2 + \hT^2\seminorm[H^1(T)]{v}^2.
\end{equation}

The following lemma is also crucial to the consistency of the scheme.

For all \(w\in\mathcal{V}^{k+2}(\Omega)\), \(F\in\Fh\), \(F\not\subset\partial \Omega\),

\begin{equation}\label{eq:gradient.regularity}
\nabla w \cdot \bmn_{T_1 F} + \nabla w \cdot \bmn_{T_2 F} = 0,
\end{equation}

where \(\{T_1,T_2\} = \Th[F]\).

The proof is analogous to that of <cit.> once we note that for all \(w\in\mathcal{V}^{k+2}(\Omega)\) we have: \(\Delta w \in L^2(\Omega)\), \(w \in H^1(\Omega)\), and \(\nabla w \cdot \norT \in L^2(\bdryT)\) for all \(T\in\Th\).

Let \( w = w_r + \psi \in \mathcal{V}^{k+2}(\Omega)\) with \(w_r\in H^{k+2}(\Th)\) and \(\psi\in  W(\Th)\). The consistency error is given by the linear form \(\Eh (w;\cdot):\UXhkzr\to\R\) defined for all \(\ulvXh\in\UXhkzr\) via

\begin{equation}
\Eh (w;\ulvXh) \defeq -\brac[\Omega]{\Delta w, \vXh} - \aXh(\IXhk w,\ulvXh). \nn
\end{equation}

The consistency error satisfies

\begin{equation}\label{eq:error.consistency}
\big| \Eh (w;\ulvXh) \big| \lesssim \energynorm{\ulvXh} h^{k+1}\seminorm[H^{k+2}(\Th)]{w_r}.
\end{equation}

\begin{align}
-\brac[\Omega]{\Delta w, \vXh} \eq \sum_{T\in\Th}-\brac[T]{\Delta w, \vXT} \nl
\eq \sum_{T\in\Th}-\brac[T]{\Delta w, \vXT} + \brac[\bdryT]{\nabla w\cdot\norT, \vXFT}, \nn
\end{align}

where we justify introducing the term

\begin{align}
\sum_{T\in\Th}\brac[\bdryT]{\nabla w\cdot\norT, \vXFT} \eq \sum_{T\in\Th}\sum_{F\in\Fh[T]} \brac[F]{\nabla w\cdot\norTF, \vXF} \nl \eq \sum_{F\in\Fh}\sum_{T\in\Th[F]} \brac[F]{\nabla w\cdot\norTF, \vXF} = 0 \nn
\end{align}

due to equation (<ref>) and the homogeneous condition (<ref>) on the discrete space.
Due to the commutation property (<ref>), the discrete form is given by

\begin{align}
\aXh(\IXhk w,\ulvXh) \eq \sum_{T\in\Th}\Big[\brac[T]{\nabla \piXTe{k+1} w, \nabla \pXT{k+1}  \ulvXT}\Big] + \sXh(\IXhk w,\ulvXh) \nl
\eq \sum_{T\in\Th}\Big[-\brac[T]{\Delta \piXTe{k+1} w, \vXT} + \brac[\bdryT]{\nabla \piXTe{k+1} w\cdot\norT, \vXFT}\Big] + \sXh(\IXhk w,\ulvXh). \nn
\end{align}

\begin{align}\label{proof.error.consistency.exact}
\Eh (w;\ulvXh) \plus \sXh(\IXhk w,\ulvXh) \nl \eq \sum_{T\in\Th}\Big[-\brac[T]{\Delta(w- \piXTe{k+1} w), \vXT} + \brac[\bdryT]{\nabla(w- \piXTe{k+1} w)\cdot\norT, \vXFT}\Big] \nl
\eq \sum_{T\in\Th}\Big[-\brac[T]{\Delta(w_r- \piXTe{k+1} w_r), \vXT} + \brac[\bdryT]{\nabla(w_r- \piXTe{k+1} w_r)\cdot\norT, \vXFT}\Big],
\end{align}

where \(\psi\) has been eliminated via the invariance (<ref>) of \(\piXTe{k+1}\). Let us denote by \(\piTe{k+1}  \) the elliptic projector [18] on the polynomial space \(\POLY{k+1}(T)\subset \POLYX{k+1}(T)\).
On each element \(T\in\Th\) we introduce the elliptic projector as follows,

\begin{align}\label{proof.error.consistency.1}
- \brac[T]{\Delta({}&w_r - \piXTe{k+1} w_r), \vXT} + \brac[\bdryT]{\nabla(w_r- \piXTe{k+1} w_r)\cdot\norT, \vXFT} \nl
\eq-\brac[T]{\Delta(w_r- \piTe{k+1}  w_r), \vXT} + \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT} \nl \minus \brac[T]{\Delta(\piTe{k+1}  w_r- \piXTe{k+1} w_r), \vXT} + \brac[\bdryT]{\nabla(\piTe{k+1}  w_r- \piXTe{k+1} w_r)\cdot\norT, \vXFT}.
\end{align}

As \(\piTe{k+1}  w_r- \piXTe{k+1} w_r\in\POLYX{k+1}(T)\) we may invoke the definition (<ref>) of the potential reconstruction to write

\begin{align}\label{proof.error.consistency.2}
- \brac[T]{\Delta(\piTe{k+1}  w_r- \piXTe{k+1} w_r), \vXT} \plus \brac[\bdryT]{\nabla(\piTe{k+1}  w_r- \piXTe{k+1} w_r)\cdot\norT, \vXFT} \nl
\eq \brac[T]{\nabla(\piTe{k+1}  w_r- \piXTe{k+1} w_r), \nabla \pXT{k+1} \ulvXT}.
\end{align}

Consider also

\begin{align}\label{proof.error.consistency.3}
-\brac[T]{\Delta(w_r \,-\, & \piTe{k+1}  w_r), \vXT} + \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT} \nl
\eq -\brac[T]{\Delta(w_r - \piTe{k+1}  w_r), \vXT - \pXT{k+1}\ulvXT} - \brac[T]{\Delta(w_r - \piTe{k+1}  w_r), \pXT{k+1}\ulvXT} \nl \plus \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT} \nl
\eq -\brac[T]{\Delta(w_r - \piTe{k+1}  w_r), \vXT - \pXT{k+1}\ulvXT} + \brac[T]{\nabla(w_r - \piTe{k+1}  w_r), \nabla \pXT{k+1}\ulvXT} \nl \plus \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT - \pXT{k+1}\ulvXT},
\end{align}

where we have introduced the term \(\pXT{k+1}\ulvXT\) and integrated by parts.
Combining (<ref>), (<ref>) and (<ref>) yields

\begin{align}\label{proof.error.consistency.4}
- \brac[T]{\Delta(w_r \,-\, & \piXTe{k+1} w_r), \vXT} + \brac[\bdryT]{\nabla(w_r- \piXTe{k+1} w_r)\cdot\norT, \vXFT} \nl
\eq -\brac[T]{\Delta(w_r- \piTe{k+1}  w_r), \vXT-\pXT{k+1} \ulvXT} + \brac[T]{\nabla(w_r- \piXTe{k+1} w_r), \nabla \pXT{k+1} \ulvXT} \nl \plus \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT-\pXT{k+1} \ulvXT}.
\end{align}

By the definition (<ref>) of the extended elliptic projector we have that

\begin{equation}\label{proof.error.consistency.5}
\brac[T]{\nabla(w_r- \piXTe{k+1} w_r), \nabla \pXT{k+1} \ulvXT} = 0.
\end{equation}

Therefore, combining (<ref>), (<ref>) and (<ref>) yields

\begin{align}
\Eh (w;\ulvXh) \eq -\sXh(\IXhk w,\ulvXh) + \sum_{T\in\Th}\Big[-\brac[T]{\Delta(w_r- \piTe{k+1}  w_r), \vXT-\pXT{k+1} \ulvXT} \nl \plus \brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT-\pXT{k+1} \ulvXT}\Big], \nn
\end{align}

and thus

\begin{align}\label{proof.error.consistency.6}
\big|\Eh (w;\ulvXh)\big| \lea \sum_{T\in\Th}\big|\sXT(\IXTk w,\ulvXT)\big| + \sum_{T\in\Th}\big|\brac[T]{\Delta(w_r- \piTe{k+1}  w_r), \vXT-\pXT{k+1} \ulvXT}\big| \nl \plus \sum_{T\in\Th}\big|\brac[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)\cdot\norT, \vXFT-\pXT{k+1} \ulvXT}\big|.
\end{align}

By a Cauchy–Schwarz inequality, the coercivity condition (<ref>), and the approximation properties of the elliptic projector <cit.>,

\begin{align}\label{proof.error.consistency.7}
\big|\brac[T]{\Delta(w_r- \piTe{k+1}  w_r), \vXT-\pXT{k+1} \ulvXT}\big| \lea \norm[T]{\Delta(w_r- \piTe{k+1}  w_r)}\norm[T]{\vXT-\pXT{k+1} \ulvXT}\nl
\les \seminorm[H^2(T)]{w_r- \piTe{k+1}  w_r}\hT\aXT(\ulvXT,\ulvXT)^\frac12 \nl
\les \hT^{k+1}\seminorm[H^{k+2}(T)]{w_r}\aXT(\ulvXT,\ulvXT)^\frac12.
\end{align}

Similarly, by a Cauchy–Schwarz inequality, the continuous trace inequality (<ref>), and (<ref>),

\begin{align}\label{proof.error.consistency.8}
\big|\brac[\bdryT]{\nabla(w_r \minus \piTe{k+1}  w_r)\cdot\norT, \vXFT-\pXT{k+1} \ulvXT}\big| \nl \lea \norm[\bdryT]{\nabla(w_r- \piTe{k+1}  w_r)}\norm[\bdryT]{\vXFT-\pXT{k+1} \ulvXT} \nl
\les \Brac{\hT^{-\frac12}\seminorm[H^1(T)]{w_r- \piTe{k+1}  w_r} + \hT^\frac12\seminorm[H^2(T)]{w_r- \piTe{k+1}  w_r}} \hT^\frac12\aXT(\ulvXT,\ulvXT)^\frac12\nl
\les \hT^{k+1}\seminorm[H^{k+2}(T)]{w_r}\aXT(\ulvXT,\ulvXT)^\frac12,
\end{align}

where we have again invoked the approximation properties of the elliptic projector. The stabilisation term is bounded using a Cauchy–Schwarz inequality and the consistency condition (<ref>),

\begin{align}\label{proof.error.consistency.9}
\big|\sXT(\IXTk w,\ulvXT)\big| \lea \sXT(\IXTk w,\IXTk w)^\frac12\sXT(\ulvXT,\ulvXT)^\frac12 \nl \les \hT^{k+1}\seminorm[H^{k+2}(T)]{w_r}\aXT(\ulvXT,\ulvXT)^\frac12.
\end{align}

Substituting (<ref>), (<ref>), and (<ref>) into (<ref>) yields

\begin{equation}
\big|\Eh (w;\ulvXh)\big| \lesssim \sum_{T\in\Th}\hT^{k+1}\seminorm[H^{k+2}(T)]{w_r}\aXT(\ulvXT,\ulvXT)^\frac12. \nn
\end{equation}

The proof follows from a discrete Cauchy–Schwarz inequality and noting that \(\hT \le h\) for all \(T\in\Th\).

The proof of the consistency error given above is significantly more detailed than that of regular HHO given in <cit.>. The reason for this is partly due the limited regularity of element unknowns, \(\vXT\in L^2(T)\). However, even if we were to assume \(\HONE\)-regularity to arrive at the equation

\begin{multline*}
\Eh (w;\ulvXh) + \sXh(\IXhk w,\ulvXh) = \sum_{T\in\Th}\Big[\brac[\T]{\nabla(w_r - \piXTe{k+1} w_r), \nabla\vXT} \\
+ \brac[\bdryT]{\nabla(w_r - \piXTe{k+1} w_r)\cdot\norT, \vXFT - \vXT}\Big],
\end{multline*}

this is still not useful. The first term does not equate to zero as $\vXT \notin \POLYX{k+1}(T)$ in general. Moreover, there is no guarantee that \(\norm[\bdryT]{\nabla(w_r - \piXTe{k+1} w_r)\cdot\norT}\) will scale appropriately. Thus, the extra and lengthy details in the proof of Theorem <ref> appear necessary.

Let \(u=u_r+\tilde{u} \in \mathcal{V}^{k+2}(\Omega)\) be the exact solution to the continuous problem (<ref>) where \(u_r\in H^{k+2}(\Th)\) and \(\tilde{u}\in W(\Th)\). Let $\uluXh\in\UXhkzr$ be the solution to the discrete problem (<ref>). The following energy error estimates hold:

\begin{equation} \label{eq:energy.error}
\energynorm{\uluXh - \IXhk u} + \seminorm[H^1(\Th)]{\pXh{k+1} \uluXh - u} \lesssim h^{k+1}\seminorm[H^{k+2}(\Th)]{u_r}.
\end{equation}

By the coercivity conditions (<ref>) and (<ref>), as well as the homogeneous conditions on the discrete space, it is clear that $\energynorm{\cdot}$ describes a norm on $\UXhkzr$. As such, we infer from the Third Strang Lemma [17] that

\begin{equation}
\energynorm{\uluXh - \IXhk u} \le \sup_{\ulvXh\ne0}\frac{|\mathcal{E}_h( u;\ulvXh)|}{\energynorm{\ulvXh}}. \nn
\end{equation}

Combining with the consistency error (<ref>) yields the estimate

\begin{equation}\label{eq:energy.error.discrete}
\energynorm{\uluXh - \IXhk u} \lesssim h^{k+1}\seminorm[H^{k+2}(\Th)]{u_r}.
\end{equation}

Consider on each element $T\in\Th$ the triangle inequality,

\begin{align}
\seminorm[H^1(T)]{\pXT{k+1}  \uluXT - u} \lea
\seminorm[H^1(T)]{\pXT{k+1}  \uluXT - \piXTe{k+1}  u} + \seminorm[H^1(T)]{\piXTe{k+1}  u - u} \nl
\eq \seminorm[H^1(T)]{\pXT{k+1}  (\uluXT - \IXTk u)} + \seminorm[H^1(T)]{\piXTe{k+1}  u_r - u_r}\nl
\lea \aXT(\uluXT - \IXTk u, \uluXT - \IXTk u)^\frac12 + \seminorm[H^1(T)]{\piTe{k+1}  u_r - u_r},\nn
\end{align}

where \(\piXTe{k+1}\) is replaced by \(\piTe{k+1}\) due to orthogonal projectors minimising their respective norms and \(\piTe{k+1}  u_r\in\POLYX{k+1}(T)\).
Squaring, summing over all \(T\in\Th\) and invoking the approximation properties of the elliptic projector yields

\begin{equation}
\seminorm[H^1(\Th)]{\pXh{k+1} \uluXh - u}^2 \lesssim \energynorm{\uluXh - \IXhk u}^2 + \big[h^{k+1}\seminorm[H^{k+2}(\Th)]{u_r}\big]^2. \nn
\end{equation}

The proof is complete by applying the estimate (<ref>).

§ STABILISATION

In this section, we give an example of a stabilisation term satisfying the coercivity and consistency properties of Assumption <ref>. We first give below an extension of <cit.> to the extended discrete spaces considered here. The difference operators \(\deltaXT{k}:\UXTk\to\POLYD{k}(T)\) and \(\deltaXFT{k}:\UXTk\to\POLYN{k}(\Fh[T])\) are defined via

\begin{equation*}
\deltaXT{k} \ulvXT \defeq \vXT-\piXTzr{k} \pXT{k+1} \ulvXT\qquad\textrm{and}\qquad \deltaXFT{k} \ulvXT \defeq \vXFT-\piXFTzr{k} \pXT{k+1} \ulvXT.
\end{equation*}

For \(\sXT\) to satisfy the consistency condition (<ref>) it is necessary (but not sufficient) that \(\sXT\) depends on its arguments only through $\deltaXT{k}$ and $\deltaXFT{k}$.

Consider a symmetric, positive semi-definite bilinear form \(\sXT\) that satisfies condition (<ref>). Thus, for all \(w=w_r+\psi\in \POLYX{k+1}(T)\) with \(w_r\in\POLY{k+1}(T)\) and \(\psi\in W(T)\),

\begin{equation}
\sXT(\IXTk w, \IXTk w) \lesssim \locerrorRHS{w_r} =  0. \nn
\end{equation}

Therefore, by a Cauchy–Schwarz inequality,

\begin{equation}
\sXT(\IXTk w, \ulvXT) = 0 \nn
\end{equation}

for all \(\ulvXT\in\UXTk\). Thus,

\begin{align}
\sXT(\uluXT, \ulvXT) \eq \sXT(\uluXT, \ulvXT - \IXTk\pXT{k+1} \ulvXT) \nl
\eq \sXT(\uluXT - \IXTk\pXT{k+1} \uluXT, \ulvXT - \IXTk\pXT{k+1} \ulvXT). \nn
\end{align}

The proof is complete by noting that

\begin{equation}
\ulvXT-\IXTk\pXT{k+1} \ulvXT = \brac{\deltaXT{k} \ulvXT, \deltaXFT{k} \ulvXT}. \nn
\end{equation}

A stabilisation term satisfying Assumption <ref> is obtained setting

\begin{equation}
\label{def:vol.stab}
\sXT(\uluXT, \ulvXT) \defeq \hT^{-2}\brac[T]{\deltaXT{k} \uluXT, \deltaXT{k} \ulvXT} + \hT^{-1}\brac[\bdryT]{\deltaXFT{k} \uluXT, \deltaXFT{k} \ulvXT}.
\end{equation}

The stabilisation term defined by (<ref>) satisfies the coercivity conditions (<ref>) and (<ref>).

Consider by a triangle inequality, and the inclusion \(\POLY{k}(T)\subset\POLY{k}_\Delta(T)\) along with the minimisation of \(\piXTzr{k}\) on \(\POLYD{k}(T)\),

\begin{align}
\hT^{-2}\norm[T]{\vXT-\pXT{k+1} \ulvXT}^2 \les \hT^{-2}\norm[T]{\vXT-\piXTzr{k} \pXT{k+1} \ulvXT}^2 + \hT^{-2}\norm[T]{\piXTzr{k} \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 \nl
\lea \sXT(\ulvXT, \ulvXT) + \hT^{-2}\norm[T]{\piTzr{k}  \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 \nl
\les \sXT(\ulvXT, \ulvXT) + \seminorm[H^1(T)]{\pXT{k+1} \ulvXT}^2 = \aXT(\ulvXT, \ulvXT), \nn
\end{align}

where the final inequality follows from the approximation properties of the \(L^2\)-orthogonal projector on \(\POLY{k}(T)\) <cit.>.
Similarly, by a triangle inequality and noting that \(\POLY{k}(T)|_{\bdryT} \subset \POLYN{k}(\Fh[T])\) (so to replace \(\piXFTzr{k}\) with \(\piTzr{k}\)),

\begin{align}
\hT^{-1}\norm[\bdryT]{&\vXFT-\pXT{k+1} \ulvXT}^2 \nl \les \hT^{-1}\norm[\bdryT]{\vXFT-\piXFTzr{k} \pXT{k+1} \ulvXT}^2 + \hT^{-1}\norm[\bdryT]{\piXFTzr{k} \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 \nl
\lea \sXT(\ulvXT, \ulvXT) + \hT^{-1}\norm[\bdryT]{\piTzr{k}  \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 \nl
\les \sXT(\ulvXT, \ulvXT) + \hT^{-2}\norm[T]{\piTzr{k}  \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 + \seminorm[H^1(T)]{\piTzr{k}  \pXT{k+1} \ulvXT-\pXT{k+1} \ulvXT}^2 \nl
\les \sXT(\ulvXT, \ulvXT) + \seminorm[H^1(T)]{\pXT{k+1} \ulvXT}^2 = \aXT(\ulvXT, \ulvXT), \nn
\end{align}

where we have used the continuous trace inequality (<ref>) and again invoked the approximation properties of the \(L^2\)-orthogonal projector.

The stabilisation term defined by (<ref>) satisfies the consistency condition (<ref>).

By the definition of \(\sXT\), the commutation property (<ref>), and the boundedness of orthogonal projectors, it holds that

\begin{align}
\sXT(\IXTk w, \IXTk w) \eq \hT^{-2}\norm[T]{\piXTzr{k} w-\piXTzr{k} \piXTe{k+1} w}^2 + \hT^{-1}\norm[\bdryT]{\piXFTzr{k} w-\piXFTzr{k} \piXTe{k+1} w}^2 \nl
\lea \hT^{-2}\norm[T]{w-\piXTe{k+1} w}^2 + \hT^{-1}\norm[\bdryT]{w-\piXTe{k+1} w}^2. \nn
\end{align}

By the continuous trace inequality (<ref>) and a Poincaré inequality due to the zero mean value of \(w-\piXTe{k+1} w\), we infer that

\begin{align}
\sXT(\IXTk w, \IXTk w) \les \hT^{-2}\norm[T]{w-\piXTe{k+1} w}^2 + \seminorm[H^1(T)]{w-\piXTe{k+1} w}^2 \nl
\les \seminorm[H^1(T)]{w-\piXTe{k+1} w}^2. \nn
\end{align}

We note the inclusion \(\POLY{k+1}(T)\subset \POLYX{k+1}(T)\) and invoke the invariance and minimisation properties of orthogonal projectors to conclude that

\begin{equation}
\seminorm[H^1(T)]{w-\piXTe{k+1} w}^2 = \seminorm[H^1(T)]{w_r-\piXTe{k+1} w_r}^2 \le \seminorm[H^1(T)]{w_r-\piTe{k+1}  w_r}^2. \nn
\end{equation}

The proof then follows from the approximation properties of the elliptic projector <cit.>.

Designing alternate stabilisation terms proves difficult due to the limited regularity of the unknowns. However, if we assume that \(\Delta W(T) \subset \POLY{k}(T)\) (so that the element unknowns are polynomials), the stabilisation terms

\begin{align}
\sXT^\nabla(\uluXT, \ulvXT) \deq \brac[T]{\nabla \deltaXT{k}\uluXT, \nabla \deltaXT{k}\ulvXT} + \hT^{-1}\brac[\bdryT]{\deltaXFT{k}\uluXT, \deltaXFT{k}\ulvXT} , \label{eq:stab.grad}\\
\sXT^\partial(\uluXT, \ulvXT) \deq \hT^{-1}\brac[\bdryT]{(\deltaXFT{k}-\deltaXT{k})\uluXT, (\deltaXFT{k}-\deltaXT{k})\ulvXT} \label{eq:stab.bdry}
\end{align}

as defined in <cit.> both satisfy Assumption <ref>. While this may at first seem contrived, it is quite natural to consider a singular enrichment function with zero Laplacian (see Section <ref>) which clearly satisfies the aforementioned condition.

§ IMPLEMENTATION AND NUMERICAL TESTS

As mentioned in the introduction, singular solutions to (<ref>) can arise from corners in the boundary of an otherwise smooth domain. In particular, at non-convex corners we cannot even assume \(H^2\)-regularity of the solution <cit.>. Moreover, singularities on the boundary can arise from irregular boundary data, or a transition from Dirichlet to Neumann data. For simplicity, we consider here a domain \(\Omega\subset \R^2\) with one re-entrant corner located at the origin. Thus, after a possible rotation of coordinates, the domain \(\Omega\) corresponds with the region \(\{(x_1, x_2) = (r \cos \theta, r \sin \theta)|\,r > 0, 0 < \theta < \omega\}\) in some neighbourhood of the origin, where \(\pi < \omega \le 2\pi\) is the angle of the re-entrant corner. For each \(j\in\N\) we define a function

\begin{equation}\label{eq:singular.functions}
u_j =
\begin{cases}
r^\frac{j\pi}{\omega}\sin(\frac{j\pi}{\omega} \theta)&\qquad\textrm{if } \frac{j \pi}{\omega} \not\in \Z \\
r^\frac{j\pi}{\omega}(\ln r\sin(\frac{j\pi}{\omega}\theta)+\theta \cos(\frac{j\pi}{\omega}\theta))&\qquad\textrm{if } \frac{j \pi}{\omega} \in \Z
\end{cases}.
\end{equation}

Given a source term \(f\in H^k(\Omega)\), there exist numbers \(c_j\) such that the solution to the homogeneous Dirichlet problem (<ref>) satisfies

\begin{equation}\label{eq:regularity.result}
u - \sum_{1 \le j < \frac{\omega}{\pi}(k+1)} c_j u_j \in H^{k+2}(\Omega_0),
\end{equation}

where \(\Omega_0\subset \Omega\) is some open neighbourhood of the origin. We refer the reader to <cit.>, which is dedicated to proving (<ref>) and equivalent results on polygonal domains with generic boundary data. Each singular function \(u_j\) clearly satisfy Assumptions <ref> and <ref> (the latter due to \(\Delta u_j = 0\)). On an edge \(F\in\Fh\) containing the singular point \(r=0\), the least regular function (\(j=1\)) satisfies \(\nabla u_1 \cdot \norF \in L^p(F)\) for all \(p < \frac{\omega}{\omega-\pi}\). Therefore, Assumption <ref> holds true provided that \(\omega < 2\pi\) (corresponding to a cracked domain).

In practice, to assure \(H^{k+2}\)-regularity in a polygonal domain, singular functions require to be defined at every corner. However, for computational simplicity, we consider only one singular function defined at the re-entrant corner. In particular, we consider here an L-shape domain \( \Omega = (-1,1)^2\,\backslash\,[0,1]^2 \) and exact solution

\begin{equation}
u = \sin(\pi x_1) \sin(\pi x_2) + \psi, \nn
\end{equation}

where \(\psi = r^\alpha \sin(\alpha (\theta - \frac{\pi}{2}))\), \(\alpha = \frac{\pi}{\omega} = \frac23\) and \(\frac{\pi}{2} \le \theta \le 2 \pi\). Naturally, we define the enrichment space as \(W(\Th) = \textrm{span}\{\psi\}\). We remark that the exact solution considered does not have homogeneous boundary conditions. However, as mentioned in the introduction, the extension of the XHHO scheme to inhomogeneous boundary data follows seamlessly. Indeed, if we consider the inhomogeneous condition $u=g_D$ on $\partial \Omega$, we define $\ul{u}_{\xh,D}=\brac{\brac[T\in\Th]{0}, \brac[F\in\Th]{u_{\xF,D}}}\in\UXhk$ where

\[
u_{\xF,D} =
\begin{cases}
\piXFzr{k}g_D&\ \textrm{if } F\subset\partial \Omega \\
0&\ \textrm{otherwise}
\end{cases}. \nn
\]

The inhomogeneous problem is given by: find $\ul{u}_{\xh,0}\in\UXhkzr$ such that

\[
\aXh(\ul{u}_{\xh,0}, \ulvXh) = \brac[\Omega]{f, \vXh} - \aXh(\ul{u}_{\xh,D}, \ulvXh) \qquad \forall\ulvXh\in\UXhkzr,
\]

and the discrete solution is given by $\uluXh = \ul{u}_{\xh,0} + \ul{u}_{\xh,D}$.

§.§ Basis Functions

The Extended Hybrid High-Order scheme designed in this paper requires the computation of spaces \(\POLYX{k+1}(T)\) and \(\POLYD{k}(T)\) on each element \(T\in\Th\), and of the space \(\POLYN{k}(F)\) on each face \(F\in\Fh\). As the Laplacian of the singular function is \(0\), i.e. \(\Delta W(\Th) = \{0\}\), it holds that \(\POLYD{k}(T) = \POLY{k}(T) \subset \POLYX{k+1}(T)\). Therefore, on each mesh element we define basis functions for the space \(\POLYX{k+1}(T)\) and consider \(\POLYD{k}(T)\) a subspace.

Consider monomial basis functions, \(\{\phi_T^j\}_{j=1}^{\textrm{dim}\{\POLY{k+1}(T)\}}\) and \(\{\phi_F^j\}_{j=1}^{\textrm{dim}\{\POLY{k}(F)\}}\), in locally scaled coordinates for the respective polynomial spaces. On each mesh element we write the extended space as

\begin{equation}\label{eq:cell.basis}
\POLYX{k+1}(T) = \textrm{span}\big\{ (\phi_T^j)_{j=1}^{\textrm{dim}\{\POLY{k+1}(T)\}}, \psi \big\},
\end{equation}

and we then \(L^2\)-orthonormalise the basis functions following a Gram-Schmidt process. It is natural to define a basis for the enriched space on the mesh faces as

\begin{equation}\label{eq:face.basis}
\POLYN{k}(F) = \textrm{span}\big\{ (\phi_F^j)_{j=1}^{\textrm{dim}\{\POLY{k}(F)\}}, \nabla \psi \cdot \norF \big\},
\end{equation}

for some choice of normal vector \(\norF\). However, there may exist certain faces \(F\in\Fh\) such that this choice of basis is not linearly independent. Consider a unit normal to a face \(F\in\Fh\) defined by \(\norF = (\cos \theta_\nor, \sin \theta_\nor)\) for some constant angle \(\theta_\nor\). Taking the scalar product with \(\nabla \psi\),

\begin{equation}\label{eq:grad.dot.n}
\nabla \psi \cdot \norF = \alpha r^{\alpha - 1} \sin(\theta_\nor - \theta + \alpha(\theta - \frac{\pi}{2})).
\end{equation}

If \(\theta \equiv const\) then \(\theta_{\nor} = \theta + (2 j + 1) \frac{\pi}{2} \), for an integer \(j \in \N \). Substituting into (<ref>) yields

\begin{equation}
\nabla \psi \cdot \norF =  \pm\alpha r^{\alpha - 1} \cos(\alpha(\theta - \frac{\pi}{2})) \nn
\end{equation}

which is identically \(0\) if \(\theta \equiv \frac{(6j-1)\pi}{4}\) for some \(j\in\N\). Thus, as we only consider \(\frac{\pi}{2} \le \theta \le 2\pi\), the basis for \(\POLYN{k}(F)\) defined by (<ref>) is linearly independent except along faces such that \(\theta = const = \frac{5\pi}{4}\) (for which \(\POLYN{k}(F) = \POLY{k}(F)\)). The bases on each face are then also \(L^2\)-orthonormalised via a Gram-Schmidt process.

§.§ Integration Rules

It is well known that classical Gauss-Legendre quadrature rules suffer from large errors and numerical instability near a singularity. To design an alternate numerical integration rule we make use of the fact that \(\psi\) is a homogeneous function of degree \(\alpha\). Thus, we can use a homogeneous integration rule (as described in [14]) to rewrite the volumetric integrals on the boundary. In particular we consider the integral \( \int_T \psi \phi^j_T \)
for some \(1 \le j \le \textrm{dim}\{\POLY{k+1}(T)\}\). As \(\phi_T^j\) is a monomial in local coordinates \(\tilde{\bmx} = \frac{\bmx - \bmx_0}{\hT}\), it is homogeneous in \(\bmx - \bmx_0\). By Euler's homogeneous function theorem, \(\nabla \psi \cdot \bmx = \alpha \psi\) and \(\nabla \phi_T^j \cdot (\bmx-\bmx_0) = q \phi_T^j\), where \(q\) is the degree of the monomial \(\phi_T^j\). Therefore,

\begin{align}
\int_{\bdryT} \psi \phi^j_T \bmx\cdot\norT \eq \int_T \psi \phi^j_T \nabla \cdot \bmx + \int_T \bmx \cdot\nabla (\psi \phi^j_T) \nl
\eq (d + \alpha + q)\int_T \psi \phi^j_T + \int_T \psi\bmx_0 \cdot\nabla \phi^j_T, \nn
%	 \label{eq:homogeneous.integration}
\end{align}

which yields

\begin{equation}
(d + \alpha + q)\int_T \psi \phi^j_T = \sum_{F\in\Fh[T]}\bmx_F\cdot\norTF\int_F\psi \phi^j_T - \int_T \psi\bmx_0 \cdot\nabla \phi^j_T, \nn
\end{equation}

where \(\bmx_F\) is an arbitrary point in the face \(F\). As each component of \(\nabla \phi^j_T\) is a monomial of degree \(q- 1\), the above process can be repeated iteratively until the whole integral is described on the boundary. For edges that do not pass through the singular point \(r=0\), the integral \(\int_F\psi \phi^j_T\) can be accurately approximated using Gauss-Legendre quadrature. On edges \(F\in\Fh\) with \(\theta \equiv const\) (which include edges passing through the origin) the integral can be computed exactly. Consider the arc-length parametrisation \(\bmx = (r \cos(\theta), r \sin(\theta))\), \(R_0 \le r \le R_1\) and write

\begin{equation}
\int_F\psi \phi^j_T = \int_{R_1}^{R_2} \psi(\bmx(r)) \phi_T^j(\bmx(r))\,dr. \nn
\end{equation}

An integration by parts yields

\begin{align}
r\psi(\bmx(r)) \phi_T^j(\bmx(r))\Big|_{r=R_1}^{R_2} \eq \int_{R_1}^{R_2} \psi  \phi_T^j \frac{dr}{dr}\,dr + \int_{R_1}^{R_2} r\frac{d}{dr}(\psi \phi_T^j)\,dr \nl
\eq \int_{R_1}^{R_2} \psi \phi_T^j\,dr + \int_{R_1}^{R_2} r(\phi_T^j \nabla \psi \cdot \frac{d\bmx}{dr} + \psi \nabla \phi_T^j \cdot \frac{d\bmx}{dr})\,dr. \label{eq:integration.on.faces}
\end{align}

Noting that \(\bmx = r \frac{d\bmx}{dr}\), and again invoking the homogeneity results \(\nabla \psi \cdot \bmx = \alpha \psi\) and \(\nabla \phi_T^j \cdot (\bmx - \bmx_0) = q \phi_T^j\), equation (<ref>) is evaluated as

\begin{equation}
\norm{\bmx}\psi(\bmx) \phi_T^j(\bmx)\Big|_{\bmx=\bmv_0}^{\bmv_1} = (1 + \alpha + q)\int_{R_1}^{R_2} \psi \phi_T^j\,dr + \int_{R_1}^{R_2} \psi \bmx_0 \cdot \nabla \phi_T^j\,dr, \nn
\end{equation}

where \(\bmv_0\) and \(\bmv_1\) are the vertices corresponding to \(r=R_0\) and \(r=R_1\) respectively. Again, we may follow an iterative procedure for each component of \(\int_{R_1}^{R_2} \psi \bmx_0 \cdot \nabla \phi_T^j\,dr\) to evaluate the entire integral \(\int_F\psi \phi^j_T\) on the vertices \(\bmv_0\), \(\bmv_1\).

The calculations of integrals involving gradients is marginally simpler as they can be transformed directly onto the boundary via

\begin{equation}
\int_T \nabla \psi \cdot \nabla \phi_T^j = \int_{\bdryT} \phi_T^j \nabla \psi \cdot \norT, \nn
\end{equation}

(due to \(\Delta \psi = 0\)). Each edge integral is calculated using the process described above, noting that each \(\nabla \psi \cdot \norTF\) is a homogeneous function of degree \(\alpha - 1\).

§.§ Local Enrichment

Far from the singular point \(r=0\), the enrichment function \(\psi\) is smooth. More precisely, for any open set \(V\) containing the point \(r=0\), \(\psi \in C^\infty(\Omega\backslash V)\). As such, away from the singularity, \(\psi\) can be efficiently approximated by a polynomial. This can cause the local matrices appearing in the orthonormalisation process of the bases (<ref>) and (<ref>) to be highly ill conditioned. Naturally, this leads us to consider a scheme such that the unknowns spaces are enriched near the singular point and are otherwise polynomial spaces. We follow a similar approach to that taken in <cit.>.

Consider a parameter \(\gamma > 0 \) and define the set

\begin{equation}
\label{eq:enriched.elements}
\mathcal{T}_{\gamma} \defeq \{ T \in \Th: \norm{{\bm{0}} - \bmx_T} < \gamma \},
\end{equation}

where we denote by \(\bmx_T\) the centroid of an element \(T\). A cut-off function \(\delta_{\gamma}:\Omega\to\R\) is defined by

\begin{equation}
\delta_{\gamma}|_T =
\begin{cases}
1&\ \textrm{if } T \in \mathcal{T}_{\gamma} \\
0&\ \textrm{otherwise}
\end{cases}. \nn
\end{equation}

If \(u=u_r + \psi\), with \(u_r \in C^\infty(\Omega)\) we can write \(u = u_r + \psi - \delta_{\gamma}\psi + \delta_{\gamma}\psi = \tilde{u}_r + \delta_{\gamma}\psi\) where \(u_r + \psi - \delta_{\gamma}\psi = \tilde{u}_r \in H^{k+2}(\Th)\). Thus, we define the enrichment space as \(W(\Th) = \textrm{span}\{\delta_{\gamma}\psi\}\). As the assumptions on the enrichment space are all made locally, \(W(\Th)\) clearly still satisfies Assumption <ref>. Definitions (<ref>) and (<ref>) of the discrete spaces then correspond to

\begin{align}
\POLYX{k+1}(T) \eq
\begin{cases}
\POLY{k+1}(T) + \textrm{span}\{\psi\}&\ \textrm{if } T \in \mathcal{T}_{\gamma} \\
\POLY{k+1}(T)&\ \textrm{otherwise}
\end{cases}\nl
\POLYN{k}(F) \eq
\begin{cases}
\POLY{k}(F) + \textrm{span}\{\nabla\psi\cdot\norF\}&\ \textrm{if } F \in \mathcal{F}_{\gamma} \\
\POLY{k}(F)&\ \textrm{otherwise}
\end{cases} \nn
\end{align}

for an arbitrary normal \(\norF\), where

\begin{equation}
\mathcal{F}_{\gamma} = \{ F \in \Fh: \Th[F]\cap\mathcal{T}_\gamma \ne \emptyset \}. \nn
\end{equation}

§.§ Tests

As mentioned previously, we consider an  L-shape domain \( \Omega = (-1,1)^2\,\backslash\,[0,1]^2 \) and exact solution \(u = \sin(\pi x_1) \sin(\pi x_2) + r^\alpha \sin(\alpha( \theta - \frac{\pi}{2}))\), $\alpha=\frac23$. The stabilisation term defined by (<ref>) is considered. As the element unknowns are polynomials, more choices of stabilisation terms exist (see Section <ref>). A comparison between the various stabilisation terms is carried out for the standard HHO method in <cit.>. The XHHO scheme is tested on a sequence of hexagonal meshes with maximum element diameter \(h\to 0\). Two members of this family are plotted in Figure <ref> and the mesh data is shown in Table <ref>.

Two members of the mesh sequence

columns/MeshSize/.style=column name=$h$,/pgf/number format/.cd,fixed zerofill,precision=4,
columns/NbCells/.style=column name=Nb. Elements,
columns/NbInternalEdges/.style=column name=Nb. Internal Edges,
every head row/.style=before row=,after row=,
every last row/.style=after row=

]Parameters of the mesh sequence

The conditioning of the system poses a significant challenge to the scheme. This is a common problem for enriched schemes and the matter is discussed in detail for the enriched NCVEM [3]. For a \(T\in\Th\), let us denote by \( \{\phi_T^j\}_{j=1}^{\dim\{\POLYX{k+1}(T)\}} \) the basis for \(\POLYX{k+1}(T)\) prior to orthonormalising. The mass matrix is defined by

\begin{equation}
(\mat{M}_T)_{i,j} = \brac[T]{\phi_T^i,\phi_T^j}. \nn
\end{equation}

The maximum condition number of the element mass matrices is defined by \(\mathcal{C}\defeq\max_{T\in\Th} \frac{\lambda_{T,max}}{\lambda_{T,min}}\), where \(\lambda_{T,max}\) and \(\lambda_{T,min}\) denote the maximum and minimum eigen values of \(\mat{M}_T\) respectively. The parameter \(\mathcal{C}\) gives a measure of how linearly independent the worst performing basis is. In Table <ref> we present the maximum condition number on the finest mesh in the sequence with various values of \(k\) and \(\gamma\) (where \(\gamma=0\) corresponds to an non-enriched scheme). It is clear that the conditioning of the scheme is significantly worse for enriched schemes than non-enriched schemes and gets progressively worse with increasing \(k\) and \(\gamma\). We note that once \(\mathcal{C}\sim 10^{16}\) the orthonormalisation process (and thus the scheme itself) fails due to division by numbers which are numerically zero.

columns/CellCondition/.style=column type/.add=|,column name=\(k\eqs0\),/pgf/number format/sci,zerofill,precision=3,
columns/Gamma/.style=string type, column type = l,column name=,
every head row/.style=after row=
columns/CellCondition/.style=column name=\(k\eqs1\),/pgf/number format/sci,zerofill,precision=3,
every head row/.style=after row=
columns/CellCondition/.style=column name=\(k\eqs2\),/pgf/number format/sci,zerofill,precision=3,
every head row/.style=after row=
columns/CellCondition/.style=column name=\(k\eqs3\),/pgf/number format/sci,zerofill,precision=3,
every head row/.style=after row=
]Maximum local condition number with \(h\approx 0.0257\)

Denote by \(\uluXh \in\UXhkzr\) the exact solution to the discrete problem (<ref>). The relative error of the scheme is determined via the following three quantities,

\begin{align}
E_{0,\Th} \deq \Big[\frac{\sum_{T\in\Th}\norm[T]{\uXT-\piXTzr{k} u}^2}{\sum_{T\in\Th}\norm[T]{\piXTzr{k} u}^2}\Big]^\frac12 + \Big[\frac{\sum_{F\in\Fh}h_F\norm[F]{\uXF-\piXFzr{k} u}^2}{\sum_{F\in\Fh}h_F\norm[F]{\piXFzr{k} u}^2}\Big]^\frac12 \nl
E_{1,\Th} \deq \Big[\frac{\sum_{T\in\Th}\seminorm[H^1(T)]{\pXT{k+1}\uluXT - \piXTe{k+1}  u}^2}{\sum_{T\in\Th}\seminorm[H^1(T)]{\piXTe{k+1}  u}^2}\Big]^\frac12 \nl
E_{\a,\xh} \deq \frac{\energynorm{\uluXh - \IXhk u}}{\energynorm{\IXhk u}}. \nn
\end{align}

Unlike [3], we do not observe saturation in error rates as the condition number gets large, but rather an instantaneous failure of the scheme. As such we would like to choose a cut-off value \(\gamma\) that is as large as possible and does not result in system failure. In Figures <ref> and <ref> we test convergence of the scheme with \(k=1\) and \(k=2\) respectively. The enriched scheme performs significantly better than the non-enriched scheme, particularly in \(H^1\)-error (\(E_{1,\Th}\)) and for higher polynomial degree \(k\). In Figure <ref>, we test two cut-off values (\(\gamma=0.075\) and \(\gamma=0.15\)). The larger cut-off value performs slightly better in \(H^1\)-error, however, the scheme fails on the final mesh.

$E_{0,\Th}$ vs $h$

[ legend columns=3, legend to name=k1 ]
table[x=MeshSize,y=L2Error] data/unenriched_k1.dat;
table[x=MeshSize,y=L2Error] data/local_k1_cut_015.dat;
Non-enriched, Locally enriched $(\gamma = 0.15)$;

$E_{1,\Th}$ vs $h$

table[x=MeshSize,y=H1Error] data/unenriched_k1.dat;
table[x=MeshSize,y=H1Error] data/local_k1_cut_015.dat;

$E_{\a,\xh}$ vs $h$

table[x=MeshSize,y=EnergyError] data/unenriched_k1.dat;
table[x=MeshSize,y=EnergyError] data/local_k1_cut_015.dat;

Error vs $h$, \(k=1\)

$E_{0,\Th}$ vs $h$

[ legend columns=3, legend to name=k2 ]
table[x=MeshSize,y=L2Error] data/unenriched_k2.dat;
table[x=MeshSize,y=L2Error] data/local_k2_cut_0075.dat;
table[x=MeshSize,y=L2Error] data/local_k2_cut_015.dat;
Non-enriched, Locally enriched ($\gamma = 0.075$), Locally Enriched ($\gamma = 0.15$);

$E_{1,\Th}$ vs $h$

table[x=MeshSize,y=H1Error] data/unenriched_k2.dat;
table[x=MeshSize,y=H1Error] data/local_k2_cut_0075.dat;
table[x=MeshSize,y=H1Error] data/local_k2_cut_015.dat;

$E_{\a,\xh}$ vs $h$

table[x=MeshSize,y=EnergyError] data/unenriched_k2.dat;
table[x=MeshSize,y=EnergyError] data/local_k2_cut_0075.dat;
table[x=MeshSize,y=EnergyError] data/local_k2_cut_015.dat;

Error vs $h$, \(k=2\)

A standard error estimate with respect to $k$ for a classical HHO scheme is given by

\begin{equation}\label{eq:error.estimate.k}
\norm[\a,h]{\uluh - \Ihk u} \lesssim \frac{h^{k+1}}{(k+1)^k} \seminorm[\HS{k+2}(\Th)]{u}
\end{equation}

where the hidden constant is independent of $k$ <cit.>. Therefore, if

\begin{equation}\label{eq:sovolev.seminorm.bound}
\frac{\seminorm[\HS{k+2}(\Th)]{u}}{(k+1)^k} \le C
\end{equation}

for some $C$ independent of $k$, then the scheme should converge exponentially with respect to $k$.
In Figure <ref>, we fix the mesh (the second mesh in Table <ref>) and test convergence as \(k\) increases. In this case, the polynomial spaces are enriched only on the element containing the singular point and on its faces. While convergence with respect to $k$ has not been proven in this work, equation (<ref>) provides a bench mark to compare the tests to. As such, we plot on a log-linear scale and include a line of slope $\log h$. Again, it is quite clear that the enriched scheme (albeit enriched on a single element) performs much better than the non-enriched scheme. However, it is apparent that even the enriched scheme does not converge exponentially. By only enriching the element containing the singular point, there exist non-enriched elements which are `close' to the singularity. As such, the Sobolev seminorms will grow quite quickly as $k$ increases. Therefore, a large $k$ may be required before the boundedness (<ref>) is apparent. This explains why exponential convergence of the enriched scheme is not observed in Figure <ref>.

$E_{0,\Th}$ vs $k$

[ legend columns=3, legend to name=ktest ]
table[x=EdgeDegree,y=L2Error] data/ktest_hexa2_unenriched.dat;
table[x=EdgeDegree,y=L2Error] data/ktest_hexa2_enriched.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/test_file.dat;

Non-enriched, Enriched on a single element, $h^k$;

$E_{1,\Th}$ vs $k$

table[x=EdgeDegree,y=H1Error] data/ktest_hexa2_unenriched.dat;
table[x=EdgeDegree,y=H1Error] data/ktest_hexa2_enriched.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/test_file.dat;

$E_{\a,\xh}$ vs $k$

table[x=EdgeDegree,y=EnergyError] data/ktest_hexa2_unenriched.dat;
table[x=EdgeDegree,y=EnergyError] data/ktest_hexa2_enriched.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/test_file.dat;

Error vs $k$, $h\approx 0.195$

In Figure <ref> we consider two local enrichment schemes on a uniform Cartesian mesh with 48 elements. The elements which are enriched with the singular functions are bordered in red. Due to the larger element diameters and `rounder' element geometries we are able to consider a much larger radius of enrichment without ill-conditioning causing the scheme to fail for large $k$. As such, the non-enriched elements are further from the singular point (particularly in scheme 2) so we expect the high-order Sobolev seminorms to be smaller. We plot convergence of the schemes with respect to $k$ in Figure <ref>. It appears that the locally enriched scheme 2 maintains exponential convergence with respect to $k$ in both $L^2$ and energy error. This is similar to results in enriched NCVEM <cit.> where exponential convergence of a globally enriched scheme on a coarse Voronoi mesh is observed. While the $H^1$ error of scheme 2 does observe a minor saturation in convergence rate, the rate is consistently better than $h^k$.  It appears that the convergence rate of scheme 1 becomes slightly sub-optimal as $k$ gets large due to the existence of non-enriched elements closer to the singular point.

Local enrichment scheme 1.

Local enrichment scheme 2.

Uniform Cartesian mesh with enriched elements bordered in red

$E_{0,\Th}$ vs $k$

[ legend columns=4, legend to name=ktest.cartesian, xtick=0,1,2,3,4,5,6,7, xticklabels=0,1,2,3,4,5,6,7 ]
table[x=EdgeDegree,y=L2Error] data/ktest_cartesian_2_nonenriched.dat;
table[x=EdgeDegree,y=L2Error] data/ktest_cartesian_2_locally_enriched_1.dat;

table[x=EdgeDegree,y=L2Error] data/ktest_cartesian_2_locally_enriched_2.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/logh_cartesian_2.dat;
Non-enriched, Locally enriched (scheme 1), Locally enriched (scheme 2), $h^k$;

$E_{1,\Th}$ vs $k$

[xtick=0,1,2,3,4,5,6,7, xticklabels=0,1,2,3,4,5,6,7]
table[x=EdgeDegree,y=H1Error] data/ktest_cartesian_2_nonenriched.dat;
table[x=EdgeDegree,y=H1Error] data/ktest_cartesian_2_locally_enriched_1.dat;

table[x=EdgeDegree,y=H1Error] data/ktest_cartesian_2_locally_enriched_2.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/logh_cartesian_2.dat;

$E_{\a,\xh}$ vs $k$

[xtick=0,1,2,3,4,5,6,7, xticklabels=0,1,2,3,4,5,6,7]
table[x=EdgeDegree,y=EnergyError] data/ktest_cartesian_2_nonenriched.dat;
table[x=EdgeDegree,y=EnergyError] data/ktest_cartesian_2_locally_enriched_1.dat;

table[x=EdgeDegree,y=EnergyError] data/ktest_cartesian_2_locally_enriched_2.dat;
[mark=none,densely dotted,blue] table[x=EdgeDegree,y=LogH_slope] data/logh_cartesian_2.dat;

Error vs $k$, Cartesian mesh, $h\approx 0.354$

§.§ Oscillatory Solution

We briefly consider here an application of the XHHO scheme to highly oscillatory solutions of problem (<ref>). If the source term consists of a highly oscillatory component, it is expected that the solution will reflect this. In standard HHO, such details can only be captured with either high-order approximations or fine meshes. The XHHO method, however, is able to obtain improved error estimates by enriching the local spaces with an oscillatory function.

Let us consider here a square domain $\Omega=(0,1)^2$ and exact solution

\begin{equation}\label{eq:oscillatory.solution}
u = \sin(\pi x_1) \sin(\pi x_2) + \sin(\frac{1}{\hat{r}^2 + \epsilon})
\end{equation}

with $\hat{r}^2=(x_1 - 0.5)^2 + (x_2 - 0.5)^2$ denoting the radial distance to the domain center and $\epsilon > 0$ a small constant. We note that equation (<ref>) does not actually define a singular solution, thus the standard HHO method is expected to converge optimally. However, as the multiplicative constants in the error depend on the $(k+2)$-th derivative of the solution, the error is expected to get very large as $\epsilon \to 0$.

The oscillatory component of (<ref>) is denoted by

\[
\psi \defeq \sin(\frac{1}{\hat{r}^2 + \epsilon}).
\]

The enrichment space is naturally defined as $W(T) = \textrm{span}\{\psi|_T\}$. We note that as $\Delta \psi \ne 0$ we are no longer able to consider $\POLYD{k}(T) = \POLY{k}(T)$. We take $\epsilon = 0.05$ and consider both non-enriched and locally enriched schemes (here, the local enrichment process is identical to that described in Section <ref>, except with the singular point taken as $\hat{r}=0$ rather than $r=0$). We consider a sequence of regular, triangular meshes and plot the results in Figures <ref> and <ref> for $k=0,1$. It is clear that while both schemes converge optimally, the absolute error of the non-enriched scheme can be several orders of magnitude worse than that of the locally enriched scheme.

$E_{0,\Th}$ vs $h$

[ legend columns=2, legend to name=oscillator_k0 ]
table[x=MeshSize,y=L2Error] data/oscillator_tri_mesh_unenriched_k0.dat;
table[x=MeshSize,y=L2Error] data/oscillator_tri_mesh_cut_05_k0.dat;
Non-enriched, Locally enriched ($\gamma=0.5$);

$E_{1,\Th}$ vs $h$

table[x=MeshSize,y=H1Error] data/oscillator_tri_mesh_unenriched_k0.dat;
table[x=MeshSize,y=H1Error] data/oscillator_tri_mesh_cut_05_k0.dat;

$E_{\a,\xh}$ vs $h$

table[x=MeshSize,y=EnergyError] data/oscillator_tri_mesh_unenriched_k0.dat;
table[x=MeshSize,y=EnergyError] data/oscillator_tri_mesh_cut_05_k0.dat;

Error vs $h$, $k=0$

$E_{0,\Th}$ vs $h$

[ legend columns=2, legend to name=oscillator_k1 ]
table[x=MeshSize,y=L2Error] data/oscillator_tri_mesh_unenriched_k1.dat;
table[x=MeshSize,y=L2Error] data/oscillator_tri_mesh_cut_05_k1.dat;
Non-enriched, Locally enriched ($\gamma=0.5$);

$E_{1,\Th}$ vs $h$

table[x=MeshSize,y=H1Error] data/oscillator_tri_mesh_unenriched_k1.dat;
table[x=MeshSize,y=H1Error] data/oscillator_tri_mesh_cut_05_k1.dat;

$E_{\a,\xh}$ vs $h$

table[x=MeshSize,y=EnergyError] data/oscillator_tri_mesh_unenriched_k1.dat;
table[x=MeshSize,y=EnergyError] data/oscillator_tri_mesh_cut_05_k1.dat;

Error vs $h$, $k=1$

§ CONCLUSION

In this paper we have introduced and analysed the Extended Hybrid High-Order method for the Poisson problem. The method is shown to be robust, and capable of handling generic singularities satisfying Assumption <ref> in arbitrary dimensions. Optimal error estimates are established in both $H^1$- and discrete energy norms. The error analysis is backed up by numerical simulations which show the Extended Hybrid High-Order method to be a viable technique for handling irregular solutions. The method is also applicable to any situation where some aspects of the solution can be (at least locally) determined. This is seen in Section <ref> where the XHHO method is applied to solutions possessing a highly oscillatory component. Indeed, the enrichment function considered is actually smooth, yet a considerable improvement in performance was observed which highlights the versatility of the method for any cases where there is prior knowledge of the solution behaviour. The scheme also has natural generalisations to more general linear elliptic problems of the form $Lu=f$.

The case of corner singularities is discussed and tested in detail in Section <ref>. The schemes are shown to converge optimally and perform significantly better than standard HHO for solutions possessing a weak singularity at a re-entrant corner in two dimensions. The XHHO method is capable of handling such irregularities for all but the limiting case of the slit domain. For such cases, further work into developing a robust Hybrid High-Order method is required. In three dimensions, the singularities arising from irregular geometries can be more complicated. Along an edge which forms a non-convex corner in the domain, the singularity can be written as the product of a smooth function and a weakly singular function whose behaviour perpendicular to the edge is described by the two-dimensional case (c.f. [28]). Again, this scheme is capable of modelling such irregularity for all but the slit domain.

§ DECLARATIONS

The author declares that they have no conflict of interest.

[1]
I. Aavatsmark, G. T. Eigestad, B. T. Mallison, and J. M. Nordbotten.
A compact multipoint flux approximation method with improved
Numer. Methods Partial Differential Equations,
24(5):1329–1360, 2008.

[2]
Bashir Ahmad, Ahmed Alsaedi, Franco Brezzi, L. Donatella Marini, and Alessandro
Equivalent projectors for virtual element methods.
Computers & Mathematics with Applications, 66(3):376–391,

[3]
Edoardo Artioli and Lorenzo Mascotto.
Enrichment of the nonconforming virtual element method with singular
Computer Methods in Applied Mechanics and Engineering,
385:114024, 2021.

[4]
Blanca Ayuso de Dios, Konstantin Lipnikov, and Gianmarco Manzini.
The nonconforming virtual element method.
ESAIM: Mathematical Modelling and Numerical Analysis,
50(3):879–904, 2016.

[5]
Ivo Babuška and Jens M Melenk.
The partition of unity method.
International journal for numerical methods in engineering,
40(4):727–758, 1997.

[6]
Laurenço Beirão da Veiga, Franco Brezzi, Andrea Cangiani, Gianmarco
Manzini, L Donatella Marini, and Alessandro Russo.
Basic principles of virtual element methods.
Mathematical Models and Methods in Applied Sciences,
23(01):199–214, 2013.

[7]
Ted Belytschko and Tom Black.
Elastic crack growth in finite elements with minimal remeshing.
International journal for numerical methods in engineering,
45(5):601–620, 1999.

[8]
E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar.
Extended virtual element method for the laplace problem with
singularities and discontinuities.
Computer Methods in Applied Mechanics and Engineering,
356:571–597, 2019.

[9]
Lorenzo Botti, Daniele Antonio Di Pietro, and Jérôme Droniou.
A Hybrid High-Order method for the incompressible Navier–Stokes
equations based on Temam's device.
J. Comput. Phys., 376:786–816, 2019.

[10]
Michele Botti, Daniele Antonio Di Pietro, and Pierre Sochala.
A hybrid high-order method for nonlinear elasticity.
SIAM Journal on Numerical Analysis, 55(6):2687–2717, 2017.

[11]
Franco Brezzi, Richard S. Falk, and L. Donatella Marini.
Basic principles of mixed virtual element methods.
ESAIM Math. Model. Numer. Anal., 48(4):1227–1240, 2014.

[12]
Andrea Cangiani, Zhaonan Dong, Emmanuil H. Georgoulis, and Paul Houston.
$hp$-version discontinuous Galerkin methods on polygonal and
polyhedral meshes.
SpringerBriefs in Mathematics. Springer, Cham, 2017.

[13]
Andrea Cangiani, Gianmarco Manzini, and Oliver J. Sutton.
Conforming and nonconforming virtual element methods for elliptic
IMA J. Numer. Anal., 37(3):1317–1354, 2017.

[14]
Eric B Chin, Jean B Lasserre, and N Sukumar.
Numerical integration of homogeneous functions on convex and
nonconvex polygons and polyhedra.
Computational Mechanics, 56(6):967–981, 2015.

[15]
B. Cockburn, B. Dong, J. Guzmán, M. Restelli, and R. Sacco.
A hybridizable discontinuous Galerkin method for steady-state
convection-diffusion-reaction problems.
SIAM J. Sci. Comput., 31(5):3827–3846, 2009.

[16]
Daniele Antonio Di Pietro and Jérôme Droniou.
A Hybrid High-Order method for Leray–Lions elliptic equations on
general meshes.
Mathematics of Computation, 86(307):2159–2191, 2017.

[17]
Daniele Antonio Di Pietro and Jérôme Droniou.
A third strang lemma and an aubin–nitsche trick for schemes in fully
discrete formulation.
Calcolo, 55(3):40, 2018.

[18]
Daniele Antonio Di Pietro and Jérôme Droniou.
The Hybrid High-Order Method for Polytopal Meshes: Design,
Analysis, and Applications, volume 19 of Modeling, Simulation and
Springer International Publishing,
<https://hal.archives-ouvertes.fr/hal-02151813>, 01 2020.

[19]
Daniele Antonio Di Pietro, Jérôme Droniou, and Alexandre Ern.
A discontinuous-skeletal method for advection-diffusion-reaction on
general meshes.
SIAM Journal on Numerical Analysis, 53, 11 2014.

[20]
Daniele Antonio Di Pietro and Alexandre Ern.
Mathematical aspects of discontinuous Galerkin methods,
volume 69.
Springer Science & Business Media, 2011.

[21]
Daniele Antonio Di Pietro and Alexandre Ern.
A hybrid high-order locking-free method for linear elasticity on
general meshes.
Computer Methods in Applied Mechanics and Engineering,
283:1–21, 2015.

[22]
Daniele Antonio Di Pietro, Alexandre Ern, and Simon Lemaire.
An arbitrary-order and compact-stencil discretization of diffusion on
general meshes based on local reconstruction operators.
Computational Methods in Applied Mathematics, 14(4):461–472,

[23]
Daniele Antonio Di Pietro, Alexandre Ern, Alexander Linke, and Friedhelm
A discontinuous skeletal method for the viscosity-dependent stokes
Computer Methods in Applied Mechanics and Engineering,
306:175–195, 2016.

[24]
Jérôme Droniou.
Interplay between diffusion anisotropy and mesh skewness in hybrid
high-order schemes.
In International Conference on Finite Volumes for Complex
Applications, pages 3–23. Springer, 2020.

[25]
Jérôme Droniou, Robert Eymard, Thierry Gallouët, and Raphaèle
A unified approach to mimetic finite difference, hybrid finite volume
and mixed finite volume methods.
Mathematical Models and Methods in Applied Sciences,
20(02):265–295, 2010.

[26]
Jérôme Droniou and Liam Yemm.
Robust hybrid high-order method on polytopal meshes with small faces.
Computational Methods in Applied Mathematics, 22(1):47–71,

[27]
Pierre Grisvard.
Elliptic Problems in Nonsmooth Domains.
Pitman Advanced Publishing Program, 1985.

[28]
Pierre Grisvard.
Singularities in boundary value problems and exact controllability of
hyperbolic systems.
In Optimization, Optimal Control and Partial Differential
Equations, pages 77–84. Springer, 1992.

[29]
Jens M Melenk and Ivo Babuška.
The partition of unity finite element method: basic theory and
Computer methods in applied mechanics and engineering,
139(1-4):289–314, 1996.

[30]
Nicolas Moës, John Dolbow, and Ted Belytschko.
A finite element method for crack growth without remeshing.
International journal for numerical methods in engineering,
46(1):131–150, 1999.

[31]
Lin Mu, Junping Wang, and Xiu Ye.
Weak Galerkin finite element methods on polytopal meshes.
Int. J. Numer. Anal. Model., 12(1):31–53, 2015.

[32]
N. Sukumar and A. Tabarraei.
Conforming polygonal finite elements.
Internat. J. Numer. Methods Engrg., 61(12):2045–2066, 2004.

[33]
Arash Zamani and M Reza Eslami.
Embedded interfaces by polytope fem.
International journal for numerical methods in engineering,
88(8):715–748, 2011.

[1]
I. Aavatsmark, G. T. Eigestad, B. T. Mallison, and J. M. Nordbotten.
A compact multipoint flux approximation method with improved
Numer. Methods Partial Differential Equations,
24(5):1329–1360, 2008.

[2]
Bashir Ahmad, Ahmed Alsaedi, Franco Brezzi, L. Donatella Marini, and Alessandro
Equivalent projectors for virtual element methods.
Computers & Mathematics with Applications, 66(3):376–391,

[3]
Edoardo Artioli and Lorenzo Mascotto.
Enrichment of the nonconforming virtual element method with singular
Computer Methods in Applied Mechanics and Engineering,
385:114024, 2021.

[4]
Blanca Ayuso de Dios, Konstantin Lipnikov, and Gianmarco Manzini.
The nonconforming virtual element method.
ESAIM: Mathematical Modelling and Numerical Analysis,
50(3):879–904, 2016.

[5]
Ivo Babuška and Jens M Melenk.
The partition of unity method.
International journal for numerical methods in engineering,
40(4):727–758, 1997.

[6]
Laurenço Beirão da Veiga, Franco Brezzi, Andrea Cangiani, Gianmarco
Manzini, L Donatella Marini, and Alessandro Russo.
Basic principles of virtual element methods.
Mathematical Models and Methods in Applied Sciences,
23(01):199–214, 2013.

[7]
Ted Belytschko and Tom Black.
Elastic crack growth in finite elements with minimal remeshing.
International journal for numerical methods in engineering,
45(5):601–620, 1999.

[8]
E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar.
Extended virtual element method for the laplace problem with
singularities and discontinuities.
Computer Methods in Applied Mechanics and Engineering,
356:571–597, 2019.

[9]
Lorenzo Botti, Daniele Antonio Di Pietro, and Jérôme Droniou.
A Hybrid High-Order method for the incompressible Navier–Stokes
equations based on Temam's device.
J. Comput. Phys., 376:786–816, 2019.

[10]
Michele Botti, Daniele Antonio Di Pietro, and Pierre Sochala.
A hybrid high-order method for nonlinear elasticity.
SIAM Journal on Numerical Analysis, 55(6):2687–2717, 2017.

[11]
Franco Brezzi, Richard S. Falk, and L. Donatella Marini.
Basic principles of mixed virtual element methods.
ESAIM Math. Model. Numer. Anal., 48(4):1227–1240, 2014.

[12]
Andrea Cangiani, Zhaonan Dong, Emmanuil H. Georgoulis, and Paul Houston.
$hp$-version discontinuous Galerkin methods on polygonal and
polyhedral meshes.
SpringerBriefs in Mathematics. Springer, Cham, 2017.

[13]
Andrea Cangiani, Gianmarco Manzini, and Oliver J. Sutton.
Conforming and nonconforming virtual element methods for elliptic
IMA J. Numer. Anal., 37(3):1317–1354, 2017.

[14]
Eric B Chin, Jean B Lasserre, and N Sukumar.
Numerical integration of homogeneous functions on convex and
nonconvex polygons and polyhedra.
Computational Mechanics, 56(6):967–981, 2015.

[15]
B. Cockburn, B. Dong, J. Guzmán, M. Restelli, and R. Sacco.
A hybridizable discontinuous Galerkin method for steady-state
convection-diffusion-reaction problems.
SIAM J. Sci. Comput., 31(5):3827–3846, 2009.

[16]
Daniele Antonio Di Pietro and Jérôme Droniou.
A Hybrid High-Order method for Leray–Lions elliptic equations on
general meshes.
Mathematics of Computation, 86(307):2159–2191, 2017.

[17]
Daniele Antonio Di Pietro and Jérôme Droniou.
A third strang lemma and an aubin–nitsche trick for schemes in fully
discrete formulation.
Calcolo, 55(3):40, 2018.

[18]
Daniele Antonio Di Pietro and Jérôme Droniou.
The Hybrid High-Order Method for Polytopal Meshes: Design,
Analysis, and Applications, volume 19 of Modeling, Simulation and
Springer International Publishing,
<https://hal.archives-ouvertes.fr/hal-02151813>, 01 2020.

[19]
Daniele Antonio Di Pietro, Jérôme Droniou, and Alexandre Ern.
A discontinuous-skeletal method for advection-diffusion-reaction on
general meshes.
SIAM Journal on Numerical Analysis, 53, 11 2014.

[20]
Daniele Antonio Di Pietro and Alexandre Ern.
Mathematical aspects of discontinuous Galerkin methods,
volume 69.
Springer Science & Business Media, 2011.

[21]
Daniele Antonio Di Pietro and Alexandre Ern.
A hybrid high-order locking-free method for linear elasticity on
general meshes.
Computer Methods in Applied Mechanics and Engineering,
283:1–21, 2015.

[22]
Daniele Antonio Di Pietro, Alexandre Ern, and Simon Lemaire.
An arbitrary-order and compact-stencil discretization of diffusion on
general meshes based on local reconstruction operators.
Computational Methods in Applied Mathematics, 14(4):461–472,

[23]
Daniele Antonio Di Pietro, Alexandre Ern, Alexander Linke, and Friedhelm
A discontinuous skeletal method for the viscosity-dependent stokes
Computer Methods in Applied Mechanics and Engineering,
306:175–195, 2016.

[24]
Jérôme Droniou.
Interplay between diffusion anisotropy and mesh skewness in hybrid
high-order schemes.
In International Conference on Finite Volumes for Complex
Applications, pages 3–23. Springer, 2020.

[25]
Jérôme Droniou, Robert Eymard, Thierry Gallouët, and Raphaèle
A unified approach to mimetic finite difference, hybrid finite volume
and mixed finite volume methods.
Mathematical Models and Methods in Applied Sciences,
20(02):265–295, 2010.

[26]
Jérôme Droniou and Liam Yemm.
Robust hybrid high-order method on polytopal meshes with small faces.
Computational Methods in Applied Mathematics, 22(1):47–71,

[27]
Pierre Grisvard.
Elliptic Problems in Nonsmooth Domains.
Pitman Advanced Publishing Program, 1985.

[28]
Pierre Grisvard.
Singularities in boundary value problems and exact controllability of
hyperbolic systems.
In Optimization, Optimal Control and Partial Differential
Equations, pages 77–84. Springer, 1992.

[29]
Jens M Melenk and Ivo Babuška.
The partition of unity finite element method: basic theory and
Computer methods in applied mechanics and engineering,
139(1-4):289–314, 1996.

[30]
Nicolas Moës, John Dolbow, and Ted Belytschko.
A finite element method for crack growth without remeshing.
International journal for numerical methods in engineering,
46(1):131–150, 1999.

[31]
Lin Mu, Junping Wang, and Xiu Ye.
Weak Galerkin finite element methods on polytopal meshes.
Int. J. Numer. Anal. Model., 12(1):31–53, 2015.

[32]
N. Sukumar and A. Tabarraei.
Conforming polygonal finite elements.
Internat. J. Numer. Methods Engrg., 61(12):2045–2066, 2004.

[33]
Arash Zamani and M Reza Eslami.
Embedded interfaces by polytope fem.
International journal for numerical methods in engineering,
88(8):715–748, 2011.