11institutetext:  University of Texas at Austin
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# You Are What You Tweet: Profiling Users by Past Tweets to Improve Hate
Speech Detection

Prateek Chaudhry and Matthew Lease

###### Abstract

Hate speech detection research has predominantly focused on purely content-
based methods, without exploiting other contextual data. We briefly critique
pros and cons of this task formulation. We then investigate profiling users by
their past utterances as an informative prior to better predict whether new
utterances constitute hate speech. To evaluate this, we augment three Twitter
hate speech datasets with additional timeline data, then embed this additional
context into a strong baseline model. Promising results suggest merit for
further investigation.

###### Keywords:

hate speech classification modeling profiles twitter

##  1 Introduction

Online hate speech is a vast and continually growing problem [10, 12, 13, 18].
The detection task is most commonly framed as purely content-based: each
utterance is classified without any additional context. However, history often
repeats itself, and “A user who is known to write hate speech messages may do
so again” [18]. Facebook researchers have also recently noted that “a user or
group that has posted violating content in the past may be prone to do so more
often in the future” [10]. This suggests that modeling of user priors may be
highly informative and complementary to purely content-based detection models.
We return to discussion of ethical considerations in Section 7.

To investigate user profiling, we augment three Twitter hate speech datasets
[1, 5, 20] with additional recent timeline Tweets by each author, then embed
this additional context into a strong baseline model [2] that was later
further refined [1]. Results show strong improvement on one dataset, little
benefit on another, and fairly consistent improvement on a third. Overall,
results across experimental conditions and metrics suggest user modeling
merits further work, though analysis is complicated by differences in
annotation schemes and processes, as well as Twitter API limitations and data
sharing policies.

##  2 Related Work

While most hate speech detection models have been content-oriented and non-
contextual, there are notable exceptions. After Hovy [11] showed modeling
value of demographics for classification tasks, Waseem and Hovy [20] applied
this to hate speech, inferring gender (by name) and location (by timezone) to
classify hate speech with modest benefit: only gender raised hate detection
accuracy, and it was statistically significant only when both gender and
location were used. A challenge with demographics is that they are sensitive
and often not widely available, if captured at all by the platform.

Rather than classify hateful utterances, other work has sought to classify
hateful users. Assuming “birds of a feather flock together”, community
detection [8] has sought to identify hater communities. Ribeiro et al. [17]
and Mathew et al. [14] analyze Twitter’s retweet graph to detect users likely
to spread hateful content. Mishra et al. [15] is the only work we are aware of
to profile hateful users, via the follower graph, and then use these profiles
to improve hate speech detection over baseline models.

Some hate speech datasets are heavily dominated by a few prolific haters. In
Waseem and Hovy [20], all racist-labeled tweets came from 9 users. In fact,
Arango et al. [1] note that a single user generates 96% of all racist tweets,
while another user produces 44% of all sexist tweets. Arango et al. aptly
critique that benchmarking on such highly skewed datasets risks overfitting a
few individuals rather than learning a broadly applicable model. Similar
concerns have been voiced about over-fitting idiosyncrasies of a few
annotators [9]. We agree that datasets should have many diverse examples of
the phenomena being modeled (i.e., hate speech examples from many users).
However, dataset composition is orthogonal to model design: how can we best
model historical context in prediction and fairly evaluate vs. context-free
models?

Some past studies (cf., [20, 21]) have partitioned train/test data by
utterance (i.e., by tweet), rather than by user. The consequence of such
experimental design is that any model trained and tested in this manner can be
expected to make better predictions on users found in both train and test
splits. However, we should not confuse desirable user modeling with
questionable experimental design in which user history is only captured
haphazardly by whichever users happen to have some tweets in the training
data, and the ratio of that user-specific history vs. training tweets from
other users. Such a scheme does not reflect intentional design or controlled
evaluation of user modeling.

An intuitive language-based approach is to model each user by a “document” of
their past utterances. Such historical data is certainly available to social
media companies, and often publicly via API or crawling. The closest work we
are aware of on content-based modeling of haters is by Dadvar et al. [4],
though they only considered the user’s past number of obscene words as a
modeling feature.

##  3 Datasets

Figure 1: Log-log distribution of tweets per user (top) and hate tweets per
user for users producing 1 or more hate tweets (bottom), for all three
datasets. The y-axis indicates the number of tweets while the x-axis simply
enumerates users in descending order of activity.

We adopt the same three datasets studied by Arango et al. [1]; see statistics
in Table 1. As is common in other problem domains, each hate speech dataset
has various limitations; by evaluating across several datasets, we test across
varying data conditions.

Waseem and Hovy [20] ( W&H ) label 16K tweets for 3 classes: Racism, Sexism
and Neither. After author labeling, tweets labeled as hate speech were
reviewed by a 25-year old female “studying gender studies and a nonactivist
feminist”. Agreement between author and reviewer was $\kappa$= 0.84, with 85%
of disagreements on sexist labels, and 98% of these changed to neither. Tweet
IDs and labels were shared, with tweets obtained via the Twitter API. Since
tweets and user accounts are often deleted, 15K tweets are found in [1]. We
find 10K tweets (and only 6 racist ones) over 1,458 users.

Davidson et al. [5] crowd-label 25K tweets for 3 classes: Hate, Offensive and
Neither. Majority voting over 3 or more workers was used for label
aggregation. 200 tweets were discarded with no majority over the 3 classes.
For hate, only 5% were labeled by 2/3 and only 1.3% by 3/3. Most were labeled
offensive (76% at 2/3, 53% at 3/3) and the rest non-offensive (16.6% at 2/3,
11.8% at 3/3). We find 18K tweets ($\sim$75%) from 741 users.

Arango et al. [1] identify three issues: 1) W&H has few haters; 2) those
haters are prolific; and 3) train/test splits by utterance mean the same hater
often appears in both train and test splits. To address this, they 1) add to
W&H all hate-labeled tweets from DAVIDSON; 2) restrict to at most 250 tweets
per user per class; and 3) perform train/test split by user rather than by
utterance. Their fused dataset had 7K tweets, collapsing labels to simple
binary hate vs. non-hate. We find 5.4K over 1746 users.

Figure 1 shows the distribution of tweets per user across datasets. For all
users who produce 1 or more hate tweet (W&H: 443 users, Davidson: 566, Arango:
734), we also plot the number of hate tweets from each such “hater”.

##  4 Profiling Users by Past Tweets

As a strong baseline, we adopt Badjatiya et al. [2]’s 2-phase model, using
[1]’s corrected version in which word embeddings are derived only from
training data. An LSTM classifier is first trained to predict the label. Each
word is converted to a dense vector representation using a word embedding
matrix initialized with pretrained GloVe embeddings [16]. The final vector
from the LSTM is followed by a fully connected layer and a softmax or sigmoid
layer for obtaining prediction probability. The network is trained with cross
entropy loss with Adam optimizer. Once training of LSTM is done, the first
layer is extracted, i.e. word embedding matrix fine tuned on the training set.
In the second phase, a tweet is converted to a fixed sized vector by averaging
the embedding vectors of its tokens using the trained embeddings from the
previous phase. Representing the tweet by this vector, a Gradient Boosted
Decision Tree (GBDT) is trained for classification.

Given a tweet ID, we use the Twitter API to retrieve the author’s
timeline111https://developer.twitter.com/en/docs/tweets/timelines/: their
latest 20 tweets. This size of 20 reflects an API maximum, but it would be
interesting in future work to model longer histories. Given this history, we
augment existing public datasets (section 3) with these timelines. W&H
includes tweet IDs, and Davidson kindly shared their tweet IDs with us upon
request. In practice, timeline tweets should precede tweets being classified;
here we use existing datasets. By augmenting existing datasets, we can assess
the relative benefit with and without user profiling on hate speech datasets
that are already familiar to the research community.

We utilize these user profiles in (Arango et al. [1]’s corrected version of)
[2]’s model as follows. Given a tweet, the author’s timeline is averaged using
the same trained word embedding used for tweet representation. Following this,
both tweet representation and timeline representation are concatenated and
used to train the GBDT classifier. In general, a history-informed prior of the
user ought to complement any content-based model of the tweet alone. We expect
future work will continue to benefit from exploiting such a prior-informed
modeling architecture, irrespective of the specific model.

##  5 Experimental Setup

We build on Badjatiya et al. [2]’s and Arango et al. [1]’s shared source code.
We use 10-fold cross validation with default hyperparameters for [2]’s model.
The LSTM word embeddings are initialized with 200 dimensional pretrained GloVe
embeddings, and the size of the LSTM representation is 200. We also add
dropout of 0.25 and 0.5 after the word embedding layer and the LSTM,
respectively. We train the LSTM based architecture for 10 epochs, and then
train the GBDT using the final word embeddings.

With nearly all racism tweets in W&H deleted, we train it as a binary
classifier: sexism vs. none. For DAVIDSON, we train a ternary classifier over
its three classes. ARANGO fuses W&H and DAVIDSON label sets by collapsing
classes into simple binary classification of hate vs. non-hate.

Tweet deletions make comparison to prior published results more difficult,
with only 2/3 of W&H and 3/4 of DAVIDSON datasets still available, and nearly
all W&H racist tweets deleted. Arango et al. [1] note that partitioning
train/test by tweet, as in past studies, results in prolific tweeters
appearing in both splits. The risk this poses is potentially overfitting to
particular users. Instead, they argue for splitting train/test by user. To be
as comprehensive as possible, we report results both ways for completeness.

##  6 Results

| Split by Tweet | Split by User
---|---|---
| Baseline | With Timeline | Baseline | With Timeline
Dataset | Label | Tweets | Found | P | R | F1 | P | R | F1 | P | R | F1 | P | R | F1
W&H | Racism | 1,851 | 6 | - | - | - | - | - | - | - | - | - | - | - | -
Sexism | 2,988 | 2,777 | 76.8 | 65.6 | 70.7 | 83.7 | 80.2 | 81.9 | 59.2 | 33.2 | 37.8 | 63.7 | 35.4 | 40.8
Neither | 10,110 | 7,320 | 87.7 | 92.4 | 89.9 | 92.6 | 94.1 | 93.3 | 82.2 | 91.4 | 84.3 | 82.6 | 93.9 | 86.5
Micro Avg |  |  | 85.1 | 85.1 | 85.1 | 90.3 | 90.3 | 90.3 | 76.0 | 76.0 | 76.0 | 79.3 | 79.3 | 79.3
| Macro Avg | 14,949 | 10,103 | 82.2 | 79.0 | 80.3 | 88.2 | 87.2 | 87.6 | 70.7 | 62.3 | 61.1 | 73.1 | 64.6 | 63.6
DAVID-SON | Hate | 1,430 | 1,118 | 55.4 | 24.3 | 33.6 | 56.4 | 24.7 | 34.1 | 53.6 | 22.1 | 31.0 | 52.2 | 21.5 | 30.3
Offensive | 19,190 | 13,931 | 91.3 | 95.9 | 93.5 | 91.2 | 95.9 | 93.5 | 91.2 | 95.9 | 93.5 | 91.2 | 95.8 | 93.4
Neither | 4,163 | 2,978 | 83.8 | 81.5 | 82.7 | 83.7 | 81.1 | 82.3 | 82.8 | 81.3 | 82.0 | 82.7 | 81.2 | 81.9
Micro Avg |  |  | 89.1 | 89.1 | 89.1 | 89.0 | 89.0 | 89.0 | 88.9 | 88.9 | 88.9 | 88.8 | 88.8 | 88.8
| Macro Avg | 24,783 | 18,027 | 76.8 | 67.3 | 69.9 | 77.1 | 67.2 | 70.0 | 75.9 | 66.4 | 68.8 | 75.3 | 66.2 | 68.6
ARANGO | Hate | 2,920 | 1,988 | 81.0 | 72.4 | 76.4 | 85.2 | 75.2 | 79.8 | 80.7 | 63.1 | 68.5 | 79.6 | 65.7 | 69.8
None | 4,086 | 3,481 | 85.1 | 90.3 | 87.6 | 86.7 | 92.5 | 89.5 | 81.5 | 87.9 | 82.8 | 82.6 | 86.1 | 82.3
Micro Avg |  |  | 83.8 | 83.8 | 83.8 | 86.1 | 86.1 | 86.1 | 79.8 | 79.8 | 79.8 | 80.1 | 80.1 | 80.1
Macro Avg | 7,006 | 5,469 | 83.1 | 81.3 | 82.0 | 86.0 | 83.8 | 84.6 | 81.1 | 75.5 | 75.7 | 81.1 | 75.9 | 76.1

Table 1: Baseline vs. timeline accuracy across datasets and train/test
partitions: by tweet or user. For each (dataset, split, metric) tuple, we bold
the higher accuracy between baseline vs. timeline.

| Baseline | With Timeline
---|---|---
Dataset | Size Timeline | Tweets | P | R | F | P | R | F
W&H | 0-5 | 58 | 78.2 | 74.1 | 75.8 | 68.8 | 67.4 | 68.0
6-10 | 18 | 47.1 | 47.1 | 47.1 | 75.0 | 97.1 | 81.8
11-15 | 34 | 60.0 | 54.7 | 53.0 | 50.5 | 50.2 | 46.0
16-20 | 9,991 | 82.2 | 79.0 | 80.4 | 88.3 | 87.4 | 87.9
DAVIDSON | 0-5 | 4,943 | 76.6 | 66.6 | 69.6 | 76.3 | 66.4 | 69.4
6-10 | 22 | 96.5 | 72.2 | 80.4 | NAN | 0 | NAN
11-15 | 287 | 66.0 | 58.8 | 61.5 | 58.5 | 54.0 | 55.8
16-20 | 12,775 | 76.9 | 67.5 | 70.1 | 77.2 | 67.6 | 70.3
ARANGO | 0-5 | 59 | 76.5 | 75.4 | 75.9 | 70.0 | 70.0 | 70.0
6-10 | 20 | 100.0 | 100.0 | 100.0 | 97.2 | 83.3 | 88.6
11-15 | 49 | 81.2 | 80.3 | 79.5 | 77.0 | 76.2 | 75.4
16-20 | 5,341 | 83.0 | 81.3 | 82.0 | 86.1 | 84.0 | 84.9
ARANGO | 0-5 | 59 | 81.9 | 77.8 | 79.4 | 85.3 | 78.9 | 81.3
(split by user) | 6-10 | 20 | 97.2 | 83.3 | 88.6 | 97.2 | 83.2 | 88.6
11-15 | 49 | 84.9 | 80.8 | 79.2 | 84.9 | 80.8 | 79.2
| 16-20 | 5,341 | 78.3 | 76.1 | 76.9 | 78.2 | 77.1 | 77.6

Table 2: Detection accuracy on user subgroups, based on amount of Twitter
timeline available per user (0-20 most recent tweets). Train/test partition of
data is by tweet unless noted otherwise.

Table 1 compares model performance with and without user timeline history
across the three datasets and the two train/test data partitions: by tweet vs.
by user. With regard to testing baseline vs. timeline results, splitting by
user is cleaner because for any users appearing in the test data, there are no
tweets from that user in the training set.

Results show strong improvement on W&H but little change on DAVIDSON. We do
see split-by-tweet shows modest improvement for DAVIDSON on the Hate category,
but not on other categories or for split-by-user. Between W&H and DAVIDSON
extremes, we see more modest but fairly consistent improvement on Arango et
al.’s dataset, which fuses W&H and DAVIDSON by adding to W&H all hate-labeled
tweets from DAVIDSON and down-sampling to at most 250 tweets per user per
class to reduce skew. Consistent with Arango et al. [1], we see much higher
results on W&H and their fused dataset when splitting by tweet rather than by
user (Arango et al. do not report on DAVIDSON).

Given the difference in model performance across the datasets, what explains
this? Table 1 and Figure 1 show important differences in data scale and
distribution across different datasets. The classes being annotated also
differ, as does the method of annotation: W&H uses traditional annotators
while DAVIDSON relies on crowd annotators. Such differences highlight several
dimensions of ongoing debate in hate speech research surrounding differing
approaches to annotation categories and processes [3, 7, 19].

We also conducted a further analysis to assess whether the amount of timeline
history per user varied significantly across datasets. Results appear in Table
2, where we bin users by the count of timeline tweets found: 0-5, 6-10, 11-15,
or 16-20. We see that the vast majority of users have 16-20 past timeline
tweets, and for this category we see consistent improvement across datasets
and train/test split conditions, driving scores.

Crucially, however, on DAVIDSON we see nearly 5K tweets ($\sim$28%) from users
with only 0-5 timeline tweets available, where we expect small to no
improvement from timelines due to lack of history. This holds across datasets
when splitting by tweet, including DAVIDSON, which indeed drives down overall
results on this dataset. When splitting by user on ARANGO, we do see increased
accuracy even with 0-5 timeline tweets, though the sample size is quite small
to draw conclusions.

Further limited qualitative analysis helped give us a sense of examples where
profiling a user created an informative prior for correctly classifying tweets
that might have been missed otherwise. We explored most prolific tweeters in
W&H for evidence of sexist timeline tweets. As an example, one user made a
possibly sexist tweet “.@USER1 @USER2 when was she good? i confuse her and ten
other women, which is why their pay is lower btw. supply vs. demand.” which
was easier to predict by their earlier timeline tweet “#FeminismIsAwful”. Use
of the timeline correctly predicts the ‘sexist’ label, and increases overall
detection accuracy for this user from 53% to 87%.

##  7 Ethical Considerations

While user modeling is common in personalized ranking and recommendation
systems (e.g., Google or Netflix), some users may wish to opt-out, and GDPR
allows “the right to be forgotten”. User modeling to prevent fraud and abuse,
while equally common in commercial systems, may raise different sorts of
ethical and legal questions. We know that social media platforms already do
monitor user accounts for terms of service violations, suspending accounts
demonstrating repeated abusive behaviors [10]. Given the vast scale and
expense of commercial content moderation, and extensive user histories
available, platforms cannot afford to ignore such a strong predictive signal.

However, what challenges may arise? Would a reformed hater have difficulty
overcoming his/her past profile a model had learned? One idea would be to
decay weight assigned to past tweets by age, or to completely restrict history
to a recent window (e.g., our current timeline of past 20 tweets). In fact,
the opposite problem is significant on platforms today: a bad actor who is
suspended often circumvents this by creating a new account, with a cycle of
continuing abusive behaviors by re-entrant bad actors.

##  8 Conclusion

Producers of hate speech are often repeat offenders, yet we know of no prior
work explicitly profiling users by their past tweets to improve hate speech
detection accuracy. We collect Twitter timeline data toward this end to
augment existing hate speech datasets. Results on several datasets, metrics,
and experimental settings are encouraging, but confounds remain. Future work
might explore better modeling (e.g., via BERT [6]), collecting more user
history, and combining with other profiling approaches [15]. Other approaches
may explore adding other features extracted from user twitter timeline, such
as hate speech tweet frequency, number of hate speech retweets, etc.

## Acknowledgments

We thank the reviewers for their valuable feedback, and the many talented
workers who provided the hate speech annotations that enable research on it.
This research was supported in part by Wipro, the Micron Foundation, and by
Good Systems222http://goodsystems.utexas.edu/, a UT Austin Grand Challenge to
develop responsible AI technologies. The statements made herein are solely the
opinions of the authors and do not reflect the views of the sponsoring
agencies.

## References

  * Arango et al. [2019] Arango, A., Pérez, J., Poblete, B.: Hate speech detection is not as easy as you may think: A closer look at model validation.  In: Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval. pp. 45–54 (2019), https://github.com/aymeam/User_distribution_experiments
  * Badjatiya et al. [2017] Badjatiya, P., Gupta, S., Gupta, M., Varma, V.: Deep learning for hate speech detection in tweets.  In: Proceedings of the 26th International Conference on World Wide Web Companion. pp. 759–760 (2017), https://github.com/pinkeshbadjatiya/twitter-hatespeech
  * Balayn and Bozzon [2019] Balayn, A., Bozzon, A.: Designing evaluations of machine learning models for subjective inference: The case of sentence toxicity.  arXiv preprint arXiv:1911.02471 (2019)
  * Dadvar et al. [2013] Dadvar, M., Trieschnigg, D., Ordelman, R., de Jong, F.: Improving cyberbullying detection with user context.  In: European Conference on Information Retrieval. pp. 693–696. Springer (2013)
  * Davidson et al. [2017] Davidson, T., Warmsley, D., Macy, M., Weber, I.: Automated hate speech detection and the problem of offensive language.  In: Eleventh international aaai conference on web and social media (2017), https://github.com/t-davidson/hate-speech-and-offensive-language
  * Devlin et al. [2018] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805 (2018)
  * Fortuna et al. [2020] Fortuna, P., Soler, J., Wanner, L.: Toxic, hateful, offensive or abusive? what are we really classifying? an empirical analysis of hate speech datasets.  In: Proceedings of The 12th Language Resources and Evaluation Conference. pp. 6786–6794 (2020)
  * Fortunato [2010] Fortunato, S.: Community detection in graphs.  Physics reports 486(3-5), 75–174 (2010)
  * Geva et al. [2019] Geva, M., Goldberg, Y., Berant, J.: Are we modeling the task or the annotator? an investigation of annotator bias in natural language understanding datasets.  In: Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP). pp. 1161–1166 (2019)
  * Halevy et al. [2020] Halevy, A., Ferrer, C.C., Ma, H., Ozertem, U., Pantel, P., Saeidi, M., Silvestri, F., Stoyanov, V.: Preserving integrity in online social networks.  arXiv preprint arXiv:2009.10311 (2020)
  * Hovy [2015] Hovy, D.: Demographic factors improve classification performance.  In: Proceedings of the 53rd annual meeting of the Association for Computational Linguistics and the 7th international joint conference on natural language processing (volume 1: Long papers). pp. 752–762 (2015)
  * Jurgens et al. [2019] Jurgens, D., Chandrasekharan, E., Hemphill, L.: A just and comprehensive strategy for using nlp to address online abuse.  arXiv preprint arXiv:1906.01738 (2019)
  * MacAvaney et al. [2019] MacAvaney, S., Yao, H.R., Yang, E., Russell, K., Goharian, N., Frieder, O.: Hate speech detection: Challenges and solutions.  PloS one 14(8) (2019)
  * Mathew et al. [2019] Mathew, B., Dutt, R., Goyal, P., Mukherjee, A.: Spread of hate speech in online social media.  In: Proceedings of the 10th ACM Conference on Web Science. pp. 173–182 (2019)
  * Mishra et al. [2018] Mishra, P., Del Tredici, M., Yannakoudakis, H., Shutova, E.: Author profiling for abuse detection.  In: Proceedings of the 27th International Conference on Computational Linguistics. pp. 1088–1098 (2018)
  * Pennington et al. [2014] Pennington, J., Socher, R., Manning, C.D.: Glove: Global vectors for word representation.  In: Empirical Methods in Natural Language Processing (EMNLP). pp. 1532–1543 (2014), http://www.aclweb.org/anthology/D14-1162
  * Ribeiro et al. [2018] Ribeiro, M.H., Calais, P.H., Santos, Y.A., Almeida, V.A., Meira Jr, W.: Characterizing and detecting hateful users on twitter.  In: Twelfth international AAAI conference on web and social media (2018)
  * Schmidt and Wiegand [2017] Schmidt, A., Wiegand, M.: A survey on hate speech detection using natural language processing.  In: Proceedings of the Fifth International Workshop on Natural Language Processing for Social Media. pp. 1–10 (2017)
  * Waseem [2016] Waseem, Z.: Are you a racist or am i seeing things? annotator influence on hate speech detection on twitter.  In: Proceedings of the first workshop on NLP and computational social science. pp. 138–142 (2016)
  * Waseem and Hovy [2016] Waseem, Z., Hovy, D.: Hateful symbols or hateful people? predictive features for hate speech detection on twitter.  In: Proceedings of the NAACL student research workshop. pp. 88–93 (2016), https://github.com/zeerakw/hatespeech
  * Xiang et al. [2012] Xiang, G., Fan, B., Wang, L., Hong, J., Rose, C.: Detecting offensive tweets via topical feature discovery over a large scale twitter corpus.  In: Proceedings of the 21st ACM international conference on Information and knowledge management. pp. 1980–1984 (2012)


# Robust interior penalty
discontinuous Galerkin methods

Zhaonan Dong  1) Inria, 2 rue Simone Iff, 75589 Paris, France, and 2) CERMICS,
Ecole des Ponts, 77455 Marne-la-Vallée 2, France<EMAIL_ADDRESS>and
Emmanuil H. Georgoulis  1) School of Computing and Mathematical Sciences,
University of Leicester, University Road, Leicester, LE1 7RH, United Kingdom,
2) Department of Mathematics, School of Applied Mathematical and Physical
Sciences, National Technical University of Athens, Zografou 15780, Greece, and
3) IACM-FORTH, Crete, Greece<EMAIL_ADDRESS>
###### Abstract.

A new variant of the IPDG method is presented, involving carefully constructed
weighted averages of the gradient of the approximate solution. The method is
shown to be robust even for the most extreme simultaneous local mesh,
polynomial degree and diffusion coefficient variation scenarios, without
resulting into unreasonably large penalization. The new IPDG method,
henceforth termed as _robust IPDG_ (RIPDG), offers typically significantly
better error behaviour and conditioning than the standard IPDG method when
applied to scenarios with strong mesh/polynomial degree/diffusion local
variation, especially when the underlying approximation space does not contain
a sufficiently rich conforming subspace. The latter is of particular
importance in the context of IPDG methods on polygonal/polyhedral meshes. On
the other hand, when using uniform meshes, constant polynomial degree for
problems with constant diffusion coefficients the RIPDG method is identical to
the classical IPDG. Numerical experiments indicate the favourable performance
of the new RIPDG method over the classical version in terms of conditioning
and error.

##  1\. Introduction

Discontinuous Galerkin (DG) methods are a staple in the cascade of Galerkin
approaches for the numerical solution of PDEs in divergence form. For the
discretization of elliptic operators, the class of interior penalty DG (IPDG)
methods is arguably the most popular choice due to its versatility and
simplicity in implementation. The consistency of the IPDG method is achieved
by the inclusion of specific terms involving integrals of solution fluxes on
element interfaces. The basic, yet revolutionary, idea in IPDG methods is the
further addition of penalization terms involving the jump of the approximate
solution (or its derivatives for higher order elliptic PDEs) to control the
departure from the PDE solution space and, therefore, to achieve stability in
a consistent fashion. The IPDG method design concept was presented in full by
Baker [4] and, briefly later, by Wheeler [28] and Arnold [2], thereby
providing a consistent extension to the classical work of Nitsche [23] for the
weak imposition of boundary conditions and to the classical inconsistent
finite element method with penalty of Babuška [3]. These classical IPDG
methods and their $hp$-version variants [25, 21] are widely used to this day.

Classical interior penalty discontinuous Galerkin (IPDG) methods for diffusion
problems require a number of assumptions on the local variation of mesh-size,
polynomial degree, and of the diffusion coefficient to determine the values of
the, so-called, discontinuity-penalization parameter and/or to perform error
analysis. Variants of IPDG methods involving weighted averages of the normal
flux functions of the approximate solution have been proposed [15, 6, 16] in
the context of high-contrast diffusion coefficients to mitigate the dependence
of the contrast in the stability and in the error analysis. More recently,
weighted averages involving local meshsizes have been proposed for treating
interface problems [1, 20] on non-matching grids on either side of the
interface.

To fix ideas, we shall focus on the basic second order elliptic boundary value
problem on an open domain $\Omega\subset\mathbb{R}^{d}$, $d\in\mathbb{N}$:
find $u\in H^{1}(\Omega)$ such that

(1.1) $\displaystyle-\nabla\cdot(a\nabla u)=$ $\displaystyle\ f\quad\text{ in
}\Omega;$ $\displaystyle u=$ $\displaystyle\ g\quad\text{ on }\partial\Omega,$

with data $f\in L_{2}(\Omega)$, $g\in L_{2}(\partial\Omega)$ and a positive
definite diffusion tensor $a\in[L_{\infty}(\Omega)]^{d\times d}$. Although
considerable extensions to the problem can be incorporated immediately to the
developments presented below, we refrain from exhausting the generalization
possibilities in the interest of clarity of exposition.

Classical IPDG methods for the problem (1.1) involve a user-defined quantity,
the so-called _penalty_ or _discontinuity-penalization_ parameter. The
judicious choice of this parameter is crucial for the stability of the IPDG
method, but also for its practical behaviour regarding approximation quality
and conditioning. In particular, it is well known that excessive penalty
parameter typically leads to increased ill-conditioning of the IPDG stiffness
matrix. At the same time, the penalty parameter has to be chosen large enough
to ensure stability of the method.

In this work, we focus on the in the behaviour of interior penalty methods in
the presence of highly locally variable meshes, local polynomial degrees,
and/or diffusion coefficients $a$. As we will see below, in such extreme,
often combined, local variation scenarios of discretization parameters and/or
coefficients, the classical IPDG method typically requires excessively large
penalization to retain stability. In fact, _it is possible that classical IPDG
penalty parameter choices degenerate to infinity while the dimension of the
approximation spaces remains finite_. This, in turn may be detrimental to the
conditioning of the problem and, in certain cases, to the quality of the
approximation.

To resolve this state of affairs, we present a new weighted-type interior
penalty discontinuous Galerkin method which is provably stable with a new
choice of the discontinuity-penalization parameter and is _robust_ with
respect to extreme local variations in the discretization parameters and in
the diffusion coefficient. The method, termed henceforth as _robust IPDG_
(RIPDG) replaces the classical average-of-the-normal-flux term(s) in the IPDG
by suitable, explicitly constructed weighted averages. This modification
results into a robust and concrete selection of the discontinuity-penalization
parameter, without resulting into excessive penalization in extreme combined
local mesh, polynomial degree and coefficient variation scenarios. We note
that a balanced choice for the discontinuity-penalization parameter is
particularly important when _no_ conforming subspace is available. Moreover,
the new RIPDG method allows for the first time in the interior penalty
literature for the proof of a priori error bounds _without_ any local-quasi-
uniformity/local bounded variation assumptions for the discretization
parameters.

The modifications required to implement RIPDG starting from an available IPDG
code are minimal and trivial. We envisage that the RIPDG method will be
effective in the numerical approximation of complex, multiscale problems
characterized by extreme local physical features necessitating highly non-
uniform Galerkin approximation spaces. In particular, as we shall see in the
numerical experiments below, RIPDG is particularly pertinent in scenarios
where _no_ conforming subspace of sufficient approximation capabilities is
available. This is often the case when employing IPDG methods on meshes
consisting of general polygonal/polyhedral (i.e., polytopic) elements,
typically with many faces per element. These meshes can arise, among other
processes, via mesh agglomeration procedures [11, 10, 9, 5]. Using elements
resulting from agglomeration of finer unstructured meshes, in an effort to
reduce computational complexity leads naturally to Galerkin spaces with
inherent significant non-conformity. In such cases, the effect of
discontinuity-penalization becomes evident and, in some extreme cases, appears
to be dominating the convergence behaviour (cf., Example 1 with the “zigzag”
mesh below).

At the other end of the spectrum, when using uniform meshes, constant
polynomial degree to solve problems with constant diffusion coefficients RIPDG
is identical to the classical IPDG. Moreover, within the usual setting of
locally quasi-uniform meshes, even with locally bounded variable polynomial
degree and small/medium contrast in the diffusion coefficient, RIPDG and IPDG
offer essentially identical numerical solutions.

Indeed, numerical experiments indicate the favourable performance of the new
RIPDG method over the classical version in terms of conditioning and error. In
particular, despite extensive testing, we have _not_ been able to identify an
example whereby the RIPDG method is inferior in terms of conditioning and/or
error in various norms when compared to the classical IPDG method. In most
cases the difference between IPDG and RIPDG is not significant, due to the
fact that a sufficiently rich conforming subspace is present in the
approximation. However, in certain extreme scenarios RIPDG offers
significantly better approximation and/or conditioning. Although we expect
that cases with IPDG outperforming RIPDG exist, we have not yet been able to
identify such example despite extensive testing.

The remainder of this work is structured as follows. In Section 2 we review
the classical IPDG method and the weighted variant from [15, 6, 16], along
with the basic argument for proving coercivity and continuity of the
respective bilinear form. In Section 3, we present the new RIPDG method and
discuss its construction and properties, as well as the motivating differences
in the choice of the penalty parameter. In Sections 4 and 6 we provide
extensions of RIPDG to polytopic meshes and to degenerate elliptic problems,
respectively. In Section 7, we present the basic Strang’s-type Lemma satisfied
by IPDG and RIPDG, showcasing that in the latter case the resulting a priori
error bound is _independent_ of any local variation in meshsize, polynomial
degree and diffusion contrast. We conclude with a series of numerical
experiments showcasing the benefits in using weighted averages in scenarios
with high such local variations.

##  2\. Interior penalty discontinuous Galerkin methods

We consider meshes $\mathcal{T}$ consisting of mutually disjoint open
simplicial and/or box-type elements $K\in\mathcal{T}$ so that
$\cup_{K\in\mathcal{T}}\bar{K}=\bar{\Omega}$. Let also
$h_{K}:=\operatorname{diam}(K)$ the diameter of $K\in\mathcal{T}$, and define
the mesh-function $\mathbf{h}:\cup_{K\in\mathcal{T}}K\to\mathbb{R}_{+}$ by
$\mathbf{h}|_{K}=h_{K}$, $K\in\mathcal{T}$. Let also $m_{K}$ denote the number
of $(d-1)$-dimensional faces of the element $K\in\mathcal{T}$, i.e.,
$m_{K}=d+1$ for a $d$-dimensional simplex and $m_{K}=2d$ for a $d$-dimensional
box-type element, and define
$\mathbf{m}:\cup_{K\in\mathcal{T}}K\to\mathbb{R}_{+}$ by
$\mathbf{m}|_{K}=m_{K}$, $K\in\mathcal{T}$. Further, we let
$\Gamma:=\cup_{K\in\mathcal{T}}\partial K$ denote the mesh skeleton and set
$\Gamma_{\text{\rm int}}:=\Gamma\backslash\partial\Omega$. The mesh skeleton
$\Gamma$ is decomposed into $(d-1)$–dimensional simplicial or quadrilateral
_faces_ $F$, shared by at most two elements. These are distinct from elemental
_interfaces_ $I$, herein defined as the simply-connected components of the
intersection between the boundaries of two neighbouring elements. As such,
hanging nodes/edges are naturally permitted: an interface between two elements
may consist of more than one faces residing on the same hyperplane in the
presence of hanging nodes. Moreover, as we will see below, we may wish to view
each quadrilateral (inter)face of a three-dimensional box-type element mesh as
the union of two simplicial faces. Let also
$\nabla_{\kern-2.27626pt\mathcal{T}}$ denote the _broken gradient_ defined as
$\nabla_{\kern-2.27626pt\mathcal{T}}v|_{K}:=\nabla(v|_{K})$,
$K\in\mathcal{T}$.

The _discontinuous Galerkin space_ $S_{\mathcal{T}}^{{\bf p}}$ with respect to
$\mathcal{T}$ is defined as

$S_{\mathcal{T}}^{{\bf p}}:=\\{v\in
L_{2}(\Omega):v|_{K}\in\mathcal{P}_{p_{K}}(K),\,K\in\mathcal{T}\\},$

with $\mathcal{P}_{p}(K)$ denoting the space of $d$-variate polynomials of
total degree up to $p$ on $K$, and
$\mathbf{p}:\cup_{K\in\mathcal{T}}K\to\mathbb{R}_{+}$ with
$\mathbf{p}|_{K}=p_{K}$, $K\in\mathcal{T}$. The local elemental polynomial
spaces employed within $S_{\mathcal{T}}^{{\bf p}}$ are defined in the
_physical coordinate system_ , i.e., without mapping from a given reference or
canonical frame (cf. [11, 10, 9]). All developments discussed below are also
valid for the more ‘classical’ case of the discontinuous Galerkin spaces
constructed through the usual mapping $\mathbf{F}_{K}:\hat{K}\to K$ from a
suitable reference element $\hat{K}$, viz.,

$S_{\mathcal{T},map}^{{\bf p}}:=\\{v\in
L_{2}(\Omega):v|_{K}\circ\mathbf{F}_{K}\in\mathcal{P}_{p_{K}}(\hat{K}),\,K\in\mathcal{T}\\}.$

Let $K_{+}$ and $K_{-}$ be two adjacent elements of $\mathcal{T}$ sharing a
face $F\subset\partial K_{+}\cap\partial K_{-}\subset\Gamma_{\text{\rm int}}$.
For $v$ and $\mathbf{q}$ element-wise continuous scalar- and vector-valued
functions, respectively, we define the _weighted average_ across $F$ by

$\\{\kern-3.1298pt\\{v\\}\kern-3.1298pt\\}_{w_{F}}|_{F}:=w_{+}v_{+}|_{F}+w_{-}v_{-}|_{F},\quad\\{\kern-3.1298pt\\{\mathbf{q}\\}\kern-3.1298pt\\}_{w_{F}}|_{F}:=w_{+}\mathbf{q}_{+}|_{F}+w_{-}\mathbf{q}_{-}|_{F},$

for $w_{F}:=(w_{+},w_{-})\in\mathbb{R}_{+}^{2}$ with $w_{+}+w_{-}=1$,
respectively, with $v_{\pm}|_{F}$ denoting the trace of $v$ from the element
$K_{\pm}$ on $F$, and correspondingly for $\mathbf{q}$. Also, we define the
_jump_ across $F$ by

$[\\![v]\\!]|_{F}:=v_{+}\mathbf{n}_{+}+v_{-}\mathbf{n}_{-},\quad[\\![\mathbf{q}]\\!]|_{F}:=\mathbf{q}_{+}\cdot\mathbf{n}_{+}+\mathbf{q}_{-}\cdot\mathbf{n}_{-},$

with $\mathbf{n}_{\pm}$ denoting the unit outward normal of $K_{\pm}$ on $F$.
We collect all weight pairs in a _weight function_
$\mathbf{w}:\Gamma_{\text{\rm int}}\to\mathbb{R}_{+}^{2}$, with
$\mathbf{w}|_{F}=w_{F}$, for each face $F\subset\Gamma_{\text{\rm int}}$. On a
boundary face $F\subset\partial\Omega\cap\partial K$, we set
$\\{\kern-3.1298pt\\{v\\}\kern-3.1298pt\\}:=v,$
$\\{\kern-3.1298pt\\{\mathbf{q}\\}\kern-3.1298pt\\}:=\mathbf{q},$
$[\\![v]\\!]:=v\mathbf{n},$ and
$[\\![\mathbf{q}]\\!]:=\mathbf{q}\cdot\mathbf{n}$, respectively, with
$\mathbf{n}$ the unit outward normal to $\partial\Omega$. In view of the
latter, we extend $\mathbf{w}$ to $\Gamma$ by setting
$\mathbf{w}|_{\partial\Omega}=(1,0)$ with $K_{-}=\emptyset$ there.

The (weighted) interior penalty discontinuous Galerkin method reads: find
$u_{h}\in S_{\mathcal{T}}^{{\bf p}}$, (or in $S_{\mathcal{T},map}^{{\bf p}}$,)
such that

(2.1) $B(u_{h},v_{h})=\ell(v_{h}),\qquad\text{for all }v_{h}\in
S_{\mathcal{T}}^{{\bf p}},$

with

$\displaystyle B(u_{h},v_{h}):=$ $\displaystyle\
\int_{\Omega}a\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}\cdot\nabla_{\kern-2.27626pt\mathcal{T}}v_{h}\,\mathrm{d}x+\int_{\Gamma}\sigma[\\![u_{h}]\\!]\cdot[\\![v_{h}]\\!]\,\mathrm{d}s$
$\displaystyle-\int_{\Gamma}\big{(}\\{\kern-3.1298pt\\{a\nabla
u_{h}\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]+\theta\\{\kern-3.1298pt\\{a\nabla
v_{h}\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![u_{h}]\\!]\big{)}\,\mathrm{d}s,$

and

$\displaystyle\ell(v_{h}):=$ $\displaystyle\
\int_{\Omega}fv_{h}\,\mathrm{d}x+\int_{\partial\Omega}g(\sigma v_{h}-\theta
a\nabla v_{h}\cdot\mathbf{n})\,\mathrm{d}s,$

with $\theta\in[-1,1]$ and $\sigma:\Gamma\to\mathbb{R}$ the, so-called,
_discontinuity-penalization_ or _penalty_ function, whose precise definition,
along with that of the weights $\mathbf{w}$ is a central concern in this work.
Note that, selecting $w_{F}=(1/2,1/2)$, for all $F\subset\Gamma_{\text{\rm
int}}$, we retrieve the classical IPDG method of Wheeler [28] and of Arnold
[2], while the case $w_{F}=(1,0)$ or $w_{F}=(0,1)$ for all
$F\subset\Gamma_{\text{\rm int}}$ gives the original method of Baker [4].
Moreover, in [15, 6], the weighted average choice
$w_{F}=(\delta_{+},\delta_{-})$ with

$\delta_{\pm}:=\frac{\alpha_{\mp}}{\alpha_{+}+\alpha_{-}},$

with $\alpha_{\pm}:=a|_{K_{\pm}}$, was proposed for the case of element-wise
constant scalar diffusion coefficients, while in [16] the selection
$\alpha_{\pm}:=\mathbf{n}_{F}^{T}a|_{K_{\pm}}\mathbf{n}_{F}$ for element-wise
constant diffusion tensors was provided, by combining the ideas from [15, 6]
and [18]. More recently, weighted averages involving also local meshsizes have
been proposed for treating interface problems [1, 20] on non-matching grids on
either side of the interface. These choices result to robust _a priori_ error
analysis with respect to the relative sizes of diffusion locally and/or the
relative meshsize across interfaces. The main contribution of this work is the
proposition of a different recipe for $\mathbf{w}$, which will yield robust
dependence not only with respect to the local variation of the diffusion, but
also with respect to the local _direction-wise_ variations in meshsize and
polynomial degree. This will, in turn, yield a robust and rigorous choice for
the discontinuity-penalization parameter $\sigma$.

The choice $\theta=1$ in (2.1) yields the _symmetric_ version, while the
choice $\theta=-1$ the non-symmetric version of the weighted IPDG method.
Although, we shall focus on the symmetric version in this work, as the most
popular choice, the developments presented below are also valid for the whole
range $\theta\in[-1,1]$, with trivial modifications of the arguments and of
the formulas.

To highlight the importance in the careful selection of $\sigma$, we consider
a $d$-dimensional simplicial triangulation $\mathcal{T}$, we set $\theta=1$,
$w_{F}=(1/2,1/2)$ for all faces $F\subset\Gamma_{\text{\rm int}}$, and we
study the coercivity of the bilinear form $B(\cdot,\cdot)$ on
$S_{\mathcal{T}}^{{\bf p}}$. To that end, for $v_{h}\in S_{\mathcal{T}}^{{\bf
p}}$, have

(2.2) $\displaystyle B(v_{h},v_{h})=$ $\displaystyle\
\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}v_{h}}\|_{L_{2}(\Omega)}^{2}+\|{\sqrt{\sigma}[\\![v_{h}]\\!]}\|_{L_{2}(\Gamma)}^{2}-2\int_{\Gamma}\\{\kern-3.1298pt\\{a\nabla
v_{h}\\}\kern-3.1298pt\\}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s,$

with
$\\{\kern-3.1298pt\\{\cdot\\}\kern-3.1298pt\\}:=\\{\kern-3.1298pt\\{\cdot\\}\kern-3.1298pt\\}_{(1/2,1/2)}$
the usual arithmetic average used in standard IPDG methods. To further
estimate the last term on the right-hand side of (2.2), we employ a _trace
inverse estimate_ , such as the one below.

######  Lemma 2.1.

Let $K\in\mathcal{T}$ a simplicial $d$-dimensional simplex and one of its
faces $F\subset\partial K$. For any $v\in\mathcal{P}_{p}(K)$, we have

(2.3)
$\|{v}\|_{L_{2}(F)}^{2}\leq\frac{(p+1)(p+d)}{d}\frac{|F|}{|K|}\|{v}\|_{L_{2}(K)}^{2}.$

###### Proof.

This follows immediately by combining a careful scaling argument with the main
result from [27]. ∎

We refer to [9] for a generalization of the above estimate to general curved
polygonal/polyhedral elements with arbitrary number of faces. Note that for
shape-regular elements, (2.3) yields the familiar version
$\|{v}\|_{L_{2}(\partial K)}^{2}\leq c_{\rm
inv}p^{2}/h_{K}\|{v}\|_{L_{2}(K)}^{2},$ with the constant $c_{\rm inv}>0$
depending on the aspect ratio of the element $K$ and the dimension $d$. If,
instead, we have $v_{h}\in S_{\mathcal{T},map}^{{\bf p}}$, then $c_{\rm inv}$
will also depend on the elemental maps $F_{K}$ if the latter are non-affine.

Returning to (2.2), elementary calculations, along with (2.3), for faces
$F\subset\partial K_{+}\cap\partial K_{-}$ give

$\displaystyle\Big{|}2\int_{\Gamma}\\{\kern-3.1298pt\\{a\nabla
u_{h}\\}\kern-3.1298pt\\}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s\Big{|}$
$\displaystyle\leq$
$\displaystyle{\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\int_{F}|a\nabla
u_{h}|_{K_{*}}[\\![v_{h}]\\!]|\,\mathrm{d}s}$ $\displaystyle\leq$
$\displaystyle\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!C_{\rm
inv,*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|{\nabla
u_{h}}\|_{L_{2}(K_{*})}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}$ $\displaystyle\leq$
$\displaystyle\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!C_{{\rm
inv},*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|a^{-\frac{1}{2}}\|_{L_{\infty}(K_{*})}\|{\sqrt{a}\nabla
u_{h}}\|_{L_{2}(K_{*})}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}$ $\displaystyle\leq$
$\displaystyle\
\frac{1}{2}\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}}\|_{L_{2}(\Omega)}^{2}+\frac{1}{2}\sum_{F\subset\Gamma}\tau_{F}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}^{2},$

with $C_{{\rm
inv},*}:=\sqrt{p_{K_{*}}(p_{K_{*}}\\!\\!\\!+d\\!-\\!1)|F|/(d|K_{*}|)}$,

(2.4) $\tau_{F}:=2\max_{*\in\\{+,-\\}}m_{K_{*}}C_{{\rm
inv},*}^{2}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}^{2}\|a^{-1}\|_{L_{\infty}(K_{*})},$

and $m_{K_{*}}$ the number of faces of the element $K_{*}$, defined above;
note that $\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}\in S_{\mathcal{T}}^{{\bf
p}-{\bf 1}}$, i.e., the discontinuous Galerkin space defined by the same mesh
having polynomial basis of one degree less than $S_{\mathcal{T}}^{{\bf p}}$.

These developments, in conjunction with (2.2), imply

$\displaystyle B(v_{h},v_{h})\geq$ $\displaystyle\
\frac{1}{2}\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}}\|_{L_{2}(\Omega)}^{2}+\sum_{F\subset\Gamma}\int_{F}\big{(}\sigma-\frac{\tau_{F}}{2}\big{)}[\\![u_{h}]\\!]^{2}\,\mathrm{d}s,$

upon setting $v_{h}=u_{h}$. Therefore, the selection $\sigma|_{F}=\tau_{F}$,
for all faces $F\subset\Gamma$ gives the coercivity estimate

(2.5)
$B(v_{h},v_{h})\geq\frac{1}{2}|\kern-0.7113pt\|{v_{h}}|\kern-0.7113pt\|^{2},\
\text{ with }\
|\kern-0.7113pt\|{v_{h}}|\kern-0.7113pt\|:=\sqrt{\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}}\|_{L_{2}(\Omega)}^{2}+\|{\sqrt{\sigma}[\\![v_{h}]\\!]}\|_{L_{2}(\Gamma)}^{2}}\,.$

The above choice of the discontinuity-penalization parameter $\sigma$ is, to
the best of our knowledge, the sharpest general estimate for simplicial
meshes. An analogous choice of $\sigma$ was proposed in [11, 9] for
significantly more complex element shapes. We note that hanging nodes are
permissible by viewing a simplex as a polytopic element with many co-planar
faces and by modifying $m_{K}$ accordingly. The case of box-type elements will
be discussed in more detail in Section 5 below.

Let us define now an extension of the bilinear form $B:S_{\mathcal{T}}^{{\bf
p}}\times S_{\mathcal{T}}^{{\bf p}}\to\mathbb{R}$ for the special case at
hand, to accept arguments from the larger space
$\mathcal{V}:=H^{1}(\Omega)+S_{\mathcal{T}}^{{\bf p}}$, viz.,
$B:\mathcal{V}\times\mathcal{V}\to\mathbb{R}$ with

$\displaystyle B({z},v):=$ $\displaystyle\
\int_{\Omega}a\nabla_{\kern-2.27626pt\mathcal{T}}{z}\cdot\nabla_{\kern-2.27626pt\mathcal{T}}v\,\mathrm{d}x+\int_{\Gamma}\sigma[\\![{z}]\\!]\cdot[\\![v]\\!]\,\mathrm{d}s$
$\displaystyle-\int_{\Gamma}\big{(}\\{\kern-3.1298pt\\{a\Pi_{\bf{p-1}}\nabla{z}\\}\kern-3.1298pt\\}\cdot[\\![v]\\!]+\\{\kern-3.1298pt\\{a\Pi_{\bf{p-1}}\nabla
v\\}\kern-3.1298pt\\}\cdot[\\![{z}]\\!]\big{)}\,\mathrm{d}s,$

for all ${z},v\in\mathcal{V}$, with $\Pi_{\bf{p-1}}:L_{2}(\Omega)\to
S_{\mathcal{T}}^{{\bf p}-{\bf 1}}$ denoting the orthogonal $L_{2}$-projection
operator onto $S_{\mathcal{T}}^{{\bf p}-{\bf 1}}$. Completely analogous
arguments to the ones presented for the proof of (2.5), along with the
stability of $\Pi_{\bf{p-1}}$ in the $L_{2}$-norm, result in the continuity
bound

(2.6)
$B({z},v)\leq\frac{3}{2}|\kern-0.7113pt\|{{z}}|\kern-0.7113pt\||\kern-0.7113pt\|{v}|\kern-0.7113pt\|,\qquad{z},v\in\mathcal{V}.$

######  Remark 2.2.

If a different Galerkin space than $S_{\mathcal{T}}^{{\bf p}}$ is used, one
has to modify the coercivity and continuity analysis and, correspondingly, the
definition of $\sigma$ accordingly. In particular, it may not be the case
anymore that $\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}\in
S_{\mathcal{T}}^{{\bf p}-{\bf 1}}$, but rather that
$\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}$ resides on the Galerkin space
itself.

##  3\. A robust interior penalty discontinuous Galerkin method

We now propose a _robust_ interior penalty discontinuous Galerkin method by
taking advantage of appropriately selected weights in (2.1). In particular, we
will select diffusion tensor, mesh and local polynomial degree dependent
weights. To that end, upon considering (2.1), we have

(3.1) $\displaystyle B(v_{h},v_{h})=$ $\displaystyle\
\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}v_{h}}\|_{L_{2}(\Omega)}^{2}+\|{\sqrt{\sigma}[\\![v_{h}]\\!]}\|_{L_{2}(\Gamma)}^{2}-2\int_{\Gamma}\\{\kern-3.1298pt\\{a\nabla
v_{h}\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s.$

Setting $C_{{\rm
inv},*}:=\sqrt{p_{K_{*}}(p_{K_{*}}\\!\\!\\!+d\\!-\\!1)|F|/(d|K_{*}|)}$, using
(2.3), and working as before, we get, for $F\subset\partial K_{+}\cap\partial
K_{-}$,

$\displaystyle\Big{|}2\int_{\Gamma}\\{\kern-3.1298pt\\{a\nabla
u_{h}\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s\Big{|}$
$\displaystyle\leq$ $\displaystyle\
2\\!\sum_{F\subset\Gamma}\int_{F}\Big{(}\sum_{*\in\\{+,-\\}}w_{*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}|\nabla
u_{h}|_{K_{*}}|\Big{)}|[\\![v_{h}]\\!]|\,\mathrm{d}s$ $\displaystyle\leq$
$\displaystyle\ 2\\!\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!C_{{\rm
inv},*}w_{*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|a^{-\frac{1}{2}}\|_{L_{\infty}(K_{*})}\|{\sqrt{a}\nabla
u_{h}}\|_{L_{2}(K_{*})}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}.$

Setting now

$\zeta_{*}:=\Big{(}2\sqrt{m_{K_{*}}}C_{{\rm
inv},*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|a^{-\frac{1}{2}}\|_{L_{\infty}(K_{*})}\Big{)}^{-1},$

for $*\in\\{+,-\\}$, and selecting

(3.2) $w_{*}=\frac{\zeta_{*}}{\zeta_{+}+\zeta_{-}},$

we deduce

$\displaystyle\Big{|}2\int_{\Gamma}\\{\kern-3.1298pt\\{a\nabla
u_{h}\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s\Big{|}\leq$
$\displaystyle\
\frac{1}{2}\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}}\|_{L_{2}(\Omega)}^{2}+\frac{1}{2}\sum_{F\subset\Gamma}(\zeta_{+}+\zeta_{-})^{-2}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}^{2}.$

A reasonable selection for the discontinuity-penalization parameter is,
therefore,

(3.3) $\sigma|_{F}=(\zeta_{+}+\zeta_{-})^{-2},$

with $\zeta_{-}=0$ in the case of a boundary face.

######  Definition 3.1.

The _robust interior penalty discontinuous Galerkin method (RIPDG)_ is the
Galerkin procedure defined by (2.1) with $\mathbf{w}=(w_{+},w_{-})$ for
$w_{*}$ given by (3.2), $*\in\\{+,-\\}$ and the discontinuity-penalization
parameter selected as in (3.3).

######  Remark 3.2.

The elementary identity
$(\zeta_{+}+\zeta_{-})^{-2}\leq(\min_{*\in\\{+,-\\}}2\zeta_{*})^{-2}$, results
into the bound

(3.4) $\sigma|_{F}\leq\min_{*\in\\{+,-\\}}m_{K_{*}}C^{2}_{{\rm
inv},*}\|a\mathbf{n}|_{K_{*}}\|^{2}_{L_{\infty}(F)}\|a^{-\frac{1}{2}}\|^{2}_{L_{\infty}(K_{*})}.$

The estimate (3.4) shows that the RIPDG discontinuity-penalization parameter
grows proportionally $\min_{*\in\\{+,-\\}}C_{{\rm inv},*}$ and is effectively
independent of the local size of the diffusion coefficient. In contrast, for
IPDG, the discontinuity-penalization parameter grows proportionally with
$\max_{*\in\\{+,-\\}}C_{{\rm inv},*}$ and is not robust with respect to the
contrast; cf., (2.4).

######  Remark 3.3.

The above argument still applies for different choices of element-wise
polynomial Galerkin spaces by replacing suitably $C_{{\rm inv},*}$ by an
available upper bound of the respective inverse estimate constant for the
space at hand.

######  Remark 3.4 (On inverse estimate constants).

For different choices of element-wise Galerkin spaces, we may end up with
unspecified/hard to estimate universal constants in the inverse estimate. We
note carefully that _any such constants cancel out in the definition of the
weights $\mathbf{w}$ and, thus, the weights are independent from such
constants_, ensuring the practical selection of weights in various
approximation space scenarios. Of course, as is the case for the standard IPDG
method, the RIPDG penalty parameter will depend proportionally on such
universal constants.

The proof of (2.6) is also then immediate upon considering the inconsistent
formulation

$\displaystyle B({z},v):=$ $\displaystyle\
\int_{\Omega}a\nabla_{\kern-2.27626pt\mathcal{T}}{z}\cdot\nabla_{\kern-2.27626pt\mathcal{T}}v\,\mathrm{d}x+\int_{\Gamma}\sigma[\\![{z}]\\!]\cdot[\\![v]\\!]\,\mathrm{d}s$
$\displaystyle-\int_{\Gamma}\big{(}{\\{\kern-3.1298pt\\{a\Pi_{\bf{p-1}}\nabla{z}\\}\kern-3.1298pt\\}_{\mathbf{w}}}\cdot[\\![v]\\!]+{\\{\kern-3.1298pt\\{a\Pi_{\bf{p-1}}\nabla
v\\}\kern-3.1298pt\\}_{\mathbf{w}}}\cdot[\\![{z}]\\!]\big{)}\,\mathrm{d}s,$

for all ${z},v\in\mathcal{V}$, resulting in the continuity bound

(3.5)
$B(w,v)\leq\frac{3}{2}|\kern-0.7113pt\|{{z}}|\kern-0.7113pt\||\kern-0.7113pt\|{v}|\kern-0.7113pt\|,\qquad{z},v\in\mathcal{V}.$

We now compare the effects on the size of $\sigma$ of the classical choice
$\mathbf{w}=(1/2,1/2)$ and the of the one proposed herein for the averaging
operator. Any significantly different behaviour is bound to occur in extreme
mesh and/or polynomial degree local variation scenarios. To that end, we
consider the problem (1.1) for $d=2$ with $a=I_{2\times 2}$, the identity
matrix. The latter’s solution is approximated via (2.1) for
$\mathbf{w}=(1/2,1/2)$ and also via (2.1) with the choice (3.2), over a mesh
consisting of two quadrilateral elements

$K_{1}=(0,1-\delta)\times(0,1),\qquad K_{2}=(1-\delta,1)\times(0,1),$

for $\delta<1/2$ and with uniform local polynomial degree $p$. The
discontinuity-penalization parameter will differ for the two variants only on
the interior face $F=\\{1-\delta\\}\times(0,1)$. For the classical IPDG method
((2.1) with $\mathbf{w}=(1/2,1/2)$), we calculate:

$\sigma_{F}\equiv\sigma_{F}^{\rm IP}=4p(p+1)/\delta,$

whereas for the choice (3.2), we get

$\sigma_{F}\equiv\sigma_{F}^{\rm
RIP}=\frac{8p(p+1)}{(\sqrt{1-\delta}+\sqrt{\delta})^{2}}.$

So, as $\delta\to 0$, we have $\sigma_{F}^{\rm IP}\to\infty$, while
$\sigma_{F}^{\rm RIP}\to 8p(p+1).$ This extreme scenario highlights vastly
different penalization pattern between the two variants. In the new choice of
the weighted version of IPDG ((2.1) with (3.2)), we observe _robust_ behaviour
with respect to the local mesh variation.

Completely analogously, we now set $\delta=1/2$, so that $K_{1}$ and $K_{2}$
are identical in terms of shape, and we set $p_{K_{1}}=1$ and $p_{K_{2}}=p>1$.
We then compute

$\sigma_{F}^{\rm IP}=8p(p+1)\to\infty,$

as $p\to\infty$, whereas

$\sigma_{F}^{\rm RIP}=16\big{(}(p(p+1))^{-1/2}+2^{-1/2}\big{)}^{-2}\to 32.$

Having the penalty parameter converging to infinity may become an issue in
terms of the conditioning of the resulting linear systems. Therefore, when a
conforming subspace of the DG space with the same approximation properties
exists, e.g., for the case of a simplicial mesh, one expects only conditioning
differences between the two variants. However, in the more practically
interesting case whereby there is _no_ conforming subspace of full order,
e.g., for the case of $S_{\mathcal{T}}^{{\bf p}}$ on a mesh $\mathcal{T}$
consisting of box-type (quadrliateral/hexahedral) elements, $\sigma_{F}^{\rm
IP}\to\infty$ may be compromising in terms of approximation also.

##  4\. Extension to general polytopic elements

The classical IPDG has been put forward as a method of choice for the
development of Galerkin methods on meshes consisting of extremely general
polygonal/polyhedral elements; see, e.g., [11, 10, 9, 5] and the references
therein. A key development for practical applications has been their ability
to admit provably stable approximations on elements $K$ with _arbitrary_
number of faces $F$, possibly with $|F|\ll|K|^{1/d}$. The weight selection in
the robust IPDG method presented above depends on the number of elemental
faces. Applying the method in the form (2.1) with the selection (3.2) is may
lead to spurious, excessive over-penalization as $m_{K_{*}}\to\infty$.

Fortunately, it is often possible to use _non-overlapping_ simplicial
subdivisions of the element $K$ to effectively arrive at a bounded $m_{K}$.
More specifically, a number of extremely mild geometric conditions are
presented in [9] (cf. also [10] for earlier ideas in this vein) so that a
polygonal/polyhedral element $K$, possibly containing curved faces, admits an
inverse estimate of the form

(4.1) $\|{v}\|_{L_{2}(F_{i})}^{2}\leq C^{2}_{{\rm
inv}}(p,F_{i},K)\|{v}\|_{L_{2}(K)}^{2},$

for any $v\in\mathcal{P}_{p}(K)$ and for $\\{F_{i}\\}_{i=1}^{m_{K}}$ a
mutually disjoint subdivision of $\partial K$, with $C_{{\rm
inv}}(p,F_{i},K)>0$ a concretely determined expression; see [9, Lemma 4.21]
for the exact formula. A key idea in [9] is that each $F_{i}$ is star-shaped
with respect to one interior point $x_{i}$ of $K$, $i=1,\dots,m_{K}$, _and_
may contain an arbitrary number of (possibly very small in size)
$(d-1)$-dimensional faces. Following the same line of argument as in Section
3, we can define a robust IPDG method on meshes consisting of general,
possibly curved, polygonal/polyhedral elements by selecting

$\zeta_{*}:=\Big{(}2\sqrt{m_{K_{*}}}C_{{\rm
inv}}(p,F_{i},K)\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|a^{-\frac{1}{2}}\|_{L_{\infty}(K_{*})}\Big{)}^{-1},$

and retaining the definitions of (3.2) and (3.3). We note carefully that in
the case of general polytopic meshes (4.1) may contain an unknown constant
which, nevertheless, cancels out in the definition of the weights as discussed
in Remark 3.4, yielding explicit choice for the weights $\mathbf{w}$.

##  5\. On box-type elements

In Section 3, the discussion focused on simplicial meshes, for which weights
and discontinuity-penalization parameters for RIPDG, that are free of unknown
constants are provided. The main technical tool has been the availability of a
trace inverse estimate with _explicitly known_ constants for the simplicial
reference element [27], also presented in Lemma 2.1 for an arbitrary simplex.

For box-type (quadrilateral, hexahedral, etc.) element meshes, classical
discontinuous Galerkin literature and popular available implementations in
software libraries propose the use of ‘$Q$-type’ element bases (i.e.,
polynomial spaces of degree $p$ in _each_ variable) in conjunction with
nonlinear element mappings. In other words, we are in the setting of
$S_{\mathcal{T},map}^{{\bf p}}$ defined above, in complete analogy with the
classical developments in conforming finite element methods. In this case, the
trace inverse estimate becomes

(5.1)
$\|{v}\|_{L_{2}(F)}^{2}\leq(p_{K}+1)^{2}\|J_{\mathbf{F}_{K}}|_{F\circ\mathbf{F}_{K}}\|_{L_{\infty}(F\circ\mathbf{F}_{K})}\|J_{\mathbf{F}_{K}}^{-1}\|_{L_{\infty}(K)}\|{v}\|_{L_{2}(K)}^{2},$

with $J_{\mathbf{F}_{K}}$ denoting the Jacobian of the mapping
$\mathbf{F}_{K}:(0,1)^{d}\to K$, for a box-type element $K\in\mathcal{T}$. The
proof follows by employing Fubini’s Theorem in conjunction with Lemma 2.1 for
$d=1$ to show the trace inverse estimate

$\|{v\circ\mathbf{F}_{K}}\|_{L_{2}(F\circ\mathbf{F}_{K})}^{2}\leq(p_{K}+1)^{2}\|{v\circ\mathbf{F}_{K}}\|_{L_{2}((0,1)^{d})}^{2},$

for any $(d-1)$-dimensional interface $F\circ\mathbf{F}_{K}$ of the reference
hypercube $(0,1)^{d}$, along with a standard scaling argument. This setting
was alluded to in Remark 3.4 above. Therefore, employing (5.1) to construct
the weights from RIPDG will involve the Jacobian of the element maps; this may
be somewhat inconvenient in certain practical scenarios.

An alternative point of view was proposed in [11] (see also [7, 10, 8, 13, 14,
9] for further results), whereby the finite element space is defined via
$S_{\mathcal{T}}^{{\bf p}}$ for box-type elements also, that is element-wise
polynomials of _total_ degree $p$ are employed even for box-type elements,
without the use of non-linear mappings. As it was shown numerically in [11, 8]
and proven in [14], this choice results into improved spectral convergence
under $p$-refinement in various settings, since less basis functions per
element are used. In this spirit, we may view a box-type element
$K\in\mathcal{T}$ as the union of $d$ non-overlapping simplicial elements and,
correspondingly each $(d-1)$-dimensional box-type interface can be viewed as a
union of more than one $(d-1)$-dimensional simplicial faces. In such a
setting, we can employ the developments of Section 4 to define the RIPDG
weights and discontinuity-penalization parameter, exactly like in the case of
a general polytopic element. The advantage of this approach is that nonlinear
mappings are not involved in the constants, at the expense of a “notional
subdivision” of the box-type element into non-overlapping, simplicial
subelements.. Note that this process is not increasing the number of degrees
of freedom: it is performed locally only for the definition of the weights and
of $\sigma$.

##  6\. Robust IPDG for degenerate elliptic problems

It is possible to have well-posed elliptic problems of the form (1.1) with
diffusion tensors $a$ for which $a\mathbf{n}=\mathbf{0}$, or even
$a=\mathbf{0}$, on lower-dimensional manifolds of $\Omega$ [24]. In such
cases, the choice of weights from (3.2) has to be revisited. To that end, we
consider the _inconsistent robust interior penalty discontinuous Galerkin
method_ reading: find $u_{h}\in S_{\mathcal{T}}^{{\bf p}}$, (or in
$S_{\mathcal{T},map}^{{\bf p}}$,) such that

(6.1) $\tilde{B}(u_{h},v_{h})=\tilde{\ell}(v_{h}),\qquad\text{for all
}v_{h}\in S_{\mathcal{T}}^{{\bf p}},$

with

$\displaystyle\tilde{B}(u_{h},v_{h}):=$ $\displaystyle\
\int_{\Omega}a\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}\cdot\nabla_{\kern-2.27626pt\mathcal{T}}v_{h}\,\mathrm{d}x+\int_{\Gamma}\sigma[\\![u_{h}]\\!]\cdot[\\![v_{h}]\\!]\,\mathrm{d}s$
$\displaystyle-\int_{\Gamma}\big{(}\\{\kern-3.1298pt\\{\sqrt{a}\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
u_{h})\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]+\theta\\{\kern-3.1298pt\\{\sqrt{a}\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
v_{h})\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![u_{h}]\\!]\big{)}\,\mathrm{d}s,$

and

$\displaystyle\ell(v_{h}):=$ $\displaystyle\
\int_{\Omega}fv_{h}\,\mathrm{d}x+\int_{\partial\Omega}g\big{(}\sigma
v_{h}-\theta\sqrt{a}\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
v_{h}\cdot\mathbf{n})\big{)}\,\mathrm{d}s,$

for $\theta\in[-1,1]$, with $\Pi:L_{2}(\Omega)\to S_{\mathcal{T}}^{{\bf p}}$
(or, $\Pi:L_{2}(\Omega)\to S_{\mathcal{T},map}^{{\bf p}}$,) denoting the
orthogonal $L_{2}$-projection operator onto the dG space. The method was first
presented in [18] for the case $\mathbf{w}=(1/2,1/2)$ and with
$\Pi_{\mathbf{p-1}}$ replaced by $\Pi_{\mathbf{p}}$. For simplicity of the
presentation, we shall again focus on the symmetric case $\theta=1$; the cases
$\theta\in[-1,1)$ follow via trivial modifications of the arguments given
below.

From (2.3), and working as before, we get, for $F\subset\partial
K_{+}\cap\partial K_{-}$,

$\displaystyle\Big{|}2\int_{\Gamma}\\{\kern-3.1298pt\\{\sqrt{a}\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
u_{h})\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s\Big{|}$
$\displaystyle\leq$ $\displaystyle\
2\sum_{F\subset\Gamma}\int_{F}\Big{(}\sum_{*\in\\{+,-\\}}w_{*}\|\sqrt{a}\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}|\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
u_{h})|_{K_{*}}|\Big{)}|[\\![v_{h}]\\!]|\,\mathrm{d}s$ $\displaystyle\leq$
$\displaystyle\ 2\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!C_{{\rm
inv},*}w_{*}\|\sqrt{a}\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|{\sqrt{a}\nabla
u_{h}}\|_{L_{2}(K_{*})}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)},$

using the stability of the orthogonal $L_{2}$-projection operator. Setting now

$\tilde{\zeta}_{*}:=\Big{(}2\sqrt{m_{K_{*}}}C_{{\rm
inv},*}\|\sqrt{a}\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\Big{)}^{-1},$

for $*\in\\{+,-\\}$, and selecting, again,

(6.2) $w_{*}=\frac{\tilde{\zeta}_{*}}{\tilde{\zeta}_{+}+\tilde{\zeta}_{-}},$

we deduce

$\displaystyle\Big{|}2\int_{\Gamma}\\{\kern-3.1298pt\\{\sqrt{a}\Pi_{\mathbf{p-1}}\sqrt{a}(\nabla
u_{h})\\}\kern-3.1298pt\\}_{\mathbf{w}}\cdot[\\![v_{h}]\\!]\,\mathrm{d}s\Big{|}\leq$
$\displaystyle\
\frac{1}{2}\|{\sqrt{a}\nabla_{\kern-2.27626pt\mathcal{T}}u_{h}}\|_{L_{2}(\Omega)}^{2}+\frac{1}{2}\sum_{F\subset\Gamma}\sigma|_{F}\|{[\\![v_{h}]\\!]}\|_{L_{2}(F)}^{2},$

with $\sigma|_{F}:=(\tilde{\zeta}_{+}+\tilde{\zeta}_{-})^{-2},$ for every face
$F\subset\Gamma$. If the limit $\sqrt{a}\mathbf{n}|_{K_{*}}\to\mathbf{0}$
(from within the element interior to the face) is valid a.e. on $F$ for
exactly one of $*\in\\{+,-\\}$, we have $w_{*}\to 1$, with the other weight in
the pair tending to zero; also, in this case we have $\sigma\to 0$. Further,
if we have $\sqrt{a}\mathbf{n}|_{K_{\pm}}\to\mathbf{0}$ from within the
element interior to the face for both traces a.e. on $F$, we adopt the
convention that both $w_{\pm}\to 0$ and $\sigma\to 0$ on $F$ (i.e., we no
longer enforce that $w_{+}+w_{-}=1$) .

##  7\. A priori error analysis

We now show that standard a priori error bounds hold for all the robust IPDG
methods described in the previous section, assuming only sufficient regularity
of the exact solution. In particular, _no_ ‘local bounded variation/local
quasi-uniformity’ mesh and polynomial degree assumptions are required, as is
standard in the error analysis of the classical IPDG method with
$\mathbf{w}=(1/2,1/2)$.

More specifically, let $u\in H^{3/2+\epsilon}(\Omega)$, $\epsilon>0$.
Elementary calculations imply

(7.1)
${B(u_{h},w_{h})=B(u,w_{h})+\int_{\Gamma}r(u)\cdot[\\![w_{h}]\\!]\,\mathrm{d}s,\qquad
w_{h}\in S_{\mathcal{T}}^{{\bf p}},}$

with $r(u):=\\{\kern-3.1298pt\\{a\nabla u-\Pi_{\mathbf{p-1}}(a\nabla
u)\\}\kern-3.1298pt\\}_{\mathbf{w}}$ for the robust IPDG from Section 3, or
$r(u):=\\{\kern-3.1298pt\\{a\nabla u-\sqrt{a}\Pi_{\mathbf{p-1}}(\sqrt{a}\nabla
u)\\}\kern-3.1298pt\\}_{\mathbf{w}}$ for the variant from Section 6.

Then, the coercivity, the continuity along with standard calculations imply

(7.2) $|\kern-0.7113pt\|{u-u_{h}}|\kern-0.7113pt\|\leq{4}\inf_{v_{h}\in
S_{\mathcal{T}}^{{\bf
p}}}|\kern-0.7113pt\|{u-v_{h}}|\kern-0.7113pt\|+2\sup_{0\neq w_{h}\in
S_{\mathcal{T}}^{{\bf
p}}}\frac{\int_{\Gamma}r(u)\cdot[\\![w_{h}]\\!]\,\mathrm{d}s}{|\kern-0.7113pt\|{w_{h}}|\kern-0.7113pt\|}.$

To estimate the inconsistency term, assuming further that $u\in
H^{2}(\Omega)$, we make use of the best approximation estimate from [12] (see
also [21, 17] for earlier results on box-type elements): there exists a
constant $C_{\rm ap}>0$, independent of $v$ and of $p$ such that

(7.3) ${\|{v-\Pi_{\mathbf{p}}v}\|_{{L_{2}(\hat{F})}}^{2}\leq C_{\rm
ap}(p+1)^{-1}\|{\nabla v}\|_{L_{2}(\hat{K})}^{2}},$

for any face $\hat{F}\subset\partial\hat{K}$ of a reference simplicial element
$\hat{K}$. Employing a standard scaling argument, (7.3) implies for an element
$K\in\mathcal{T}$ the following best approximation estimate:

(7.4) $\|{v-\Pi_{\mathbf{p}}v}\|_{{L_{2}(F)}}^{2}\leq C_{\rm
ap}\frac{|F|h_{K}^{2}}{|K|(p+1)}\|{\nabla v}\|_{L_{2}(K)}^{2}.$

Upon observing the identity $w_{*}\sigma^{-1/2}=\zeta_{*}$ (or
$w_{*}\sigma^{-1/2}=\tilde{\zeta}_{*}$ for the method of Section 6), we have,
respectively, for any $\mathbf{v}_{h}\in[S_{\mathcal{T}}^{{\bf p}-{\bf
1}}]^{d}$,

$\displaystyle\int_{\Gamma}r(u)\cdot[\\![w_{h}]\\!]\,\mathrm{d}s$
$\displaystyle\leq$
$\displaystyle\sum_{F\subset\Gamma}\\!\sum_{*\in\\{+,-\\}}\\!\\!\\!\zeta_{*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}\|{\\!\big{(}\nabla
u-\mathbf{v}_{h}-\Pi_{\mathbf{p-1}}(\nabla
u-\mathbf{v}_{h})\big{)}\\!|_{K_{*}}\\!}\|_{{L_{2}(F)}}\|{\sqrt{\sigma}[\\![w_{h}]\\!]}\|_{{L_{2}(F)}}$
$\displaystyle\leq$
$\displaystyle\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!\frac{\sqrt{C_{{\rm
ap}}|F|}h_{K_{*}}}{\sqrt{|K_{*}|p_{K_{*}}}}\zeta_{*}\|a\mathbf{n}|_{K_{*}}\|_{L_{\infty}(F)}|\nabla
u-\mathbf{v}_{h}|_{H^{1}(K_{*})}\|{\sqrt{\sigma}[\\![w_{h}]\\!]}\|_{{L_{2}(F)}}$
$\displaystyle\leq$
$\displaystyle\sum_{F\subset\Gamma}\sum_{*\in\\{+,-\\}}\\!\\!\\!\frac{\|\sqrt{a}\|_{L_{\infty}(K_{*})}\sqrt{dC_{{\rm
ap}}}}{{2}\sqrt{m_{K_{*}}}\sqrt{p_{K_{*}}\\!\\!\\!+d\\!-\\!1}}\frac{h_{K_{*}}}{p_{K_{*}}}|\nabla
u-\mathbf{v}_{h}|_{H^{1}(K_{*})}\|{\sqrt{\sigma}[\\![w_{h}]\\!]}\|_{{L_{2}(F)}}$
$\displaystyle\leq$ $\displaystyle\ \Big{(}\frac{dC_{{\rm
ap}}}{{2}}\sum_{K\in\mathcal{T}}\|\sqrt{a}\|_{L_{\infty}(K)}^{2}\frac{h^{2}_{K}}{p_{K}^{3}}|\nabla
u-\mathbf{v}_{h}|_{H^{1}(K)}^{2}\Big{)}^{\frac{1}{2}}\|{\sqrt{\sigma}[\\![w_{h}]\\!]}\|_{{L_{2}(\Gamma)}},$

using (7.4) in the first step and the definition of $\zeta_{*}$ in the
penultimate step. The last estimate concerns the robust IPDG method from
Section 3 and holds under the regularity assumption $u\in H^{2}(\Omega)$.

Correspondingly, for the method of Section 6, assuming that the $a$ is such
that $\sqrt{a}\nabla u\in[H^{1}(\Omega)]^{d}$, (_instead of_ $u\in
H^{2}(\Omega)$ as before,) we arrive at

$\displaystyle\int_{\Gamma}r(u)\cdot[\\![w_{h}]\\!]\,\mathrm{d}s\leq$
$\displaystyle\ \Big{(}\frac{dC_{{\rm
ap}}}{{2}}\sum_{K\in\mathcal{T}}\frac{h^{2}_{K}}{p_{K}^{3}}|\sqrt{a}\nabla
u-\mathbf{v}_{h}|_{H^{1}(K)}^{2}\Big{)}^{\frac{1}{2}}\|{\sqrt{\sigma}[\\![w_{h}]\\!]}\|_{{L_{2}(\Gamma)}}.$

Combining the above developments with (7.2), we deduce

(7.5) $\displaystyle|\kern-0.7113pt\|{u-u_{h}}|\kern-0.7113pt\|\leq$
$\displaystyle\ {4}\inf_{v_{h}\in S_{\mathcal{T}}^{{\bf
p}}}|\kern-0.7113pt\|{u-v_{h}}|\kern-0.7113pt\|$ $\displaystyle\
+\sqrt{{2}dC_{\rm ap}}\inf_{\mathbf{v}_{h}\in[S_{\mathcal{T}}^{{\bf p}-{\bf
1}}]^{d}}\Big{(}\sum_{K\in\mathcal{T}}\|\sqrt{a}\|_{L_{\infty}(K)}^{2}\frac{h^{2}_{K}}{p_{K}^{3}}|\nabla
u-\mathbf{v}_{h}|_{H^{1}(K)}^{2}\Big{)}^{\frac{1}{2}},$

for the robust IPDG of Section 3, or

$|\kern-0.7113pt\|{u-u_{h}}|\kern-0.7113pt\|\leq{4}\inf_{v_{h}\in
S_{\mathcal{T}}^{{\bf
p}}}|\kern-0.7113pt\|{u-v_{h}}|\kern-0.7113pt\|+\sqrt{{2}dC_{\rm
ap}}\inf_{\mathbf{v}_{h}\in[S_{\mathcal{T}}^{{\bf p}-{\bf
1}}]^{d}}\Big{(}\sum_{K\in\mathcal{T}}\frac{h^{2}_{K}}{p_{K}^{3}}|\sqrt{a}\nabla
u-\mathbf{v}_{h}|_{H^{1}(K)}^{2}\Big{)}^{\frac{1}{2}},$

for the method of Section 6, respectively. These basic estimates, together
with standard $hp$-version best approximation results, give rise to standard a
priori error bounds. Although not considered in detail here in the interest of
brevity, following essentially the same steps as above, it is also possible to
prove a priori error bounds for the robust IPDG method on polygonal/polyhedral
elements from Section 4, with the possible exception of considering simplicial
coverings of the elements $K\in\mathcal{T}$ along with extension operators; we
refer to [11, 10, 9] for details.

The crucial development in the a priori error bounds above, compared to all
respective bounds we are aware of in the literature, is that the best
approximation estimates are _completely localized element-wise_. This is in
contrast to standard IPDG error bounds requiring local quasi-
uniformity/bounded variation assumptions for the mesh size and local
polynomial degree. This property is expected to be of importance on
approximation space scenarios admitting extreme local variation.

######  Remark 7.1.

In the above best approximation results we have taken a somewhat “classical”
viewpoint of assuming sufficient regularity of the exact solution. More recent
developments e.g. the so-called medius analysis [19], or the use of recovery
operators [26], can provide quasi-optimality under minimal regularity
assumptions on the exact solution. These approaches, however, are currently
available only for the $h$-version IPDG under standard local quasi-uniformity
mesh and diffusion coefficient bounded local variation assumptions for which
the robust IPDG method presented in this work offers little or no advantage.
Nonetheless, the extension of such “optimal” error analysis frameworks to the
RIPDG setting is certainly of interest.

##  8\. Numerical Experiments

We now present a series of numerical experiments showcasing the comparative
performance between the classical IPDG and RIPDG. All numerical experiments
below use the locally _total degree $p$_ approximation space per element
$S_{\mathcal{T}}^{{\bf p}}$ on box-type or general polytopic element shapes,
noting carefully that no conforming subspace of $S_{\mathcal{T}}^{{\bf p}}$ of
the same degree is available. This choice of approximation spaces has been
advocated in [11, 10, 9] leading to optimal approximation yet reduced
complexity for box-type elements.

Given that the effect of high contrast diffusion coefficient has been
discussed extensively in [15, 6, 16], we confine our study to the effects of
extreme mesh and polynomial degree variation. Nonetheless, we expect that the
inclusion of high contrast diffusion coefficient in the spirit of [6, 16] in
conjunction with extreme local variation scenarios will only augment the
behaviours presented below. We stress that in _all_ comparisons presented
below both methods have exactly the same numerical degrees of freedom as they
are applied to the same approximation space; the two respective
implementations differ _only_ on the choice of the weights $\bf{w}$ and,
correspondingly, to the choice of the penalty parameter. Therefore, any
differences observed in performance are due to the different methods
implemented in the otherwise identical program, and not on any other
algorithmic factors (e.g., choice of quadrature or condition number estimation
routines, etc.).

###  8.1. Example 1: singularly perturbed reaction-diffusion problem

$l$$l$

Figure 1. Example 1. The $9$ element mesh.

We begin by testing the method on anisotropic rectangular elements in the
context of layer adapted $hp$-version IPDG/RIPDG method. Consider the problem

$-\epsilon\Delta u+u=f,\qquad\text{on }\Omega:=(-1,1)^{2},$

with homogeneous Dirichlet boundary conditions and $f$ an analytic function
chosen so that the exact solution reads

(8.1)
$u(x,y):=\Big{(}1-\frac{\cosh(x/\sqrt{\epsilon})}{\cosh(1/\sqrt{\epsilon})}\Big{)}\Big{(}1-\frac{\cosh(y/\sqrt{\epsilon})}{\cosh(1/\sqrt{\epsilon})}\Big{)}.$

This example has been studied extensively in [22, 17]. The solution exhibits
boundary layers of thickness $\mathcal{O}(\sqrt{\epsilon})$.

To resolve the layers, we first use a layer-adapted anisotropic $9$-element
mesh with characteristic width $l:=\min\\{\lambda p\sqrt{\epsilon},0.5\\}$,
and some constant $\lambda>0$; see Figure 1 for an illustration. It was shown
in [22] that conforming Galerkin methods on such meshes offer a robust
exponential convergence under $p$ refinement. Note that the maximum meshsize
local variation constant $r$ for the $9$-element mesh is given by
$r:={2(1-l)}/{l}.$ For $\epsilon\ll 1$, we have $r\approx 2/(\lambda
p\sqrt{\epsilon}).$

We compare the symmetric versions of the classical IPDG against the RIPDG
method on the same meshes for polynomial degrees $p=1,\dots,7$. In Figure 2,
it is observed that both methods converge exponentially in both dG-norm and
broken $H^{1}$-seminorm
$|\cdot|_{H^{1}(\Omega,\mathcal{T})}:=\|\nabla_{\kern-2.27626pt\mathcal{T}}\cdot\|_{L_{2}(\Omega)}$
against the square root of the total numerical degrees of freedom ($DoFs$)
under $p$ refinement for $\epsilon=10^{-5}$ and $\lambda=0.9$. The RIPDG
appears to outperform IPDG in both measures.

Figure 2. Example 1. Convergence in dG-norm (left) and broken $H^{1}$-seminorm
(right) for $\epsilon=10^{-5}$ and $p=1,\dots,7$.

Figure 3. Example 1. Maximum of the penalty parameter (left) and the condition
number of the linear system (right) for $\epsilon=10^{-5}$ and $p=1,\dots,7$.

In Figure 3 (left), we compare the magnitudes of the global maximum of the
penalty parameter values $\max_{F\subset\Gamma}\sigma_{F}^{\rm IP}$ and
$\max_{F\subset\Gamma}\sigma_{F}^{\rm RIP}$, respectively, as well as the
condition numbers for each method based on the same choice of basis. As
expected, RIPDG requires a far smaller penalty parameter in theory for
stability compared to the standard IPDG. For $p=1$,
$\max_{F\subset\Gamma}\sigma_{F}^{\rm IP}$ is $120$ times larger than
$\max_{F\subset\Gamma}\sigma_{F}^{\rm RIP}$. This difference in the definition
of the penalty parameter explains also the improvement in the condition number
for RIPDG compared to IPDG by a factor  about 1.3 for $p=1,\cdots,7$; see
Figure 3 (right).

In the above experiment, we observe only a modest improvement on error and the
conditioning by using RIPDG as opposed to IPDG. We highlight that the
$hp$-version spaces used for convergence admit a conforming subspace of the
same approximation capabilities.

Figure 4. Example 1. “Zigzag” meshes for $p=2$ (top left), $p=3$ (top right),
$p=5$ (bottom left), and $p=7$ (bottom right).

We now employ a “zigzag” polygonal version of the $9$ element mesh shown in
Figure 4 for $p=2,3,5,7$. In particular, we replace the interior faces from
the Figure 1 by “zigzag” curves of amplitude $l/6$ with respect to the length
scale $l=\min\\{\lambda p\sqrt{\epsilon},0.5\\}$; the latter still
characterises the distance of each interface to the non-intersecting portion
of the boundary. On this $9$ element polytopic mesh, we solve the same problem
using the RIPDG and IPDG methods, on physical (i.e., unmapped) element-wise
polynomial spaces of uniform degree $p$ for ${p=2},\dots,8$ and
$\epsilon=10^{-3}$, with the penalty choice informed by the inverse estimate
(4.1). We observe that in this mesh _no_ conforming subspace of sufficient
approximation capabilities is available.

In Figure 5, we compare the magnitudes of the global maximum of the penalty
parameter values and the condition numbers for each method based on the same
choice of basis. In Figure 6, we compare the respect errors for $p=1,\dots,8$.
In Figure 6, it is observed that both methods converge exponentially in dG-
norm, and broken $H^{1}$-seminorm against the square root of the total
numerical degrees of freedom ($DoFs$) under $p$ refinement for
$\epsilon=10^{-3}$ and $\lambda=0.9$. The RIPDG appears to significantly
outperform IPDG for higher $p$; the IPDG appears to stagnate for $p=7,8$. We
believe that the reason behind this behaviour is the IPDG’s penalization
magnitude in conjunction with the essential non-conformity of the
approximation space. On the other hand, RIPDG’s $p$-convergence appears not to
be affected.

Figure 5. Example 1. “Zigzag” mesh. Maximum of the penalty parameter (left)
and the condition number of the linear system (right) for $\epsilon=10^{-3}$
and $p=1,\dots,8$.

Figure 6. Example 1. “Zigzag” mesh. Convergence in dG-norm (left), broken
$H^{1}$-seminorm (right), for $\epsilon=10^{-3}$ and $p=1,\dots,8$.

###  8.2. Example 2: Poisson problem with Gaussian solutions.

We now test the behaviour of each method against rapidly changing local
polynomial degree $p$. To that end, in $\Omega:=(-1,1)^{2}$, we consider the
Poisson problem, with $a=I_{2\times 2}$ the $2\times 2$-identity matrix, with
Dirichlet boundary conditions, and we choose $f$ such that the exact solution
$u$ is given by

(8.2) $u(x,y):=\exp(-\alpha(x^{2}+y^{2})),$

for some $\alpha>0$, which is analytic; in Figure 7 (left) $u$ for
$\alpha=100$ is depicted.

Figure 7. Example 2. $u$ for $\alpha=100$ (left) and mesh and local polynomial
degree distribution for the results of Table 1 (right).

For the first numerical experiment, we employ a $9$-element uniform square
mesh polynomial degree $p=2$ in the eight boundary elements and polynomial
degree $p=30$ on the interior $\tilde{K}:=(-1/3,1/3)^{2}$; this approximation
space is presented in Figure 7 (right). In Table 1, we present the comparison
between IPDG and RIPDG in terms of the maximum penalty used, condition number
of the respective stiffness matrices and three different error measures. As we
can see, RIPDG is outperforming in all respects, with the most striking
improvement being in the stiffness matrix condition number.

method | $\max_{F\subset\Gamma}\sigma_{F}$ | condition no. | $\|u-u_{h}\|$ | $|u-u_{h}|_{H^{1}(\Omega,\mathcal{T})}$ | $|||u-u_{h}|||$
---|---|---|---|---|---
IPDG | 5.5800e+03 | 5.1229e+06 | 7.1923e-06 | 1.9711e-04 | 2.1681e-04
RIPDG | 61.69 | 5.1148e+05 | 6.5842e-06 | 1.9169e-04 | 2.1305e-04
Table 1. Example 2. $\alpha=100$. Comparison on the mesh and polynomial degree
distribution given in Figure 7 (right).

The significant difference in the maximum of the penalty parameter confirms
the respective discussion in Section 3. As a consequence, we can see that
IPDG’s condition number is about $10$ times larger than RIPDG’s for exactly
the same problem and approximation space.

Finally, we demonstrate the qualitative behaviour between the two methods
resulting from the significant differences in the penalization magnitude. To
that end, we consider the same problem with $\alpha=10$, resulting to a
significantly less sharp Gaussian bump at the centre of the domain. This
problem is approximated by the mesh given in Figure 8, with each elemental
polynomial degree also given.

Figure 8. Mesh used for the solutions in Figure 9 with $p=1$ on all ‘small’
elements and $p=5,8$ in the central ‘large’ element.

In Figure 9 we provide the solution profiles for RIPDG (left) and IPDG (right)
for $\alpha=10$, $p=1$ on all ‘small’ elements and $p=3$ (first line), $p=5$
(second line), and $p=8$ (last line) in the central ‘large’ element. We
observe significantly different behaviour in the ‘large’ central element
interfaces, due to the substantially smaller penalization required by the
RIPDG method. Interestingly, however, we also observe modestly different
behaviour further away from the ‘large’ element interfaces also, especially
for $p=3$.

Figure 9. Solution profiles for RIPDG (left) and IPDG (right) for $\alpha=10$,
$p=1$ on all ‘small’ elements and $p=3$ (first line), $p=5$ (second line), and
$p=8$ (last line) in the central ‘large’ element.

Of course, one may not use encounter or use approximations spaces as the ones
considered in this example. Nonetheless, such mesh and polynomial degree
distributions may arise locally in the context of $hp$-adaptive procedures. In
any case, this study is revealing in how extreme approximation space local
variation behaviour may result into possibly inferior performance of classical
approaches.

###  8.3. Example 3. Highly agglomerated meshes.

We now consider the Poisson problem with Dirichlet boundary conditions
admitting the smooth solution $u:=\sin(\pi x)\sin(\pi y)$ on highly
agglomerated meshes: we start from a fine simplicial mesh with $131,072$
triangles which is further agglomerated in a random fashion using standard
partitioning tools into $37$ complicated polygonal elements shown in Figure
10. The discontinuity penalization parameter is selected as described briefly
in Section 4; we refer to [9] for details. On this fixed mesh, we compute the
numerical solution for each method for uniform $p=1,\dots,5$; these are
presented in Figure 11. As expected, the errors in all norms decay
exponentially. Even though the polynomial degree is uniform and neighbouring
elements are of ‘similar’ sizes, we can still observe an improvement in the
errors for RIPDG against IPDG of about 25%. In Figure 12, we can see that the
maximum penalty and condition numbers for IPDG compared to RIPDG are
approximately $4$ and $2.5$ times larger, respectively.

Figure 10. $37$ agglomerated elements made by $131072$ triangular meshes

Figure 11. Exponential convergence in dG–norm (Left) and broken $H^{1}$–norm
(Right) with $p=1,\dots,5$

Figure 12. The maximum penalty parameters (Left) and the condition number of
the linear system(Right) $p=1,\dots,5$

## References

  * [1] C. Annavarapu, M. Hautefeuille, and J. E. Dolbow, A robust Nitsche’s formulation for interface problems, Comput. Methods Appl. Mech. Engrg., 225/228 (2012), pp. 44–54.
  * [2] D. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal., 19 (1982), pp. 742–760.
  * [3] I. Babuška, The finite element method with penalty, Math. Comp., 27 (1973), pp. 221–228.
  * [4] G. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comp., 31 (1977), pp. 45–59.
  * [5] F. Bassi, L. Botti, A. Colombo, D. Di Pietro, and P. Tesini, On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations, J. Comput. Phys., 231 (2012), pp. 45–65.
  * [6] E. Burman and P. Zunino, A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems, SIAM J. Numer. Anal., 44 (2006), pp. 1612–1638.
  * [7] A. Cangiani, Z. Dong, E. Georgoulis, and P. Houston, $hp$–Version discontinuous Galerkin methods for advection–diffusion–reaction problems on polytopic meshes, ESAIM: M2AN, 50 (2016), pp. 699–725.
  * [8] A. Cangiani, Z. Dong, and E. H. Georgoulis, $hp$-version space-time discontinuous Galerkin methods for parabolic problems on prismatic meshes, SIAM J. Sci. Comput., 39 (2017), pp. A1251–A1279.
  * [9] , $hp$-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements, Math. Comp., 91 (2022), pp. 1–35.
  * [10] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston, $hp$-version discontinuous Galerkin methods on polygonal and polyhedral meshes, SpringerBriefs in Mathematics, Springer, Cham, 2017.
  * [11] A. Cangiani, E. Georgoulis, and P. Houston, $hp$–Version discontinuous Galerkin methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci., 24 (2014), pp. 2009–2041.
  * [12] A. Chernov, Optimal convergence estimates for the trace of the polynomial $L^{2}$-projection operator on a simplex, Math. Comp., 81 (2012), pp. 765–787.
  * [13] Z. Dong, Discontinuous galerkin methods for the biharmonic problem on polygonal and polyhedral meshes, Int. J. Numer. Anal. Model., 16 (2019), pp. 825–846.
  * [14] Z. Dong, On the exponent of exponential convergence of $p$-version FEM spaces, Adv. Comput. Math., 45 (2019), pp. 757–785.
  * [15] M. Dryja, On discontinuous Galerkin methods for elliptic problems with discontinuous coefficients, vol. 3, 2003, pp. 76–85.  Dedicated to Raytcho Lazarov.
  * [16] A. Ern, A. F. Stephansen, and P. Zunino, A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity, IMA J. Numer. Anal., 29 (2009), pp. 235–256.
  * [17] E. Georgoulis, Discontinuous Galerkin methods on shape-regular and anisotropic meshes, D.Phil. Thesis, University of Oxford, (2003).
  * [18] E. Georgoulis and A. Lasis, A note on the design of $hp$-version interior penalty discontinuous Galerkin finite element methods for degenerate problems, IMA J. Numer. Anal., 26 (2006), pp. 381–390.
  * [19] T. Gudi, A new error analysis for discontinuous finite element methods for linear elliptic problems, Math. Comp., 79 (2010), pp. 2169–2189.
  * [20] T. Gustafsson, R. Stenberg, and J. Videman, On Nitsche’s method for elastic contact problems, SIAM J. Sci. Comput., 42 (2020), pp. B425–B446.
  * [21] P. Houston, C. Schwab, and E. Süli, Discontinuous $hp$-finite element methods for advection-diffusion-reaction problems, SIAM J. Numer. Anal., 39 (2002), pp. 2133–2163.
  * [22] J. Melenk and C. Schwab, An $hp$–finite element method for convection-diffusion problems in one dimension, IMA J. Numer. Anal., 19 (1999), pp. 425–453.
  * [23] J. Nitsche, Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind, in Abhandlungen aus dem mathematischen Seminar der Universität Hamburg, vol. 36, Springer, 1971, pp. 9–15.
  * [24] O. Oleinik and E. Radkevič, Second Order Equations with Nonnegative Characteristic Form, American Mathematical Society, 1973.
  * [25] B. Rivière, M. Wheeler, and V. Girault, A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal., 39 (2001), pp. 902–931.
  * [26] A. Veeser and P. Zanotti, Quasi-optimal nonconforming methods for symmetric elliptic problems. III—Discontinuous Galerkin and other interior penalty methods, SIAM J. Numer. Anal., 56 (2018), pp. 2871–2894.
  * [27] T. Warburton and J. S. Hesthaven, On the constants in $hp$-finite element trace inverse inequalities, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 2765–2773.
  * [28] M. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal., 15 (1978), pp. 152–161.