# Conditional Idempotent Generative Networks

Niccolò Ronchetti

(June 3, 2024)

###### Abstract

We propose Conditional Idempotent Generative Networks (CIGN), a novel approach
that expands upon Idempotent Generative Networks (IGN) to enable conditional
generation. While IGNs offer efficient single-pass generation, they lack the
ability to control the content of the generated data. CIGNs address this
limitation by incorporating conditioning mechanisms, allowing users to steer
the generation process towards specific types of data.

We establish the theoretical foundations for CIGNs, outlining their scope,
loss function design, and evaluation metrics. We then present two potential
architectures for implementing CIGNs: channel conditioning and filter
conditioning. Finally, we discuss experimental results on the MNIST dataset,
demonstrating the effectiveness of both approaches. Our findings pave the way
for further exploration of CIGNs on larger datasets and with more powerful
computing resources to determine the optimal implementation strategy.

##  1 Introduction

In the paper [10], the authors propose a novel approach to generation which
they call Idempotent Generative Network (IGN). One major advantage of this
approach compared to diffusion (or autoregressive) models is that generation
is done in a single forward pass, and is therefore several times faster.

However, the paper [10] does not touch into how to _condition_ that
generation: a well-trained idempotent generative network will generate
synthetic data resembling examples seen during training, but the authors do
not describe how to nudge the IGN towards generating certain examples (for
example, examples belonging to a certain class) over others.

Our paper starts answering that question: we describe how to expand the theory
of idempotent generative networks to include a conditioning mechanism, and we
call these networks _Conditional Idempotent Generative Networks_. We further
describe two possible architectures for implementing conditional idempotent
generative networks, and we discuss experimental results.

Our main contributions are then:

  1. 1.

We set the theoretical foundations for Conditional Idempotent Generative
Networks, describing their scope and intended use, the format of the loss
function and what metrics can be used to evaluate them.

  2. 2.

We describe two possible implementations of the conditioning mechanism of a
Conditional Idempotent Generative Network - we call these mechanisms _channel
conditioning_ and _filter conditioning_.

  3. 3.

We discuss results of our experiments comparing the two implementations of
CIGN on the MNIST dataset.

Our experiments show that both channel conditioning and filter conditioning
are effective approaches in implementing a Conditional Idempotent Generative
Networks. Further experiments on larger datasets and with more powerful
computing power are needed to identify which mechanism is more effective.

Our code is available at
https://github.com/niccronc/conditional-idempotent-generative-networks.

##  2 Idempotent generative networks

We start by summarizing the ideas leading to the concept of idempotent
generative networks, after [10].

The authors start from the observation that typically the true data
distribution lives in a low-dimensional submanifold $D$ of the feature space
$X$. _Generating_ a synthetic data point corresponds then to choosing a point
on this low-dimensional submanifold $D$. One possible way to choose a point on
this low-dimensional submanifold $D$ is to have a function $F$ whose image is
$D$, and then apply $F$ to a random input.

The authors’ idea is to construct $F$ to be an idempotent function111A
function $f:X\longrightarrow X$ is called idempotent if, when applied multiple
times, it produces the same output as applying it once: $(f(f(x))=f(x))$. on
the feature space, with image constrained to be _as close as possible_ to $D$.

A machine learning network (called _idempotent generative network_) is then
trained to learn the parameters of the idempotent function $F:X\longrightarrow
X$ with image $D$.

The difficulty of creating synthetic data points (in other words - the
difficulty of describing an accurate parametrization of the submanifold $D$
starting from a training dataset) is then subsumed into the difficulty of
identifying optimal parameter for the function $F$. A well-trained idempotent
generative network $F$ will then generate data points $F(x)$ when applied to
any randomly selected point $x$ in the feature space $X$. In particular,
generation happens in a single forward pass and is therefore considerably
faster than autoregressive or diffusion generative methods.

We summarize the benefits of using idempotent generative networks:

  * •

Single-step generation: unlike diffusion or autoregressive models, IGNs
generate outputs in one step.

  * •

Optional refinements: similar to diffusion models, IGNs allow for optional
sequential refinements (meaning: we can pass the output through the network a
second time, and since the trained $\hat{F}$ is not truly idempotent, we may
obtain a slightly higher quality output - this was discussed in [10]).

  * •

No feature/latent space conundrum: unlike diffusion models, IGNs work only on
the feature space, facilitating manipulations and interpolations of inputs and
outputs.

  * •

Output enhancement: an IGN could be used to project degraded data back onto
the original data distribution.

###  2.1 Loss function

While the initial version of the paper [10] does not provide insights on the
model architectures used for the experiments the authors describe, the paper
does discuss at length the loss function of an IGN from a theoretical
perspective.

Let again $X$ be the feature space and $D$ be the data manifold. Let
$F:X\longrightarrow X$ be the idempotent function we are trying to construct
and whose parameters are being learned by the IGN.

The loss function is then composed of the following three terms:

  * •

(reconstruction) Real data samples $d$ from the data distribution $D$ should
remain unchanged by the model: $F(d)=d$.

  * •

(idempotent) Any further application of the model beyond the first one should
act like the identity function: $F(F(z))=F(z)$ for every $z$ in $X$.

  * •

(tightness) The previous two conditions should be satisfied while ensuring $F$
has as small a range as possible.

We encourage the interested reader to check out additional details in the
original paper [10].

###  2.2 Network architecture

The theoretical framework of idempotent generative network is fundamentally
architecture-agnostic. And yet, one needs to choose a model architecture in
order to test the feasibility of the idea in practice.

In [10], the authors do not explicitly describe the architecture used in their
experiments. Inspired by [8], we will leverage a GAN-like architecture.

In the classic setup of Generative Adversarial Networks, a generator and a
discriminator compete in a zero-sum game: the generator creates synthetic data
points resembling as much as possible a real dataset, while the discriminator
attempts to distinguish the real data points from the synthetic ones. In
particular, in the GAN framework the discriminator is fed synthetic data
points created by the generator.

In the IGN setting, instead, the two models will work in the opposite order:
the discriminator $d$ processes input data into a latent tensor, and the
generator $g$ creates a synthetic data point from that latent tensor. Our
idempotent generative network is then $F=g\circ d$.

##  3 Conditioning idempotent generative networks

Idempotent generative networks are then both simple and elegant from a
theoretical perspective, as well as powerful from a practical one: in
particular IGNs can generate data in a single step, without the multi-steps
denoising process implemented by diffusion models or the ‘single-token-
generation’ approach of autoregressive text generation models.

However, the approach taken in [10] leaves one crucial question unanswered:
can we condition the model to generate new samples with prescribed
characteristics?

We present multiple approaches to answer affirmatively the previous question.

###  3.1 Theoretical perspective

Let $X$ be again the feature space and $D$ be the data manifold, immersed into
$X$. Suppose that every data point in $D$ comes with a _condition_. This could
for example be a label or a caption (in case $D$ represent images). Let $C$ be
the ‘condition’ space, the space where conditions $c$ belong to.

We can now extend the paradigm behind idempotent generative network to
consider conditioning: we have an augmented data manifold $\widetilde{D}$
immersed in $X\times C$, every point $(d,c)$ consisting of a data point $d$
and its condition $c$. We want to find an idempotent map

$F:X\times C\longrightarrow X\times C$

such that:

  * •

(reconstruction) $F(d,c)=(d,c)$ for every augmented data point
$(d,c)\in\widetilde{D}$.

  * •

(idempotence) $F(F(x,c))=F(x,c)$ for every $(x,c)$ in $X\times C$.

  * •

(tightness) the range of $F$ is _as small as possible_ , i.e. $F(X\times C)$
is a low-dimensional submanifold of $X\times C$.

As explained in [10], the _tightness_ condition is crucial to avoid the model
learning to replicate the identity function (which satisfies the _idempotence_
and the _reconstruction_ requirements).

We plan on using the model as follows: suppose we want to generate a sample
with condition $c$. We randomly select a noisy sample $z\in X$, we calculate
$F(z,c)\in X\times C$ and we extract the first component $\tilde{d}$ of
$F(z,c)$.

###  3.2 Loss function

Even a perfect model $F$ satisfying the three conditions in section 3.1 above
is not guaranteed to return a data point $\tilde{d}$ related to condition $c$.
Indeed, the idempotent function $F$ is not being explicitly constrained to
return a data point $\tilde{d}$ related to condition $c$, just a data point in
the extended data manifold.

In order to help the model learn the conditioning, we treat _mismatched_ data
points like noisy examples. A _mismatched_ data point is a pair $(d,\hat{c})$
where $(d,c)$ belongs to the augmented data manifold $\widetilde{D}$ and
$c\neq\hat{c}$. In other words, we tell the model that $(d,\hat{c})$ should
not be part of the range. By treating mismatched data points like noisy
samples, we teach the network to _not_ act like the identity function on them.

Let $D_{c}=\left\\{d\in D\,|\,(d,c)\in\widetilde{D}\right\\}$ be the subset of
the data manifold $D$ consisting of points with condition $c$. We are then
teaching the model $F$ to have range contained in $\widetilde{D}$ (tightness
and reconstruction conditions), but that $D_{c}\times{\hat{c}}$ should not be
part of the range for any $\hat{c}\neq c$ (mismatched conditions). By
exclusion this forces the model $F$ to further restrict its range to the
unions of $D_{c}\times{c}$ as $c$ varies in $C$.

The loss function for conditional idempotent generative networks is then made
of five terms:

  1. 1.

(reconstruction loss $L_{rec}$) Real data samples $(d,c)$ from the augmented
data manifold $\widetilde{D}$ should remain unchanged by the model:
$F(d,c)=(d,c)$.

  2. 2.

(idempotent on noisy samples $L_{idem}^{noise}$) For any noisy sample $z$ and
any condition $c$, any further application of the model beyond the first one
should act like the identity function: $F(F(z,c))=F(z,c)$.

  3. 3.

(tightness on noisy samples $L_{tight}^{noise}$) The idempotence condition on
noisy samples and the reconstruction condition should be satisfied while
ensuring $F$ has as small a range as possible.

  4. 4.

(idempotent on mismatched samples $L_{idem}^{mism}$) For a data point $(d,c)$
on the augmented data manifolds $\widetilde{D}$ and a condition $\hat{c}$
different from $c$, any further application of the model beyond the first one
should act like the identity function: $F(F(d,\hat{c}))=F(d,\hat{c})$.

  5. 5.

(tightness on mismatched samples $L_{tight}^{mism}$) The idempotence condition
on mismatched samples and the reconstruction condition should be satisfied
while ensuring $F$ has as small a range as possible.

Each one of these conditions is implemented via the choice of a distance
function $\delta$ on the space $X\times C$, which is used to compare a data
point before and after the application of $F$.222For simplicity one usually
picks the same distance function $\delta$ for all loss terms, although this is
not theoretically necessary. For example the reconstruction loss is
implemented as

$L_{rec}(d,c)=\delta\left((d,c),F(d,c)\right)\quad\forall(d,c)\in\widetilde{D}.$

The total loss is then obtained as a weighted sum of these five terms:

$L=w_{rec}\cdot L_{rec}+w_{idem}^{noise}\cdot
L_{idem}^{noise}+w_{tight}^{noise}\cdot L_{tight}^{noise}+w_{idem}^{mism}\cdot
L_{idem}^{mism}+w_{tight}^{mism}\cdot L_{tight}^{mism}.$

Each one of the five weights is a hyperparameter that can be tuned during
training.

During training, every batch will then be composed of three types of data
points:

  * •

true data points $(d,c)$ in the augmented data manifold $\widetilde{D}$.

  * •

noisy samples $(z,c)$ where $z$ is a randomly sampled from $X$ following a
prior distribution.

  * •

mismatched samples $(d,\hat{c})$ for $(d,c)\in\widetilde{D}$ and $\hat{c}\neq
c$.

The relative proportion of the three types of data points within a batch is a
training hyperparameter.

###  3.3 Metrics

We discuss some theoretical considerations regarding what metrics should be
used to evaluate Conditional Generative Idempotent Networks.

At a high level, an appropriate metric should measure the following two
dimensions:

  * •

(realism) The synthetic data point $F(z,c)$ generated from condition $c$
should indeed look like a sample from $D_{c}$.

  * •

(diversity) Synthetic data points $F(z,c)$ generated from the same condition
$c$ as we vary the noisy input $z$ should space over as large a subset of $X$
as possible.

There is a fairly standard set of metrics typically used to measure these two
dimensions for Generative Adversarial Networks, for example Inception Score
(IS) and Fréchet Inception Distance (FID). These metrics have been adapted to
the setup of conditional generative adversarial networks (see for example
[1]).

More complex metrics (see for example [4]) have been suggested for conditional
image generation models, especially if the generation is implemented via a
diffusion process. In this latter scenario, human evaluation is most of the
time still the most accurate performance metric.

All of these metrics are reasonable choices and the one to be preferred mostly
depends on the specific use case. If the condition space $C$ is a finite set
of classes, metrics based on the inception score and the Fréchet inception
distance are probably adequate. If the condition space $C$ is a continuous
manifold (for example, $c$ can be any English sentence), then human evaluation
is likely to be preferred.

##  4 Two proposed implementations of CIGN

In this section, we described two possible implementations of the Conditional
Idempotent Generative Networks.

For both implementations, we assume that the data has a channel dimension (for
example audio, images or video data).

Let $(x,c)$ be an element in $X\times C$. We denote by $(n,S)$ the shape of
the tensor $x$, where $n$ is the channel dimension and
$S=(s_{1},\ldots,s_{k})$ represents any number of additional dimensions (for
example, if $x$ is an image then $S=(h,w)$ is a pair of height and width
dimensions).

###  4.1 Channel conditioning

Let $E:C\longrightarrow\mathbb{R}^{S}$ be an embedding function, where we
denote by
$\mathbb{R}^{S}=\mathbb{R}^{s_{1}}\times\ldots\times\mathbb{R}^{s_{k}}$ the
subspace corresponding to the non-channel dimensions of the feature space $X$.

_Channel conditioning_ is implemented by concatenating $E(c)$ and $x$ along
the channel dimension.

This is a very generic and flexible idea, with many different variations
available:

  * •

We can build different embeddings $E_{i}(c)$ whose target space consists of
the non-channel dimensions of the input $x_{i}$ to layer $i$ of the network,
and then concatenate $E_{i}(c)$ to $x_{i}$.

  * •

We can build embeddings whose target space is multiple copies of
$\mathbb{R}^{S}$, so that when concatenating along the channel diemnsion we
give the model additional flexibility on learning from the condition.

In section 5.3 we describe our experiments on training a conditional
idempotent generative network with channel conditioning on the MNIST dataset.

###  4.2 Filter conditioning

Let $S^{\prime}=(s_{1}^{\prime},\ldots,s_{k}^{\prime})$ be any $k$-tuple with
$0<s_{i}$ for any $1\leq i\leq k$. Let
$E:C\longrightarrow\mathbb{R}^{S^{\prime}}$ be an embedding function.

_Filter conditioning_ is implemented by calculating cross correlation (or
transposed cross correlation) between $x$ and $E(c)$ along the non-channel
dimensions. This is inspired by the concept of correlation filters for object
detection (see for example [2]).

One can think of filter conditioning as replacing a standard or transposed
convolutional layer (where the model learns the kernel’s weights and those
weights are shared across all inputs) by a simpler cross-correlation layer,
where each input is (transpose) cross-correlated to its condition’s embedding,
and the model learns the embedding layer’s weights.

This also is a very generic and flexible idea, with many different variations
available:

  * •

We can choose different filter dimensions $S^{\prime}$ for the embedding’s
target space - in fact, $S^{\prime}$ can be considered as a hyperparameter of
the model.

  * •

We can replace any number of convolutional layers (of a base model
architecture that we start with) with filter conditioning layers where the
filter dimensions $S^{\prime}$ are the dimensions of the original
convolutional layer.

  * •

We can concatenate the result of a standard convolutional layer with the
result of a filter conditioning layer, assuming that two layers’ parameters
(kernel’s shape, padding, stride, dilation) are identical.

##  5 Experiments on the MNIST dataset

In this section we discuss our experiments implementing CIGNs on the MNIST
dataset. All experiments were run on consumer hardware - specifically one
single free GPU offered by Sagemaker Studio Lab, with a maximum session time
of 4 hours.

While this severely limited our experiments (in terms of the model size and in
terms of the number of hyperparameter runs we experimented with), we are still
able to successfully train well-performing Conditional Idempotent Generative
Networks (CIGN), able to generate different good quality images from the
_same_ noisy data point when prompted with different conditions.

###  5.1 Preprocessing steps and high-level architecture

The MNIST dataset from [12] consists of images with a single channel and a
size of $28\times 28$ pixels, valued in $[0,1]$. We preprocess the dataset by
transforming images into torch tensors, and then linearly rescaling pixels to
be valued between $-1$ and $1$. We add some small random noise to the training
data, to fight overfitting.

In the setup of the previous sections, we have then the feature space

$X=\mathbb{R}^{1}\times\mathbb{R}^{28}\times\mathbb{R}^{28},$

with the condition space

$C=\\{0,...,9\\}$

being the label set of the MNIST dataset.

Noisy data points are sampled from $X$ according to a normal distribution
centered at $0$ and with standard deviation $1$, and then clipped so that they
belong to $\left[-1,1\right]^{1\times 28\times 28}$.

The distance function $\delta$ used to implement the loss function is

$\delta\big{(}(x_{1},c_{1}),(x_{2},c_{2})\big{)}=L^{1}\big{(}x_{1}-x_{2}\big{)}=\left|x_{1}-x_{2}\right|_{1},$

the $L^{1}$ norm of the difference between the feature space components. While
this is not strictly speaking a distance on $X\times C$, it works for our
purposes since the models we implement pass through the condition $c$ ‘as-is’
(meaning that $F(x,c)=\left(F_{1}(x,c),c\right)$), and for the purpose of loss
calculation we only need to calculate the distance between $(x,c)$ and
$F(x,c)$.

At a high level, the architecture of our Conditional Idempotent Generative
Networks is composed of:

  * •

a discriminator with domain $X\times C$, outputting a pair of a latent tensor
and the condition $c$;

  * •

a generator, with input the latent tensor and the condition $c$, and range
$X\times C$.

The CIGN is then obtained by running the discriminator first and the generator
second.

The architectures of the discriminator and the generator are heavily inspired
by the DCGAN framework, in its pytorch implementation ([6]).

The discriminator is defined as a sequence of convolutional layers with:

  * •

an increasing number of channels, starting from a latent dimension $l_{dim}$
and doubling each layer, with the final latent tensor having size
$(i_{dim},1,1)$.333We denote by $i_{dim}$ the intermediate dimension. We call
this _intermediate_ because this is obtained exactly halfway through the
forward pass of the entire CIGN, when the latent tensor is outputted by the
discriminator and passed to the generator.

  * •

decreasing non-channel dimensions, until the latent tensor has non-channel
dimensions $(1,1)$ \- that is to say, it is a 1-dimensional tensor.

The generator is defined by mirroring exactly the layers of the discriminator,
but in the opposite order: a sequence of transpose convolutional layers with
decreasing number of channels and increasing non-channel dimensions. The
parameters of each convolutional layer of the discriminator (kernel size,
stide, padding, dilation) are equal to those of the corresponding transpose
convolutional layer of the generator.

###  5.2 Metrics

Our main metric is the average Fréchet Inception Distance across all 10
classes. We remark that this is identical to the metric named WCFID in section
3.2 of [1].

For each class, we compare the 1000 samples available in the validation
dataset from [12] to 1000 synthetic samples generated by the model under
consideration.

This comparison is implemented via the pytorch-fid library (see [7]): this
library leverages the ImageNet-v3 model to get intermediate embeddings of the
images, which are used to estimate the parameters of the data distribution
under a normality assumption.

The aforementioned implementation [7] allows the user to choose among four
different embeddings (each one passing the image through a partial
inception-v3 model [11] up to a certain layer), but one of them is not
available to us due to the small size of the distributions we are comparing
(just 1000 samples). The other three embeddings output respectively a
$64$-dimensional tensor, a $192$-dimensional tensor and a $768$-dimensional
tensor. We denote the corresponding metrics by $\mathrm{FID}_{64}$,
$\mathrm{FID}_{192}$ and $\mathrm{FID}_{768}$.

###  5.3 Channel conditioning

We implement the channel conditioning approach from section 4.1 in its
simplest form:

  1. 1.

we train one shared embedding layer

$E:C\longrightarrow\mathbb{R}^{\mathrm{emb_{dim}}},$

  2. 2.

and then for each (possibly transpose) convolutional layer where the input
$x_{i}$ has non-channel dimensions $(h_{i},w_{i})$, we train one linear layer

$L_{i}:\mathbb{R}^{\mathrm{emb_{dim}}}\longrightarrow\mathbb{R}^{h_{i}\cdot
w_{i}},$

and we concatenate $L_{i}(E(c))$ to $x_{i}$ along the channel dimension,
before the forward pass of the (possibly transpose) convolutional layer.

####  5.3.1 Training details

We run two channel conditioning experiments where the models have different
hyperparameters, as in table 1.

Table 1: Parameters for channel conditioning experiments Model size | Embedding dim | Latent dim | Intermediate dim | Num parameters
---|---|---|---|---
small | 5 | 32 | 128 | 1.55m
large | 5 | 64 | 128 | 5.6m

For each real data point in the augmented data manifold, we create nine
mismatched data points by pairing the image with all other labels except the
correct one. The number of noisy images in each batch is equal to the number
of real data points. Our proportion of real:noise:mismatched data points is
then 1:1:9.

We train the model for 50 epochs, with a learning rate of $0.001$ and a batch
size of $128$ (this means that each backward pass is actually calculated on
$128\times 11=1408$ data points, due to the 1:1:9 ratio mentioned in the
previous paragraph).

We experiment with a few different choices of weights for the five summands of
the loss function described in 3.2, and ultimately we choose the following
values:

  * •

$w_{rec}=20.0$

  * •

$w_{idem}^{noise}=20.0$

  * •

$w_{tight}^{noise}=2.5$

  * •

$w_{idem}^{mism}=3.0$

  * •

$w_{tight}^{mism}=1.0$

The detailed architecture of generator and discriminator are described in
appendix A.1, including figures 1 and 2.

###  5.4 Filter conditioning

We implement the filter conditioning approach from section 4.2 in its simplest
form:

  1. 1.

we train one shared embedding layer

$E:C\longrightarrow\mathbb{R}^{\mathrm{emb_{dim}}},$

  2. 2.

for each (possibly transpose) convolutional layer where the standard DCGAN
architecture in [6] has input $x_{i}$ with $c_{i}$ channels and a kernel of
size $(h_{i},w_{i})$, we train one linear layer

$L_{i}:\mathbb{R}^{\mathrm{emb_{dim}}}\longrightarrow\mathbb{R}^{c_{i}\cdot
h_{i}\cdot w_{i}},$

  3. 3.

we calculate channel-wise (possibly transpose) cross-correlation between
$x_{i}$ and $L_{i}(E(c))$,

  4. 4.

for each (possibly transpose) convolutional layer where the standard DCGAN
architecture has $c_{in}$ input channels and $c_{out}$ output channels, we
train a _channel mixer_ linear layer

$C_{i}:\mathbb{R}^{c_{in}}\longrightarrow\mathbb{R}^{c_{out}},$

which we apply to the output of the cross-correlation from the previous step.

####  5.4.1 Training details

We run two filter conditioning experiments where the models have different
hyperparameters, as in table 2.

Table 2: Parameters for channel conditioning experiments Model size | Embedding dim | Latent dim | Intermediate dim | Num parameters
---|---|---|---|---
small | 5 | 64 | 128 | 616k
large | 5 | 96 | 256 | 1.38m

For each real data point in the augmented data manifold, we create nine
mismatched data points by pairing the image with all other labels except the
correct one. The number of noisy images in each batch is equal to the number
of real data points. Our proportion of real:noise:mismatched data points is
then 1:1:9.

We train the model for 50 epochs, with a learning rate of $0.0001$ and a batch
size of $64$ (this means that each backward pass is actually calculated on
$64\times 11=704$ data points, due to the 1:1:9 ratio mentioned in the
previous paragraph).

We experiment with a few different choices of weights for the five summands of
the loss function described in 3.2, and ultimately we choose the following
values:

  * •

$w_{rec}=20.0$

  * •

$w_{idem}^{noise}=20.0$

  * •

$w_{tight}^{noise}=2.5$

  * •

$w_{idem}^{mism}=8.0$

  * •

$w_{tight}^{mism}=1.0$

The detailed architecture of generator and discriminator is described in
appendix A.2, including figures 3 and 4.

###  5.5 Results

In appendices B.1 and B.2 we display samples of the images generated by the
best experiments.

In this section we detail the metrics calculated on our trained experiments.

As discussed in section 5.2, our metric is the average Fréchet Inception
Distance (FID) across the ten classes. For each trained experiment and each
class, we compare the $1000$ samples on the MNIST validation dataset with
$1000$ samples generated by the trained experiment. We also report in table 3
minimum and maximum FID across all ten classes.

Table 3: FID metrics for all experiments |  | filter_large | filter_small | channel_large | channel_small
---|---|---|---|---|---
$\mathrm{FID}_{64}$ | mean | 0.8278 | 0.3601 | 0.1149 | 0.6586
min | 0.0837 | 0.0959 | 0.0455 | 0.1126
max | 3.25 | 0.7516 | 0.272 | 2.1216
$\mathrm{FID}_{192}$ | mean | 3.3095 | 1.792 | 0.6598 | 2.7955
min | 0.4912 | 0.5992 | 0.3735 | 0.5908
max | 12.5165 | 5.1655 | 1.1886 | 8.2649
$\mathrm{FID}_{768}$ | mean | 0.4423 | 0.3384 | 0.2852 | 0.4619
min | 0.2485 | 0.1869 | 0.1576 | 0.3557
max | 0.8765 | 0.4641 | 0.5184 | 0.7127

Essentially all metrics point to the large experiment trained with channel
conditioning as the best experiment of the four. On the other hand, that
experiment has about 4x as many trainable parameters as the next two
experiments (the small experiment trained with channel-conditioning and the
large experiment trained with filter conditioning).

The metrics for the small experiment with filter conditioning (the smallest
model of all, with only 616k trainable parameters) are clearly better than the
small experiment with channel conditioning.

The metrics for the large experiment with filter conditioning are unusually
off. Human review of the synthetic images generated by this model suggests
that it tends to create ‘thick font’ images (see appendix B.2) which is likely
the reason for the large FID values and therefore poor results. We have not
investigated why this experiment tends to generate ‘thick font’ images.

In conclusion, we believe that the jury is still out regarding which of the
conditioning mechanism is better - the best performing experiment is the large
model with channel conditioning, but the small filter conditioning model
clearly outperforms a slightly larger model trained with channel conditioning.

## References

  * [1] Benny, Y., Galanti, T., Benaim, S., & Wolf, L. (2021). Evaluation metrics for conditional image generation. International Journal of Computer Vision, 129, 1712-1731.
  * [2] Bertinetto, L., Valmadre, J., Henriques, J. F., Vedaldi, A., & Torr, P. H. (2016). Fully-convolutional siamese networks for object tracking. In Computer Vision-ECCV 2016 Workshops: Amsterdam, The Netherlands, October 8-10 and 15-16, 2016, Proceedings, Part II 14 (pp. 850-865). Springer International Publishing.
  * [3] Engel, J., Agrawal, K. K., Chen, S., Gulrajani, I., Donahue, C., Roberts, A. (2019). Gansynth: Adversarial neural audio synthesis. arXiv preprint arXiv:1902.08710.
  * [4] Ku, M., Li, T., Zhang, K., Lu, Y., Fu, X., Zhuang, W., & Chen, W. (2023). Imagenhub: Standardizing the evaluation of conditional image generation models. arXiv preprint arXiv:2310.01596.
  * [5] Netron web application at https://netron.app/
  * [6] Official Pytorch implementation of [9]. https://github.com/pytorch/examples/tree/main/dcgan
  * [7] Pytorch implementation of Fréchet inception distance, available at https://github.com/mseitzer/pytorch-fid
  * [8] Pulfer, B. implementation of IGN from original paper [10], available at https://github.com/brianpulfer/idempotent-generative-network
  * [9] Radford, A., Metz, L., Chintala, S. (2015). Unsupervised representation learning with deep convolutional generative adversarial networks. arXiv preprint arXiv:1511.06434.
  * [10] Shocher A, Dravid A, Gandelsman Y, Mosseri I, Rubinstein M, Efros AA. Idempotent generative network. arXiv preprint arXiv:2311.01462. 2023 Nov 2.
  * [11] Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., … & Rabinovich, A. (2015). Going deeper with convolutions. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 1-9).
  * [12] MNIST dataset from torchvision, available at https://pytorch.org/vision/stable/generated/torchvision.datasets.MNIST.html

##  Appendix A Model architectures

In this appendix, we describe in details the model architectures of our
experiments from section 5. As before, $l_{dim}$ is the latent dimension and
$i_{dim}$ is the intermediate dimension.

We also report figures displaying the architectures of our experiment (all
these figures were obtained with Netron [5]).

###  A.1 Channel conditioning architectures

The discriminator is defined as five consecutive convolutional layers. Each
convolutional layer $i$ consists of the following five steps:

  1. 1.

batch normalization (except for the initial layer)

  2. 2.

concatenation of the input tensor and the embedding of the condition
$L_{i}(E(c))$ (passed through a hyperbolic tangent activation) along the
channel dimension

  3. 3.

convolutional layer

  4. 4.

dropout layer with rate $0.15$

  5. 5.

activation function (LeakyReLU with slope $0.2$ except for the last layer
which uses the sigmoid function)

Table 4 describes the specific parameters of each convolutional layer.

Table 4: Convolutional layer parameters for discriminator with channel conditioning Layer | Input size | Output size | Kernel size | Stride | Padding
---|---|---|---|---|---
1 | (1, 28, 28) | ($l_{dim}$, 14, 14) | (4, 4) | 2 | 1
2 | ($l_{dim}$, 14, 14) | ($l_{dim}*2$, 7, 7) | (4, 4) | 2 | 1
3 | ($l_{dim}*2$, 7, 7) | ($l_{dim}*4$, 4, 4) | (3, 3) | 2 | 1
4 | ($l_{dim}*4$, 4, 4) | ($l_{dim}*8$, 2, 2) | (4, 4) | 2 | 1
5 | ($l_{dim}*8$, 2, 2) | ($i_{dim}$, 1, 1) | (2, 2) | 1 | 0

The generator is defined as five consecutive transpose convolutional layers.
Each transpose convolutional layer consists of the following five steps:

  1. 1.

batch normalization

  2. 2.

concatenation of the input tensor and the embedding of the condition, along
the channel dimension

  3. 3.

transpose convolutional layer

  4. 4.

dropout layer with rate $0.15$ (except the last layer)

  5. 5.

activation function (ReLU except for the last layer which uses the hyperbolic
tangent function)

Table 5 describes the specific parameters of each transpose convolutional
layer.

Table 5: Convolutional layer parameters for generator with channel conditioning Layer | Input size | Output size | Kernel size | Stride | Padding
---|---|---|---|---|---
1 | ($i_{dim}$, 1, 1) | ($l_{dim}*8$, 2, 2) | (2, 2) | 1 | 0
2 | ($l_{dim}*8$, 2, 2) | ($l_{dim}*4$, 4, 4) | (4, 4) | 2 | 1
3 | ($l_{dim}*4$, 4, 4) | ($l_{dim}*2$, 7, 7) | (3, 3) | 2 | 1
4 | ($l_{dim}*2$, 7, 7) | ($l_{dim}$, 14, 14) | (4, 4) | 2 | 1
5 | ($l_{dim}$, 14, 14) | (1, 28, 28) | (4, 4) | 2 | 1
Figure 1: Architecture of the generator with channel conditioning Figure 2:
Architecture of the discriminator with channel conditioning

###  A.2 Filter conditioning architectures

The discriminator is defined as five consecutive ‘convolutional’ layers. Each
convolutional layer consists of the following five steps:

  1. 1.

batch normalization (except for the initial layer)

  2. 2.

channel-wise cross-correlation between the input tensor $x_{i}$ and the
embedding of the condition $L_{i}(E(c))$ (passed through a hyperbolic tangent
activation)

  3. 3.

channel mixer layer $C_{i}$, applied on the channel dimension

  4. 4.

dropout layer with rate $0.15$

  5. 5.

activation function (LeakyReLU with slope $0.2$ except for the last layer
which uses the sigmoid function)

Table 6 describes the specific parameters of each convolutional layer.

Table 6: Convolutional layer parameters for discriminator with filter conditioning Layer | Input size | Kernel $L_{i}(E(c))$ size | Stride | Padding
---|---|---|---|---
1 | (1, 28, 28) | (1, 4, 4) | 2 | 1
2 | ($l_{dim}$, 14, 14) | ($l_{dim}$, 4, 4) | 2 | 1
3 | ($l_{dim}*2$, 7, 7) | ($l_{dim}*2$, 3, 3) | 2 | 1
4 | ($l_{dim}*4$, 4, 4) | ($l_{dim}*4$, 4, 4) | 2 | 1
5 | ($l_{dim}*8$, 2, 2) | ($l_{dim}*8$, 2, 2) | 1 | 0

The generator is defined as five consecutive ‘transpose convolutional’ layers.
Each transpose convolutional layer consists of the following five steps:

  1. 1.

batch normalization

  2. 2.

channel-wise transpose cross-correlation between the input tensor $x_{i}$ and
the embedding of the condition $L_{i}(E(c))$ (passed through a hyperbolic
tangent activation)

  3. 3.

channel mixer layer $C_{i}$, applied on the channel dimension

  4. 4.

dropout layer with rate $0.15$ (except the last layer)

  5. 5.

activation function (ReLU except for the last layer which uses the hyperbolic
tangent function)

Table 7 describes the specific parameters of each transpose convolutional
layer.

Table 7: Convolutional layer parameters for generator with filter conditioning Layer | Input size | Kernel $L_{i}(E(c))$ size | Stride | Padding
---|---|---|---|---
1 | ($i_{dim}$, 1, 1) | ($i_{dim}$, 2, 2) | 1 | 0
2 | ($l_{dim}*8$, 2, 2) | ($l_{dim}*8$, 4, 4) | 2 | 1
3 | ($l_{dim}*4$, 4, 4) | ($l_{dim}*4$, 3, 3) | 2 | 1
4 | ($l_{dim}*2$, 7, 7) | ($l_{dim}*2$, 4, 4) | 2 | 1
5 | ($l_{dim}$, 14, 14) | ($l_{dim}$, 4, 4) | 2 | 1
Figure 3: Architecture of the generator with filter conditioning Figure 4:
Architecture of the discriminator with filter conditioning

##  Appendix B Generated images

In this section, we display synthetic images generated by our experiments on
the MNIST dataset. Each group of image is generated by one experiment,
starting from the same noisy sample but with different conditions.

###  B.1 Channel conditioning experiments

The following images were generated by the best experiments trained with the
channel conditioning mechanism described in section 5.3.

Figure 5: Images generated by small model with channel conditioning Figure 6:
Images generated by large model with channel conditioning

###  B.2 Filter conditioning experiments

The following images were generated by the best experiments trained with the
filter conditioning mechanism described in section 5.4.

Figure 7: Images generated by small model with filter conditioning Figure 8:
Images generated by large model with filter conditioning


11institutetext: Department of Mathematics,
University of Cape Coast, Ghana.
11email<EMAIL_ADDRESS>
# MULTIPATCH DISCONTINUOUS GALERKIN IGA FOR THE BIHARMONIC PROBLEM ON SURFACES

Stephen E. Moore

###### Abstract

We present the analysis of interior penalty discontinuous Galerkin
Isogeometric Analysis (dGIGA) for the biharmonic problem on orientable
surfaces $\Omega\subset\mathbb{R}^{3}.$ Here, we consider a surface consisting
of several non-overlapping patches as typical in multipatch dGIGA. Due to the
non-overlapping nature of the patches, we construct NURBS approximation spaces
which are discontinuous across the patch interfaces via a penalty scheme. By
an appropriate discrete norm, we present a priori error estimates for the non-
symmetric, symmetric and semi-symmetric interior penalty methods. Finally, we
confirm our theoritical results with numerical experiments.

###### Keywords:

discontinuous Galerkin methods, surface biharmonic equation, isogeometric
analysis, a priori error estimates, lapalce-beltrami.

##  1 Introduction

Partial Differential Equations (PDEs) defined on surfaces embedded in
$\mathbb{R}^{3}$ arise in many fields of application including material
science, fluid mechanics, electromagnetics, biology and image processing.
Fourth-order partial differential equations (PDEs) are particularly important
in several areas of applied mechanics, the theory of elasticity, mechanics of
elastic plates, and the slow flow of viscous fluids, see e.g. [23]. Some
examples of physical flows modelled by fourth order PDE include fluids on the
lungs [8], ice formation [18], imaging [14], designing special curves on
surfaces [9], modeling of interfaces in multiphase fluid flows and the
modeling of surface active agents (surfactants), see e.g. [24, 19].

In this article, we consider the fourth-order boundary value problem: find
$u:\overline{\Omega}\rightarrow\mathbb{R}$ such that

$\displaystyle\Delta_{\Omega}^{2}u+u=f\quad\text{in}\quad\Omega,$ (1.1)

where $f$ is a square integrable load vector defined on a compact smooth and
orientable surface $\Omega\subset\mathbb{R}^{3}$ with boundary
$\partial\Omega$ consisting of Dirichlet $\Gamma_{D}$ and Neumann $\Gamma_{N}$
boundaries, i.e. $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ with the following
considtions

$\displaystyle
u=g_{0},\quad\mathbf{n}\cdot\nabla_{\Omega}u=g_{1},\,\,\text{on}\,\,\Gamma_{D},\quad\Delta_{\Omega}u=g_{2},\quad$
$\displaystyle\mathbf{n}\cdot\nabla_{\Omega}\Delta_{\Omega}u=g_{3}\,\,\text{
on }\,\,\Gamma_{N},$ (1.2)

where $\mathbf{n}$ is the outward directed normal vector to the boundary
$\partial\Omega.$ We define the bi-Laplacian operator
$\Delta_{\Omega}^{2}:=\Delta_{\Omega}\Delta_{\Omega}$ with $\Delta_{\Omega}$
as the Laplace-Beltrami operator, and the boundary data $g_{0},g_{1},g_{2}$
and $g_{3}$ are smooth functions.

Surface PDEs are usually treated with surface FEM which is a very popular
discretization method. However, it is well known that surface FEM has major
drawbacks due to the discrete variational formulation of the PDE that is
constructed on a triangulated surface which contains the finite elements space
[7]. Several works concerning the discretization of fourth order PDEs on
surfaces using surface FEM including an application of interior penalty
Galerkin (IPG) methods have been presented, see e.g. [13]; which is an
extension of the second order PDEs version [6].

Fourth order PDE require continuously differentiably $(i.e.C^{1}-)$ piece-wise
polynomial basis functions which are known to be practically difficult to
construct as well as computationally expensive. However, in recent times, a
new approximation method has been proposed that has $(p-1)-$continuously
differentiable basis functions i.e. $C^{(p-1)},$ with degree $p\geq 1$ which
makes it ideal towards the approximation of higher order PDEs including the
biharmonic problem. This method is known as the isogeometric analysis (IGA).
Moreover, IGA uses the same class of basis functions for both representing the
geometry of the domain and approximating the solution of the PDEs [2, 22].

Multipatch discontinuous Galerkin IGA has been introduced and analyzed for
second order elliptic problems on surfaces with matching and non-matching
meshes, see e.g. [12, 11, 20, 17]. Here, the computational domain consists of
several conforming non-overlapping subdomains. By applying interior penalty
methodology, we construct discrete spaces on the patches allowing for
discontinuity along the patch interfaces. In this regards, the results
presented in this article is an extention of the dGIGA to the biharmonic
problem [16]. In this article, we will present a priori error estimate for
multipatch discontinuous Galerkin isogeometric analysis (dGIGA) for biharmonic
problem on conforming patches with matching meshes on orientable surfaces.

We organize the article as follows; In Section 2, the function spaces, weak
formulation and the isogeometric analysis framework, NURBS surfaces,
geometrical mappings and isogeometric analysis. The derivation of the interior
penalty discontinuous Galerkin scheme is presented in Section 3. Then, in
Section 4, we present a discrete NURBS space $V_{h},$ the discrete norms
$\|\cdot\|_{h}$ and $\|\cdot\|_{h,*}$ then subsequently prove the coercivity
of the bilinear forms. The boundedness of the bilinear forms is asserted in a
product space $V_{h,*}\times V_{h},$ where we will need another discrete norm
$\|\cdot\|_{h,*}$ defined on the vector space $V_{h,*}.$ By using the idea of
equivalence of norms, we are able to present coercivity and boundedness on the
norm $\|\cdot\|_{h,*}.$ The error analysis of the dGIGA scheme is presented in
Section 5. In Section 6, we present and discuss numerical experiments to
confirm our theoretical results. Finally, we draw some conclusions and discuss
future works in Section Conclusion.

##  2 Preliminaries

Let the computational domain $\Omega$ be a compact smooth and oriented surface
with boundary $\partial\Omega.$ We introduce the Sobolev space
$H^{s}(\Omega):=\\{v\in L_{2}(\Omega)\,:\,D^{\alpha}_{\Omega}v\in
L_{2}(\Omega),\,\,\text{for}\,\,0\leq|\alpha|\leq s\\}$, where $L_{2}(\Omega)$
denote the space of square integrable functions and let
$\alpha=(\alpha_{1},\ldots,\alpha_{d})$ be a multi-index with non-negative
integers $\alpha_{1},\ldots,\alpha_{d}$,
$|\alpha|=\alpha_{1}+\ldots+\alpha_{d},$
$D^{\alpha}_{\Omega}:=\partial^{|\alpha|}/\partial x^{\alpha},$ and associate
with the sobolev space $H^{s}(\Omega)$ the norm
$\|v\|_{H^{s}(\Omega)}=\left(\sum_{0\leq|\alpha|\leq
s}\|D^{\alpha}_{\Omega}v\|^{2}_{L_{2}(\Omega)}\right)^{1/2}$ see, e.g. [1].

The weak variational formulation of the biharmonic problem (1.1) reads: find
$u\in V_{D}$ such that

$\displaystyle a(u,v)$ $\displaystyle=\ell(v),\quad\forall v\in V_{0},$ (2.1)

where the bilinear and linear forms are given by

$\displaystyle a(u,v)$
$\displaystyle=\int_{\Omega}\Delta_{\Omega}u\Delta_{\Omega}v\,dx+uv\,dx\quad\text{and}\quad$
$\displaystyle\ell(v)$
$\displaystyle=\int_{\Omega}fv\,dx+\int_{\Gamma_{N}}(\mathbf{n}\cdot\nabla_{\Omega}vg_{2}+vg_{3})ds\,$
(2.2)

and the hyperplane and test space given by $V_{D}:=\\{v\in
H^{2}(\Omega):v|_{\Gamma_{D}}=g_{0},\quad\mathbf{n}\cdot\nabla_{\Omega}v|_{\Gamma_{D}}=g_{1}\\}$
and $V_{0}:=\\{v\in
H^{2}(\Omega):v|_{\Gamma_{D}}=0,\quad\mathbf{n}\cdot\nabla_{\Omega}v|_{\Gamma_{D}}=0\\}.$

The existence and uniqueness of the variational problem (2.1) follows the
well-known Lax-Milgram lemma see e.g. [4].

###  2.1 NURBS Geometrical Mapping and Surfaces

For the parameter domain $\widehat{\Omega},$ we define a vector-valued
independent variable in the parameter domain $\widehat{\Omega}$ by
$\xi=(\xi_{1},\xi_{2})\in\mathbb{R}^{2}.$ By means of a smooth and invertible
geometrical mapping $\mathrm{\Phi}$, the computational domain
$\Omega\subset\mathbb{R}^{3}$ is defined as

$\mathrm{\Phi}:\widehat{\Omega}\rightarrow\Omega\subset\mathbb{R}^{3},\quad\mathbf{\xi}\rightarrow
x=\mathrm{\Phi}(\mathbf{\xi}),$ (2.3)

where $\widehat{\Omega}\subset\mathbb{R}^{2}$ is the parameter domain as
illustrated in Figure 1.

Figure 1: Illustration of the geometrical map
$\Phi:\widehat{\Omega}\subset\mathbb{R}^{2}\rightarrow\Omega\subset\mathbb{R}^{3}$
for a patch.

We introduce briefly some important mathematical tools necessary for the
analysis of surface PDEs. The following objects are obtained by means of the
geometrical mapping (2.3) in the parameter domain. The Jacobian $\widehat{J},$
first fundamental form $\widehat{F}$ and the determinant $\widehat{g}$ of the
geometrical mapping are respectively given by

$\displaystyle\widehat{J}$
$\displaystyle:=\left[\frac{\partial\mathrm{\Phi}_{k}}{\partial\mathbf{\xi}_{l}}\right]\in\mathbb{R}^{3\times
2},\quad k=1,2,3,\,l=1,2,$ (2.4) $\displaystyle\widehat{F}(\xi)$
$\displaystyle=\left(\widehat{J}(\xi)\right)^{T}\widehat{J}(\xi)\in\mathbb{R}^{3\times
3}\quad\text{and}\quad\widehat{g}(\xi)=\sqrt{\det\left({\widehat{F}(\xi)}\right)}\in\mathbb{R}.$
(2.5)

Next, we present some differential operators by using notations in the
parameter domain. We consider a smooth function $\phi$ defined on the manifold
$\Omega,$ by using the invertible geometrical mapping (2.3) to obtain

$\phi(x)=\widehat{\phi}(\xi)\circ\Phi^{-1}(x),\quad x\in\Omega,$ (2.6)

where $\widehat{\phi}(\xi)=\phi(\mathrm{\Phi}(\xi)).$ Using the gradient
operator in the parameter space $\nabla\widehat{\phi}$, the tangential
gradient of the manifold is given by

$\nabla_{\Omega}\phi(x):=\widehat{J}(\mathbf{\xi})\widehat{F}^{-1}(\mathbf{\xi})\nabla\widehat{\phi}(\mathbf{\xi})\circ\Phi^{-1}(x).$
(2.7)

The divergence operator for the vector-valued function can be written as

$\nabla_{\Omega}\cdot\varphi(x):=\dfrac{1}{\widehat{g}\mathbf{\xi})}\nabla\cdot\left[\widehat{g}(\mathbf{\xi})\widehat{F}^{-1}(\xi)\widehat{J}^{T}(\xi)\widehat{\varphi}(\mathbf{\xi})\right]\circ\Phi^{-1}(x).$
(2.8)

The Laplace-Beltrami operator on the manifold $\Omega$ is defined for a twice
continuously differentiable function $\phi:\Omega\rightarrow\mathbb{R}$ as

$\Delta_{\Omega}\phi(x)=\dfrac{1}{\widehat{g}(\mathbf{\xi})}\nabla\cdot\left[\widehat{g}(\mathbf{\xi})\widehat{F}^{-1}(\mathbf{\xi})\nabla\widehat{\phi}(\mathbf{\xi})\right]\circ\Phi^{-1}(x).$
(2.9)

The surface gradient of the Laplace-Beltrami operator on the manifold $\Omega$
is defined for a thrice continuously differentiable function
$\phi:\Omega\rightarrow\mathbb{R}$ as

$\nabla_{\Omega}\Delta_{\Omega}\phi(x)=\left[\widehat{J}(\mathbf{\xi})\widehat{F}^{-1}(\mathbf{\xi})\nabla\left(\dfrac{1}{\widehat{g}(\mathbf{\xi})}\nabla\cdot\left[\widehat{g}(\mathbf{\xi})\widehat{F}^{-1}(\mathbf{\xi})\nabla\widehat{\phi}(\mathbf{\xi})\right]\right)\right]\circ\Phi^{-1}(x).$
(2.10)

The unit normal vector on the manifold $\Omega\subset\mathbb{R}^{3}$ is
obtained by the geometrical mapping of

$\widehat{\mathbf{n}}(\xi):=\dfrac{\widehat{t}_{1}(\xi)\times\widehat{t}_{2}(\xi)}{\|\widehat{t}_{1}(\xi)\times\widehat{t}_{2}(\xi)\|},$
(2.11)

where $\widehat{t}_{l}(\xi):=\partial\Phi(\xi)/\partial\xi_{l}$ is the tangent
vector to a curve in $\mathbb{R}^{3}$ with $l=1,2.$ The manifold $\Omega$ has
a tangent plane at $\xi$ if the tangent vectors are linearly independent.

Finally, by means of the geometrical mapping (2.3), we can write the Jacobian,
first fundamental form and the determinant on the computational domain
$\Omega$ as follows

$\displaystyle J(x)=\widehat{J}(\xi)\circ\mathrm{\Phi}^{-1}(x),\quad
F(x)=\widehat{F}(\xi)\circ\mathrm{\Phi}^{-1}(x)\quad\text{and}\quad
g(x)=\widehat{g}(\xi)\circ\mathrm{\Phi}^{-1}(x).$ (2.12)

###  2.2 B-Spline, NURBS and Isogeometric Analysis

For a comprehensive understanding of isogeometric analysis, we refer the
reader to [5] and the references therein. However, for the purpose of
completion, we present in this article briefly some vital information
necessary for the formulation and discussion of multipatch isogeometric
analysis. For positive integers $p$ and $n,$ let us define a vector
$\mathrm{\Xi}=\left\\{0=\xi_{1},\ldots,\xi_{n+p+1}=1\right\\}$ with a non-
decreasing sequence of real numbers in the parameter domain
$\widehat{\Omega}=[0,1]$ called a knot vector on the unit interval
$\widehat{\Omega}=[0,1].$ Given $\mathrm{\Xi}$ with $p\geq 1$ and $n$ as the
number of basis functions, the univariate B-spline basis functions are defined
by the Cox -de Boor recursion formula

$\displaystyle\widehat{B}_{i,0}(\xi)$ $\displaystyle=\left\\{\begin{aligned}
&1&\text{if}\quad&\xi_{i}\leq\xi<\xi_{i+1},\\\ &0&\text{else},&\\\
\end{aligned}\right.$ $\displaystyle\widehat{B}_{i,p}(\xi)$
$\displaystyle=\frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}\widehat{B}_{i,p-1}(\xi)+\frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}\widehat{B}_{i+1,p-1}(\xi),$
(2.13)

where any division by zero is defined to be zero. We note that a basis
function of degree $p$ is $(p-m)$ times continuously differentiable across a
knot value with the multiplicity $m$. For example, if all internal knots have
the multiplicity $m=1$, then B-splines of degree $p$ are globally
$(p-1)-$continuously differentiable.

In general for higher-dimensional problems, the B-spline basis functions are
tensor products of the univariate B-spline basis functions. We define tesor
product basis functions as follows: let
$\mathrm{\Xi}_{\alpha}=\left\\{\xi_{1,\alpha},\ldots,\xi_{n_{\alpha}+p_{\alpha}+1,\alpha}\right\\}$
be the knot vectors for every direction $\alpha=1,\ldots,d.$ Let
$\mathbf{i}:=(i_{1},\ldots,i_{d}),\mathbf{p}:=(p_{1},\ldots,p_{d})$ and the
set
$\overline{\mathcal{I}}=\\{\mathbf{i}=(i_{1},\ldots,i_{d}):i_{\alpha}=1,2,\ldots,n_{\alpha};\;\alpha=1,2,\ldots,d\\}$
be multi-indicies. Then the tensor product B-spline basis functions are
defined by

$\widehat{B}_{\mathbf{i},\mathbf{p}}(\xi):=\prod\limits_{\alpha=1}^{d}\widehat{B}_{i_{\alpha},p_{\alpha}}(\xi_{\alpha}),$
(2.14)

where $\xi=(\xi_{1},\ldots,\xi_{d})\in\widehat{\Omega}=(0,1)^{d}.$ The
univariate and multivariate B-spline basis functions are defined in the
parametric domain by means of the corresponding B-spline basis functions
$\\{\widehat{B}_{\mathbf{i},\mathbf{p}}\\}_{\mathbf{i}\in\overline{\mathcal{I}}}.$

The distinct values $\xi_{i},i=1,\ldots,n$ of the knot vectors $\mathrm{\Xi}$
provides a partition of $(0,1)^{d}$ creating a mesh
$\widehat{\mathcal{K}}_{h}$ in the parameter domain where $\widehat{K}$ is a
mesh element. The computational domain is described by means of a geometrical
mapping $\mathbf{\Phi}$ such that $\Omega=\mathbf{\Phi}(\widehat{\Omega})$ and

$\displaystyle\mathbf{\Phi}(\xi):=\sum_{\mathbf{i}\in\overline{\mathcal{I}}}C_{\mathbf{i}}\widehat{B}_{\mathbf{i},\mathbf{p}}(\xi),$
(2.15)

where $C_{\mathbf{i}}$ are the control points. Next, we describe NURBS bassis
functions. These basis functions are usually prefered in industry due to their
ability to exactly represent most shapes and particularly conic families. The
NURBS basis functions are obtained from the B-spline basis functions by means
of a geometrical mapping as follows

$\displaystyle\mathbf{\Phi}(\xi):=\sum_{i=1}^{n}\widehat{N}_{\mathbf{i},\mathbf{p}}(\xi)\mathbf{P}_{i}\quad\text{and}\quad\widehat{N}_{\mathbf{i},\mathbf{p}}(\xi):=\dfrac{w_{i}}{\sum_{\mathbf{j}=1}^{n}\widehat{B}_{\mathbf{j},\mathbf{p}}(\xi)}\widehat{B}_{\mathbf{i},\mathbf{p}}(\xi),$
(2.16)

where $\mathbf{P}_{i},i=1,\ldots,n$ are the control points in the physical
space and $\widehat{N}_{\mathbf{i},\mathbf{p}}(\xi)$ are the NURBS basis
functions obtained by projective transformation of the B-spline basis
functions with weight $w_{i}\in\mathbb{R}$ and $n$ the number of basis
functions. The basis functions in the computational domain are defined by
means of the geometrical mapping as
$N_{\mathbf{i},\mathbf{p}}:=\widehat{N}_{\mathbf{i},\mathbf{p}}\circ\mathbf{\Phi}^{-1}$
and the discrete function spaces by

$\mathbb{V}_{h}=\text{span}\\{N_{\mathbf{i},\mathbf{p}}:\mathbf{i}\in\overline{\mathcal{I}}\\}.$
(2.17)

In several real life applications, the computational domain $\Omega$ is
usually decomposed into $N$ non-overlapping sub-domains $\Omega_{i}$ called
patches denoted by $\mathcal{T}_{h}:=\\{\Omega_{i}\\}_{i=1}^{N}$ such that
$\overline{\Omega}=\bigcup_{i=1}^{N}\overline{\Omega}_{i}$ and
$\Omega_{i}\cap\Omega_{j}=\emptyset$ for $i\neq j.$ Each patch is the image of
an associated geometrical mapping $\mathbf{\Phi}_{i}$ such that
$\mathbf{\Phi}_{i}(\widehat{\Omega})=\Omega_{i},i=1,\ldots,N,$ see Figure 2.

Figure 2: Illustration of the multi-patch isogeometric analysis map
$\Phi_{i}(\widehat{\Omega})=\Omega_{i}$ and
$\Phi_{j}(\widehat{\Omega})=\Omega_{j},i\neq j.$

We denote by $F_{ij}=\partial\Omega_{i}\cap\partial\Omega_{j},i\neq j,$ the
interior facets of two patches as illustrated in Figure 3. The collection of
all such interior facets is denoted by $\mathcal{F}_{I},$ the set of Dirichlet
facets $F_{i}$ on $\Gamma_{D}$ by $\mathcal{F}_{D}$ and the set of Neumann
facets $F_{i}$ on $\Gamma_{N}$ by $\mathcal{F}_{N}.$ Furthermore, the
collection of all internal, Dirichlet and Neumann facets is denoted by
$\mathcal{F}:=\mathcal{F}_{I}\cup\mathcal{F}_{D}\cup\mathcal{F}_{N}.$

Figure 3: Illustration of the underlying mesh of the multi-patch isogeometric

We assume that for each patch $\Omega_{i},i=1,\ldots,N,$ the underlying mesh
$\mathcal{K}_{h,i}$ is quasi-uniform i.e.

$\displaystyle h_{K}\leq h_{i}\leq C_{u}h_{K},\quad\text{for all}\quad
K\in\mathcal{K}_{h,i},\quad i=1,\ldots,N,$ (2.18)

where $C_{u}\geq 1$ and $h_{i}=\max\\{h_{K},K\in\mathcal{K}_{h,i}\\}$ is the
mesh size of $\Omega_{i}$ and $h_{K}$ is the diameter of of the mesh element
$K$.

##  3 Interior Penalty Variational Formulation

We introduce function spaces necessary for the derivation of interior penalty
Galerkin schemes. Also, we assign to each patch $\Omega_{i}$ an integer
$s_{i}$ and collect them in the vector $\mathbf{s}=\\{s_{1},\ldots,s_{N}\\}.$
We define the broken Sobolev space

$H^{\mathbf{s}}(\Omega,\mathcal{T}_{h}):=\\{v\in
L_{2}(\Omega):v|_{\Omega_{i}}\in H^{s_{i}}(\Omega_{i}),\;i=1,\ldots,N\\},$
(3.1)

and the corresponding broken Sobolev norm and semi-norm

$\displaystyle\|v\|_{H^{\mathbf{s}}(\Omega,\mathcal{T}_{h})}:=\left(\sum_{i=1}^{N}\|v\|^{2}_{H^{s_{i}}(\Omega_{i})}\right)^{1/2}\quad\text{and}\quad|v|_{H^{\mathbf{s}}(\Omega,\mathcal{T}_{h})}:=\left(\sum_{i=1}^{N}|v|^{2}_{H^{s_{i}}(\Omega_{i})}\right)^{1/2},$
(3.2)

respectively. We define the jump and average of the normal derivatives across
the interior facets $F_{ij}\in\mathcal{F}_{I}$ of $v\in
H^{2}(\Omega,\mathcal{T}_{h})$ by

$\displaystyle\llbracket\nabla_{\Omega}v\rrbracket:=\mathbf{n}_{i}\cdot\nabla_{\Omega}v_{i}+\mathbf{n}_{j}\cdot\nabla_{\Omega}v_{j},\quad
and\quad\\{\nabla_{\Omega}v\\}:=\frac{1}{2}\left(\mathbf{n}_{i}\cdot\nabla_{\Omega}v_{i}+\mathbf{n}_{j}\cdot\nabla_{\Omega}v_{j}\right),$
(3.3)

whereas the jump and average functions on the Dirichlet facets
$F_{i}\in\mathcal{F}_{D}$ are given by

$\displaystyle\llbracket\nabla_{\Omega}v\rrbracket:=\mathbf{n}_{i}\cdot\nabla_{\Omega}v_{i},\quad
and\quad\\{\nabla_{\Omega}v\\}:=\mathbf{n}_{i}\cdot\nabla_{\Omega}v_{i}.$
(3.4)

The following equality on the interior facets $F_{ij}\in\mathcal{F}_{I}$ is
obtained by using the definitions of jumps and averages as

$\displaystyle\llbracket{uv}\rrbracket$
$\displaystyle=\\{a\\}\llbracket{b}\rrbracket+\\{b\\}\llbracket{a}\rrbracket,\quad\forall
a,b\in\mathbb{R}.$ (3.5)

The interior penalty variational scheme reads: find $u\in
H^{4}(\Omega,\mathcal{T}_{h})$ such that,

$a_{h}(u,v)=\ell_{h}(v),\quad\forall v\in H^{4}(\Omega,\mathcal{T}_{h}),$
(3.6)

where the bilinear form is given by

$\displaystyle a_{h}(u,v)$
$\displaystyle=\sum_{i=1}^{N}\int_{\Omega_{i}}(\Delta_{\Omega}u\Delta_{\Omega}v+uv)\,dx$
$\displaystyle-\sum_{F\in\mathcal{F}}\int_{F}\llbracket\nabla_{\Omega}v\rrbracket\\{\Delta_{\Omega}u\\}\,ds-\beta_{0}\sum_{F\in\mathcal{F}}\int_{F}\llbracket\nabla_{\Omega}u\rrbracket\\{\Delta_{\Omega}v\\}\,ds$
$\displaystyle+\sum_{F\in\mathcal{F}}\int_{F}\\{\nabla_{\Omega}\Delta_{\Omega}u\\}\llbracket
v\rrbracket\,ds+\beta_{1}\sum_{F\in\mathcal{F}}\int_{F}\\{\nabla_{\Omega}\Delta_{\Omega}v\\}\llbracket
u\rrbracket\,ds$
$\displaystyle+\sum_{F\in\mathcal{F}}\frac{\delta_{1}}{h_{i}^{3}}\int_{F}\llbracket
u\rrbracket\llbracket
v\rrbracket\,ds+\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\int_{F}\llbracket\nabla_{\Omega}u\rrbracket\llbracket\nabla_{\Omega}v\rrbracket\,ds,$
(3.7)

and the linear form reads as

$\displaystyle\ell_{h}(v)$
$\displaystyle=\sum_{i=1}^{N}\int_{\Omega_{i}}fv\,dx+\sum_{F_{i}\in\mathcal{F}_{D}}\int_{F_{i}}(\delta_{1}v+\beta_{1}\mathbf{n}_{i}\cdot\nabla_{\Omega}\Delta_{\Omega}v)g_{0}\,ds$
$\displaystyle+\sum_{F_{i}\in\mathcal{F}_{D}}\int_{F_{i}}(\delta_{0}\mathbf{n}_{i}\cdot\nabla_{\Omega}v+\beta_{0}\Delta_{\Omega}v)g_{1}\,ds+\sum_{F_{i}\in\mathcal{F}_{N}}\int_{F_{i}}(\mathbf{n}_{i}\cdot\nabla_{\Omega}v\,g_{2}+vg_{3})\,ds.$
(3.8)

The parameters $\beta_{0},\beta_{1}\in\\{-1,1\\}$ and determine the interior
penalty Galerkin (IPG) scheme. Here, we will describe the four main schemes a
follows;

  1. 1.

$\beta_{0}=\beta_{1}=-1$ is the non-symmetric interior penalty Galerkin (NIPG)

  2. 2.

$\beta_{0}=\beta_{1}=1$ is the symmetric interior penalty Galerkin (SIGP)

  3. 3.

$\beta_{0}=-1,\beta_{1}=1$ is the semi-symmetric interior penalty Galerkin
(SSIGP1)

  4. 4.

$\beta_{0}=1,\beta_{1}=-1$ is the semi-symmetric interior penalty Galerkin
(SSIGP2)

###### Remark 1

Although, there is no current literature on the choice of the penalty
parameters for the bilinear form (3), we choose
$\delta_{0}=\delta_{1}=(p+1)(p+d)/d$ where $p$ is the degree of the NURBS and
$d$ is the dimension of the computational surface i.e. $d=3$ and
$\Omega\subset\mathbb{R}^{3}.$ This choice of the penalty parameters yielded
accurate simulations for the NIPG case presented in [16].

For the weak continuity of the fluxes on the interior facets, the exact
solution $u$ must satisfy the following

$\displaystyle\llbracket
u\rrbracket=0,\llbracket\nabla_{\Omega}u\rrbracket=0,\llbracket\Delta_{\Omega}u\rrbracket=0,\llbracket\nabla_{\Omega}\Delta_{\Omega}u\rrbracket=0,\quad\forall
F\in\mathcal{F}.$

This enables us to show that the interior penalty Galerkin scheme is
consistent, i.e. if $u\in H^{4}(\Omega)$ is the solution of (2.1). Then $u$ is
the solution to dGIGA variational identity (3.6), see [16] and [21, Lemma 3].

##  4 Analysis of the dGIGA Scheme

In this section, we consider the existence and uniqueness of the bilinear form
in the discrete setting. Thus, we require discrete spaces for the whole
computational domain $\Omega.$ Since the domain consists of several subdomains
or patches, we associate with each subdomain a discrete space as follows: Let
us consider the B-Spline space $V_{h}\subset H^{4}(\Omega,\mathcal{T}_{h})$
defined as

$\displaystyle V_{h}:=\\{v\in
L_{2}(\Omega):v|_{\Omega_{i}}\in\mathbb{V}_{h,i},i=1,\ldots,N\\},$ (4.1)

where the B-Spline space $\mathbb{V}_{h,i}$ corresponds to the patch
$\Omega_{i},i=1,\ldots,N$ for B-splines of degree $p_{i}\geq 3.$ The discrete
dGIGA scheme then reads as: find $u_{h}\in V_{h}$ such that

$\displaystyle a_{h}(u_{h},v_{h})=\ell(v_{h}),\quad\forall v_{h}\in V_{h}.$
(4.2)

An immediate consequence of the consistency results as discussed in the latter
part of section 3 is the Galerkin orthogonality property, i.e.,

$\displaystyle a_{h}(u-u_{h},v_{h})=0,\quad\forall v_{h}\in V_{h}.$ (4.3)

Next, for $v\in H^{2}(\Omega,\mathcal{T}_{h}),$ we show that the bilinear form
$a_{h}(\cdot,\cdot),$ is coercive with respect to the following norm

$\displaystyle\|v\|_{h}^{2}:=\sum_{i=1}^{N}\|\Delta_{\Omega}v\|^{2}_{L_{2}(\Omega_{i})}+\sum_{i=1}^{N}\|v\|_{L_{2}(\Omega_{i})}^{2}+\sum_{F\in\mathcal{F}}\frac{\delta_{1}}{h_{i}^{3}}\|\llbracket
v\rrbracket\|^{2}_{L_{2}(F)}+\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\|\llbracket\nabla_{\Omega}v\rrbracket\|^{2}_{L_{2}(F)}.$
(4.4)

###### Remark 2

For some function $v\in H^{2}(\Omega,\mathcal{T}_{h}),$ if $\|v\|_{h}=0,$ then

$\displaystyle\Delta_{\Omega}v$
$\displaystyle=0\quad\text{in}\quad\Omega_{i}\quad\text{for all}\quad
i=1,\ldots,N$ $\displaystyle\llbracket v\rrbracket$
$\displaystyle=0,\llbracket\nabla_{\Omega}v\rrbracket=0\quad on\quad
F\in\mathcal{F}.$ (4.5)

Thus using the theory of elliptic interface problems $v=0$ in the whole domain
$\overline{\Omega}$ since it is a weak solution to (2).

The analysis of the dGIGA scheme requires the patch-wise inverse and trace
inequalities given by the following lemmata, see [15, chapter 2].

###### Lemma 1

Let $K_{i}\in\mathcal{K}_{h,i},i=1,\ldots,N$. Then the inverse inequalities,

$\displaystyle\|\partial^{j}v\|_{L_{2}(\Omega_{i})}$ $\displaystyle\leq
C_{inv,1,u}h_{i}^{-1}\|\partial^{j-1}v\|_{L_{2}(\Omega_{i})},$ (4.6)
$\displaystyle\|v\|_{L_{2}(\partial\Omega_{i})}$ $\displaystyle\leq
C_{inv,0,u}h_{i}^{-1/2}\|v\|_{L_{2}(\Omega_{i})},$ (4.7)

hold for all $j\geq 1,v\in V_{h},$ where $C_{inv,1,u}$ and $C_{inv,0,u}$ are
positive constants, which are independent of $h_{i}$ and $\Omega_{i}$.

###### Lemma 2

Let $a_{h}(\cdot,\cdot):V_{h}\times V_{h}\rightarrow\mathbb{R}$ be the
discrete bilinear form defined in (4.2) with
$\beta_{0},\beta_{1}\in\\{-1,1\\}$ and let $c_{0}$ and $c_{1}$ be nonnegative
constants to be determined in the proof. Assume $\sigma_{0}>0$ and
$\sigma_{1}>0$ and suppose that $\delta_{0}\geq c_{0}$ and $\delta_{1}\geq
c_{1}.$ Then, there exists a positive constant $\mu_{c}$ such that

$\displaystyle a_{h}(v_{h},v_{h})\geq\mu_{c}\|v_{h}\|_{h}^{2},\quad\forall
v_{h}\in V_{h}.$ (4.8)

###### Proof

By setting $u_{h}=v_{h}$ in (3), we proceed as follows

$\displaystyle a_{h}(v_{h},v_{h})$
$\displaystyle=\sum_{i=1}^{N}\int_{\Omega_{i}}(\Delta_{\Omega}v_{h})^{2}\,dx+\sum_{F\in\mathcal{F}}\frac{\delta_{1}}{h_{i}^{3}}\int_{F}\llbracket
v_{h}\rrbracket^{2}\,ds+\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\int_{F}\llbracket\nabla_{\Omega}v_{h}\rrbracket^{2}\,ds$
$\displaystyle-(1+\beta_{0})\sum_{F\in\mathcal{F}}\int_{F}\llbracket\nabla_{\Omega}v_{h}\rrbracket\\{\Delta_{\Omega}v_{h}\\}\,ds+(1+\beta_{1})\sum_{F\in\mathcal{F}}\int_{F}\llbracket
v_{h}\rrbracket\\{\nabla_{\Omega}\Delta_{\Omega}v_{h}\\}\,ds$
$\displaystyle=\sum_{i=1}^{N}\|\Delta_{\Omega}v_{h}\|^{2}_{L_{2}(\Omega_{i})}+\sum_{F\in\mathcal{F}}\frac{\delta_{1}}{h_{i}^{3}}\|\llbracket
v_{h}\rrbracket\|^{2}_{L_{2}(F)}+\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\|\llbracket\nabla_{\Omega}v_{h}\rrbracket\|^{2}_{L_{2}(F)}$
$\displaystyle-(1+\beta_{0})\sum_{F\in\mathcal{F}}\|\llbracket\nabla_{\Omega}v_{h}\rrbracket\|_{L_{2}(F)}\|\\{\Delta_{\Omega}v_{h}\\}\|_{L_{2}(F)}\,$
$\displaystyle+(1+\beta_{1})\sum_{F\in\mathcal{F}}\|\llbracket
v_{h}\rrbracket\|_{L_{2}(F)}\|\\{\nabla_{\Omega}\Delta_{\Omega}v_{h}\\}\|_{L_{2}(F)}.$
(4.9)

Using Cauchy-Schwarz’s inequality, Lemma 1 and Young’s inequality, we obtain

$\displaystyle(1+\beta_{0})\sum_{F\in\mathcal{F}}\|\llbracket\nabla_{\Omega}v_{h}\rrbracket\|_{L_{2}(F)}\|\\{\Delta_{\Omega}v_{h}\\}\|_{L_{2}(F)}$
$\displaystyle\leq(1+\beta_{0})\left(\sum_{F\in\mathcal{F}}\dfrac{C^{2}_{inv,0}}{2h_{i}\sigma_{0}}\|\llbracket\nabla_{\Omega}v_{h}\rrbracket\|_{L_{2}(F)}^{2}+\dfrac{\sigma_{0}}{2}\sum_{i=1}^{N}\|\Delta_{\Omega}v_{h}\|_{L_{2}(\Omega_{i})}^{2}\right).$
(4.10) $\displaystyle(1+\beta_{1})\sum_{F\in\mathcal{F}}\|\llbracket
v_{h}\rrbracket\|_{L_{2}(F)}\|\\{\nabla_{\Omega}\Delta_{\Omega}v_{h}\\}\|_{L_{2}(F)}$
$\displaystyle\leq(1+\beta_{1})\left(\sum_{F\in\mathcal{F}}\dfrac{C^{2}_{inv,0}}{2h_{i}^{3}\sigma_{1}}\|\llbracket
v_{h}\rrbracket\|_{L_{2}(F)}^{2}+\dfrac{\sigma_{1}}{2}\sum_{i=1}^{N}\|\Delta_{\Omega}v_{h}\|_{L_{2}(\Omega_{i})}^{2}\right).$
(4.11)

Substituting (4.10) and (4.11) into (4) yields

$\displaystyle a_{h}(v_{h},v_{h})$
$\displaystyle\geq\sum_{i=1}^{N}\bigg{(}1-\frac{\sigma_{0}}{2}(1+\beta_{0})-\frac{\sigma_{1}}{2}(1+\beta_{1})\bigg{)}\|\Delta_{\Omega}\|_{L_{2}(\Omega_{i})}^{2}$
$\displaystyle+\sum_{F\in\mathcal{F}}\bigg{(}\delta_{0}-\frac{C^{2}_{inv,0,1}}{2\sigma_{0}}(1+\beta_{0})\bigg{)}\|\llbracket\nabla_{\Omega}v_{h}\rrbracket\|_{L_{2}(F)}^{2}$
$\displaystyle+\sum_{F\in\mathcal{F}}\bigg{(}\delta_{1}-\frac{C^{2}_{inv,0}}{2\sigma_{1}}(1+\beta_{1})\bigg{)}\|\llbracket
v_{h}\rrbracket\|_{L_{2}(F)}^{2}.$ (4.12)

The ellipticity constant is given by

$\mu_{c}=1-\frac{\sigma_{0}}{2}(1+\beta_{0})-\frac{\sigma_{1}}{2}(1+\beta_{1}),$

and determined by the choices of $\sigma_{0}$ and $\sigma_{1}$ as well as
$\beta_{0}$ and $\beta_{1}.$ The last two terms (4) hold for $\delta_{0}\geq
c_{0}$ and $\delta_{1}\geq c_{1}$ where

$c_{0}=\frac{(1+\beta_{0})}{2\sigma_{0}}C^{2}_{inv,0,1}\quad\text{and}\quad
c_{1}=\frac{(1+\beta_{1})}{2\sigma_{1}}C^{2}_{inv,0},$

which completes the proof. ∎

Finally, we obtain the coercivity of the bilinear forms and the corresponding
penalty parameters $\delta_{0}$ and $\delta_{1}$ in the next theorem.

###### Theorem 4.1

Let $a_{h}(\cdot,\cdot)$ be the discrete bilinear form defined in (4.2) with
positive constants $\delta_{0}$ and $\delta_{1}$ and let $c_{0}$ and $c_{1}$
be nonnegative constants as in Lemma 2. Suppose that for Then there exists a
positive constant $\mu_{c}$ such that

$\displaystyle a_{h}(v_{h},v_{h})\geq\mu_{c}\|v_{h}\|_{h}^{2},\quad\forall
v_{h}\in V_{h}.$

###### Proof

The proof follows from Lemma 2 and the method is determined by the choice of
$\beta_{0}$ and $\beta_{1}.$ ∎

From Theorem 4.1, we obtain the uniqueness of the solution of the discrete
variational problem (4.2). Since the discrete variational problem is in the
finite dimensional space $V_{h},$ the uniqueness therefore yields the
existence of the solution $u_{h}\in V_{h}$ of (4.2).

###### Lemma 3

Let $K\in\mathcal{K}_{h,i},i=1,\ldots,N$ and $\widehat{K}=\Phi_{i}^{-1}(K)$.
Then the scaled trace inequality

$\displaystyle\|v\|_{L_{2}(\partial\Omega_{i})}\leq
C_{tr,u}h_{i}^{-1/2}\left(\|v\|_{L_{2}(\Omega_{i})}+h_{i}|v|_{H^{1}(\Omega_{i})}\right),$
(4.13)

holds for all $v\in H^{1}(\Omega_{i}),$ where $h_{i}$ denotes the global mesh
size of patch $\Omega_{i}$ in the physical domain, and $C_{tr,u}$ is a
positive constant that only depends on the quasi-uniformity and shape
regularity of the mapping $\Phi_{i}$.

To enable us derive uniform boundedness of the bilinear form
$a_{h}(\cdot,\cdot):V_{h,*}\times V_{h}\rightarrow\mathbb{R},$ where
$V_{h,*}:=V_{D}\cap H^{\mathbf{s}}(\Omega,\mathcal{T}_{h})+V_{h}$ with
$\mathbf{s}\geq 4$ and equipped with the norm

$\displaystyle\|v\|_{h,*}^{2}=\|v\|_{h}^{2}+\sum_{F\in\mathcal{F}}\frac{h_{i}^{3}}{\delta_{1}}\|\\{\nabla_{\Omega}\Delta_{\Omega}v\\}\|^{2}_{L_{2}(F)}+\sum_{F\in\mathcal{F}}\frac{h_{i}}{\delta_{0}}\|\\{\Delta_{\Omega}v\\}\|^{2}_{L_{2}(F)}.$
(4.14)

Indeed, the norm (4.4) is also a norm on $H^{4}(\Omega,\mathcal{T}_{h})$ since
$H^{4}(\Omega,\mathcal{T}_{h})\subset H^{2}(\Omega,\mathcal{T}_{h}).$

###### Lemma 4

Let $a_{h}(\cdot,\cdot):V_{h,*}\times V_{h}$ be the bilinear form defined in
(3) with $\beta_{0},\beta_{1}\in\\{-1,1\\}$ and $\delta_{0},\delta_{1}>0.$
Then there exists a positive constant $\mu_{b},$ such that

$\displaystyle|a_{h}(u,v_{h})|\leq\mu_{b}\|u\|_{h,*}\|v_{h}\|_{h},\quad\forall
u\in V_{h,*},v_{h}\in V_{h}.$ (4.15)

###### Proof

The proof follows by using the he Cauchy-Schwarz inequality to estimate the
terms in the bilinear form (3). However, concerning the concerning the third
term, we apply the inverse inequality (4.7) for $v_{h}\in V_{h}$ as follows

$\displaystyle\bigg{|}\sum_{F\in\mathcal{F}}\int_{F}\\{\Delta_{\Omega}v_{h}\\}\llbracket\nabla_{\Omega}u\rrbracket\,ds\bigg{|}$
$\displaystyle\leq\left(\sum_{i=1}^{N}\frac{C_{inv,0,u}^{2}}{\delta_{0}}\|\Delta_{\Omega}v_{h}\|^{2}_{L_{2}(\Omega_{i})}\right)^{\frac{1}{2}}\left(\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\|\llbracket\nabla_{\Omega}u\rrbracket\|^{2}_{L_{2}(F)}\right)^{\frac{1}{2}}.$

The fifth term is also estimated by using the inverse inequalities (4.6) and
(4.7) for $v_{h}\in V_{h},$ to obtain

$\displaystyle\bigg{|}\sum_{F\in\mathcal{F}}\int_{F}\\{\nabla_{\Omega}\Delta_{\Omega}v_{h}\\}\llbracket
u\rrbracket\,ds\bigg{|}$
$\displaystyle\leq\left(\sum_{i=1}^{N}\frac{C_{inv,0,1}^{2}}{\delta_{1}}\|\Delta_{\Omega}v_{h}\|^{2}_{L_{2}(\Omega_{i})}\right)^{\frac{1}{2}}\left(\sum_{F\in\mathcal{F}}\int_{F}\frac{\delta_{1}}{h_{i}^{3}}\|\llbracket
u\rrbracket\|^{2}_{L_{2}(F)}\right)^{\frac{1}{2}},$ (4.16)

where $C_{inv,0,1}^{2}=C_{inv,0,u}^{2}C_{inv,1,u}^{2}.$ By puttting all the
terms together and using Cauchy-Schwarz’s inequality, we complete the proof
with the boundedness constant given as
$\mu_{b}=2\sqrt{\max\\{1,(1+C_{inv,0,1}^{2}|\beta_{1}|/\delta_{1}+C^{2}_{inv,0,u}|\beta_{0}|/\delta_{0})\\}}.$
∎

It is also possible to show the existence and uniqueness results for the norm
$\|\cdot\|_{h,*}$ due to the uniform equivalence of norms on $V_{h}.$

###### Lemma 5

The norms $\|\cdot\|_{h}$ and $\|\cdot\|_{h,*}$ are uniformly equivalent on
the discrete space $V_{h}$ such that

$C^{*}\|v_{h}\|_{h,*}\leq\|v_{h}\|_{h}\leq\|v_{h}\|_{h,*},\quad\forall
v_{h}\in V_{h},$ (4.17)

where $C^{*}$ is a mesh independent positive constant.

###### Proof

The upper bound follows immediately. The lower bound is derived by applying
the inverse inequalities of Lemma 1 where we coplete the proof with
$C^{*}=\left(1+C_{inv,0,1}^{2}/\delta_{0}+C_{inv,0,u}^{2}/\delta_{1}\right)^{-1}$
and $C_{inv,0,1}^{2}=C_{inv,0,u}^{2}C_{inv,1,u}^{2}$. ∎

Due to Lemma 4.17, we can derive coercivity and boundedness results for the
bilinear form $a_{h}(\cdot,\cdot)$ in the norm $\|\cdot\|_{h,*}.$ By using the
results from Theorem 4.1, the coercivity on $V_{h}$ yields

$a_{h}(v_{h},v_{h})\geq\widetilde{\mu}_{c}\|v_{h}\|_{h,*}^{2},\quad\forall
v_{h}\in V_{h},$ (4.18)

where $\widetilde{\mu}_{c}=C^{*}\mu_{c}$ is a nonnegative constant independent
of $h.$ Also, the boundedness of the bilinear form following from Lemma 4 is
given by

$\displaystyle|a_{h}(u_{h},v_{h})|\leq\tilde{\mu}_{b}\|u_{h}\|_{h,*}\|v_{h}\|_{h,*},\quad\forall(u,v_{h})\in
V_{h,*}\times V_{h},$ (4.19)

with $\tilde{\mu}_{b}$ independent of $h.$

##  5 Error Estimates for dGIGA discretization scheme

To obtain a priori error estimates in both $\|\cdot\|_{h}$ and
$\|\cdot\|_{h,*}$ norms, we require approximation estimates by means of a
quasi-intepolant. We denote by $\Pi_{h}:L_{2}(\Omega_{i})\rightarrow V_{h,i}$
such a quasi-interpolant that yields optimal approximation results for each
patch or subdomain $\Omega_{i},i=1,\ldots,N,$ see [2, 3] for the proof.

###### Lemma 6

Let $l_{i}$ and $s_{i}$ be integers with $0\leq l_{i}\leq s_{i}\leq p_{i}+1$
and $K\in\mathcal{K}_{h,i}$. Then there exist an interpolant $\Pi_{h}v\in
V_{h,i}$ for all $v\in H^{s_{i}}(\Omega_{i})$ and a constant $C_{s}>0$ such
that the following inequality holds

$\displaystyle\sum_{K\in\mathcal{K}_{h,i}}|v-\Pi_{h}v|^{2}_{H^{l}(K)}$
$\displaystyle\leq
C_{s}h_{i}^{2(s_{i}-l_{i})}\|v\|^{2}_{H^{s_{i}}(\Omega_{i})},$ (5.1)

where $h_{i}$ is the mesh size in the physical domain, and $p$ denotes the
underlying polynomial degree of the B-spline or NURBS.

If the multiplicity of the inner knots is not larger than $p_{i}+1-l_{i}$ and
$\Pi_{h}v\in V_{h}\cap H^{l_{i}}(\Omega_{i})$ for each patch
$\Omega_{i},i=1,\ldots,N,$ then the local estimate (5.1) yields a global
estimate.

###### Proposition 1

Let us assume that the multiplicity of the inner knots is not larger than
$p_{i}+1-l_{i}.$ Given the integers $l_{i}$ and $s_{i}$ such that $0\leq
l_{i}\leq s_{i}\leq p_{i}+1,$ there exist a positive constant $C_{s}$ such
that for a function $v\in H^{s_{i}}(\Omega_{i})$

$|v-\Pi_{h}v|_{H^{l_{i}}(\Omega_{i})}\leq
C_{s}h_{i}^{(s_{i}-l_{i})}\|v\|_{H^{s_{i}}(\Omega_{i})},$ (5.2)

where $h_{i}$ denotes the maximum mesh-size parameter in the physical domain
and the generic constant $C_{s}$ only depends on $l_{i},s_{i}$ and $p_{i}$,
the shape regularity of the physical domain $\Omega_{i}$ described by the
mapping $\Phi$ and, in particular, $\nabla_{\Omega}\mathrm{\Phi}.$

###### Proof

See [22, Proposition 3.2]. ∎

We consider that the quasi-interpolant is the same for each patch, i.e.
$\Pi_{h}:V_{D}\cap H^{\mathbf{s}}(\Omega,\mathcal{T}_{h})\rightarrow V_{h}$
with $\mathbf{s}\geq 4.$

###### Lemma 7

Let $v\in V_{D}\cap H^{\mathbf{s}}(\Omega,\mathcal{T}_{h})$ with a positive
integer $\mathbf{s}\geq 4$ and let $F\in\mathcal{F}_{I}\cup\mathcal{F}_{D}$ be
the facets. Also, let $\delta_{0}$ and $\delta_{1}$ be chosen as in Theorem
4.1. By assuming quasi-uniform meshes, then there exists a quasi-interpolant
$\Pi_{h}$ such that $\Pi_{h}v\in V_{h}$ and the following estimates hold;

$\displaystyle\|\nabla_{\Omega}^{q}(v-\Pi_{h}v)\|_{L_{2}(\partial\Omega_{i})}^{2}$
$\displaystyle\leq
c_{2}h^{2(r_{i}-q)-1}_{i}\|v\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.3)
$\displaystyle\sum_{F\in\mathcal{F}}\frac{\delta_{1}}{h_{i}^{3}}\|\llbracket
v-\Pi_{h}v\rrbracket\|_{L_{2}(F)}^{2}$ $\displaystyle\leq
c_{3}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|v\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.4)
$\displaystyle\sum_{F\in\mathcal{F}}\frac{\delta_{0}}{h_{i}}\|\llbracket\nabla_{\Omega}(v-\Pi_{h}v)\rrbracket\|_{L_{2}(F)}^{2}$
$\displaystyle\leq
c_{4}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|v\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.5)
$\displaystyle\sum_{F\in\mathcal{F}}\frac{h_{i}}{\delta_{0}}\|\\{\Delta_{\Omega}(v-\Pi_{h}v)\\}\|_{L_{2}(F)}^{2}$
$\displaystyle\leq
c_{5}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|v\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.6)
$\displaystyle\sum_{F\in\mathcal{F}}\frac{h_{i}^{3}}{\delta_{1}}\|\\{\nabla_{\Omega}\Delta_{\Omega}(v-\Pi_{h}v)\\}\|_{L_{2}(F)}^{2}$
$\displaystyle\leq
c_{6}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|v\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.7)

where $q$ is a positive integer, $r_{i}=\min\\{s_{i},p_{i}+1\\}$ and the
generic constants $c_{2},c_{3},c_{4},c_{5}$ and $c_{6}$ are independent of the
mesh size.

###### Proof

By using the trace inequality (4.13) and the approximation estimate (5.2), we
obtain the proof, see [16]. ∎

In the next lemma, we present estimates in the discrete norms $\|\cdot\|_{h}$
and $\|\cdot\|_{h,*}$ necessary for deriving the a priori error estimates.

###### Lemma 8

Let $v\in V_{D}\cap H^{\textbf{s}}(\Omega,\mathcal{T}_{h})$ for
$\mathbf{s}:=(s_{1},s_{2},\ldots,s_{N})\geq 4$ and $\Pi_{h}v\in V_{h}$ be a
projection. Then, for $\mathbf{p}:=(p_{1},p_{2},\ldots,p_{N})\geq 3,$ we have

$\displaystyle\|v-\Pi_{h}v\|_{h}^{2}$ $\displaystyle\leq
c_{7}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|u\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.8)
$\displaystyle\|v-\Pi_{h}v\|_{h,*}^{2}$ $\displaystyle\leq
c_{8}\sum_{i=1}^{N}h^{2(r_{i}-2)}_{i}\|u\|^{2}_{H^{r_{i}}(\Omega_{i})},$ (5.9)

where $r_{i}:=\min\\{s_{i},p_{i}+1\\}$, $p_{i}$ is the degree of the B-spline
and the constants $c_{7}$ and $c_{8}$ are independent of mesh size $h_{i}.$

###### Proof

The proof follows by using Lemma 7 together with the definitions of the norms
(4.4) and (4.14). ∎

Finally, we present the a priori error estimate in the norm $\|\cdot\|_{h}$
and $\|\cdot\|_{h,*}$ as follows

###### Theorem 5.1

Let $u\in V_{D}\cap H^{\textbf{s}}(\Omega,\mathcal{T}_{h})$ with
$\mathbf{s}=\\{s_{i},i=1,\ldots,N\\},s_{i}\geq 4$ be the solution of (3.6) for
non-negative real numbers $\delta_{0}$ and $\delta_{1}$ chosen as in Theorem
4.1. Let $u_{h}\in V_{h}$ be the solution of (4.2), then there exists $c>0$
independent of $h_{i}$ and $N$ such that the following holds:

$\displaystyle\|u-u_{h}\|_{h}$ $\displaystyle\leq
c\sum_{i=1}^{N}h^{r_{i}-2}_{i}\|u\|_{H^{r_{i}}(\Omega_{i})}$ (5.10)

where $r_{i}=\min\\{s_{i},p_{i}+1\\}.$

###### Proof

Using the coercivity result of Lemma 2, the Galerkin orthogonality (4.3) and
the boundedness of Lemma 4, we can derive the following estimates

$\displaystyle\mu_{c}\|\Pi_{h}u-u_{h}\|_{h}^{2}$ $\displaystyle\leq
a_{h}(\Pi_{h}u-u_{h},\Pi_{h}u-u_{h})=a_{h}(\Pi_{h}u-u,\Pi_{h}u-u_{h})$
$\displaystyle\leq\mu_{b}\|\Pi_{h}u-u\|_{h,*}\|\Pi_{h}u-u_{h}\|_{h}.$

Thus, we have

$\displaystyle\mu_{c}\|\Pi_{h}u-u_{h}\|_{h}$
$\displaystyle\leq\mu_{b}\|\Pi_{h}u-u\|_{h,*}.$ (5.11)

Using (5.11) and the estimates from Lemma 8, we obtain by using triangle
inequality the following estimate

$\displaystyle\|u-u_{h}\|_{h}$
$\displaystyle\leq\|u-\Pi_{h}u\|_{h}+\mu_{b}\|\Pi_{h}u-u\|_{h,*}$
$\displaystyle\leq\left(c_{7}^{1/2}+(\mu_{b}/\mu_{c})c_{8}^{1/2}\right)\sum_{i=1}^{N}h^{r_{i}-2}_{i}\|u\|_{H^{r_{i}}(\Omega_{i})},$
(5.12)

with $c=(c_{7}^{1/2}+(\mu_{b}/\mu_{c})c_{8}^{1/2}).$ ∎

###### Corollary 1

Following the hypothesis and assumptions of Theorem 5.1, then there exists a
constant $c_{9}$ such that the following estimate holds

$\displaystyle\|u-u_{h}\|_{h,*}$ $\displaystyle\leq
c_{9}\sum_{i=1}^{N}h^{r_{i}-2}_{i}\|u\|_{H^{r_{i}}(\Omega_{i})},$ (5.13)

where $r_{i}=\min\\{s_{i},p_{i}+1\\}.$

###### Proof

The proof follows similar argument as in Theorem 5.1 using (4.18) and (4.19)
with
$c_{9}=(c_{7}^{1/2}+(\widetilde{\mu_{b}}/\widetilde{\mu_{c}})c_{8}^{1/2}).$ ∎

The estimates from Theorem 5.1 yields a sharper bound since the estimate from
corollary 1 relies on the norm equivalence i.e. Lemma 5.

##  6 Numerical Results

In this section, we present numerical results for the model problem (1.1). The
numerical experiments are carried out in G+Smo; an open source object-oriented
simulation tool developed solely for IGA, see, [10]. The penalty parameters
$\delta_{0}=\delta_{1}=(p+1)(p+3)/3$ where $p$ is the NURBS degree, see Remark
1. We consider for the open surface, a quarter cylinder and for the closed
surface, a torus as computational domains. Each of the surfaces considered
consists of four $(4)$ patches with matching underlying meshes. The resulting
linear system from the discrete dGIGA scheme (4.2) has been solved using the
SuperLU on an Intel Core (TM) i5-4300 CPU<EMAIL_ADDRESS>The convergence rate is
computed using the formula $rate=\log_{2}\left(e_{i}/e_{i+1}\right),$ where
$e_{i+1}=\|u-u_{h,i+1}\|_{h}$ and $e_{i}=\|u-u_{h,i}\|_{h}$ to study the
discrete solution of the model problem.

###  6.1 Open Surface

We consider a Dirichlet biharmonic problem on an open surface $\Omega$ that is
given by a quarter cylinder in the first quadrant i.e. $x\geq 0$ and $y\geq 0$
with unitary radius and height $L=4.$ The computational domain $\Omega$ is
decomposed into 4 patches, with each of the patches having a unit height of
one $(1)$ and depicted by a different color as seen on the left-hand side of
Figure 4. The knot vectors representing the geometry of each patch are given
by $\Xi_{1}=\\{0,0,0,1,1,1\\}$ and $\Xi_{2}=\\{0,0,1,1\\}$ in the
$\xi_{1}-$direction and $\xi_{2}-$direction respectively. The exact solution
is chosen as $u(\phi,z)=\varrho g_{\phi,1}(\phi)g_{z}(z),$ where
$\phi:=\arctan\left(\frac{x}{y}\right),g_{\phi,1}(\phi)=(1-\cos(\phi))(1-\sin(\phi))$
and $g_{z}(z)=\sin(\sigma\pi z/L)$ are in cylindrical coordinates $(r,\phi,z)$
which yields the source function

$f(\phi,z)=\frac{\varrho
g_{z}(z)\bigg{(}2\pi^{4}\sigma^{4}+(\pi^{2}\sigma^{2}+4L^{2})^{2}\sin(2\phi)-2(\pi^{2}\sigma^{2}+L^{2})^{2}(\sin(\phi)+\cos(\phi))\bigg{)}}{2L^{4}}.$

In the numerical experiments, we set
$\sigma=3,\varrho=1/\left(3/2-\sqrt{2}\right)$ and present the contours of the
solution, see, Figure 4 (right). The rate of convergence for all four $(4)$
dGIGA schemes with respect to the discrete norm $\|\cdot\|_{h}$ is presented
in Table 1 by successive mesh refinement for NURBS degree $2\leq p\leq 6.$ We
observe the optimal convergence rate $\mathcal{O}(h^{(p-1)})$ as theoretically
predicted in Theorem 5.1. Since the discrete norm $\|\cdot\|_{h}$ is indeed a
norm on $H^{2}(\Omega,\mathcal{T}_{h}),$ we require at least a continuously
differentiable basis function i.e. $C^{p-1},p\geq 2$ to obtain optimal
convergence rates.

Figure 4: $2-$dimensional computational domain of a quarter cyliinder (left) consisting of four patches illustrated with different colors and the solution contours (right). Method | SIPG | SSIPG1 | SSIPG2 | NIPG
---|---|---|---|---
| $0.95$ | $0.99$ | $0.94$ | $1.14$
$p=2$ | $1.00$ | $1.00$ | $1.00$ | $1.00$
| $1.00$ | $1.00$ | $1.00$ | $1.00$
| $2.01$ | $2.05$ | $2.02$ | $2.14$
$p=3$ | $2.00$ | $2.00$ | $2.00$ | $2.04$
| $2.00$ | $2.00$ | $2.00$ | $2.00$
| $3.03$ | $3.05$ | $3.15$ | $3.43$
$p=4$ | $3.00$ | $3.00$ | $3.09$ | $3.11$
| $3.00$ | $3.00$ | $3.00$ | $3.00$
| $4.07$ | $4.26$ | $4.04$ | $4.13$
$p=5$ | $4.02$ | $4.02$ | $4.00$ | $4.01$
| $4.00$ | $4.00$ | $4.00$ | $4.00$
| $5.06$ | $5.13$ | $5.01$ | $5.01$
$p=6$ | $5.01$ | $5.03$ | $5.00$ | $5.00$
| $5.00$ | $5.00$ | $5.00$ | $5.00$
Table 1: Convergence rates for the open surfaace problem $\|\cdot\|_{h}.$

###  6.2 Closed Surface

We consider for the closed surface, a torus,

$\Omega=\\{(x,y)\in(-3,3)^{2},z\in(-1,1):\;r^{2}=z^{2}+(\sqrt{x^{2}+y^{2}}-R^{2})\\},$

that is decomposed into 4 patches as depicted on the left-hand side of Figure
5. Each of the NURBS patches is described by the knot vectors
$\Xi_{2}=\\{0,0,0,1,1,1\\}$ and
$\Xi_{1}=\\{0,0,0,0.25,0.25,0.50,0.50,0.75,0.75,1,1,1\\}.$ For the surface
biharmonic problem, we consider an exact solution given by
$u(\phi,\theta)=\sin(3\phi)\cos(3\theta+\phi),$ see also [13]. We chose the
exact solution $u$ and the force term $f$ such that the zero mean
compatibility condition holds. The rate of convergence for all four $(4)$
dGIGA schemes with respect to the discrete norm $\|\cdot\|_{h}$ is presented
in Table 2 by successive mesh refinement for NURBS degree $2\leq p\leq 6.$ We
observe the optimal convergence rate $\mathcal{O}(h^{(p-1)})$ as theoretically
predicted in Theorem 5.1. Since the discrete norm $\|\cdot\|_{h}$ is indeed a
norm on $H^{2}(\Omega,\mathcal{T}_{h}),$ we require at least a continuously
differentiable basis function i.e. $C^{p-1},p\geq 2$ to obtain optimal
convergence rates.

Figure 5: $2-$dimensional computational domain of a torus (left) consisting of four patches illustrated with different colors and the solution contours (right). Method | SIPG | SSIPG1 | SSIPG2 | NIPG
---|---|---|---|---
| $1.09$ | $1.12$ | $1.07$ | $1.11$
$p=2$ | $1.00$ | $1.01$ | $1.00$ | $1.00$
| $1.00$ | $1.00$ | $1.00$ | $1.00$
| $2.02$ | $2.09$ | $2.05$ | $2.11$
$p=3$ | $2.00$ | $2.00$ | $2.00$ | $2.00$
| $2.00$ | $2.00$ | $2.00$ | $2.00$
| $3.01$ | $3.07$ | $3.03$ | $3.10$
$p=4$ | $3.00$ | $3.00$ | $3.00$ | $3.00$
| $3.00$ | $3.00$ | $3.00$ | $3.00$
| $4.02$ | $4.07$ | $4.10$ | $4.09$
$p=5$ | $4.00$ | $4.00$ | $4.00$ | $4.00$
| $4.00$ | $4.00$ | $4.00$ | $4.00$
| $5.02$ | $5.03$ | $5.10$ | $5.07$
$p=6$ | $5.00$ | $5.00$ | $5.00$ | $5.00$
| $5.00$ | $5.00$ | $5.00$ | $5.00$
Table 2: Convergence rates for the closed surfaace problem $\|\cdot\|_{h}.$

## Conclusion

We have presented a priori error estimates for the multipatch discontinuous
Galerkin isogeometric analysis (dGIGA) for the surface biharmonic problem on
orientable computational domains. We assumed non-overlapping subdomains
usually referred to as patches such that the solution could be discontinuous
on the internal facets and applied interior penalty Galerkin techniques. We
arrived at four bilinear forms namely; symmetric (SIPG), non-symmetric (NIPG),
semi-symmetric 1 (SSIPG1) and semi-symmetric 2 (SSIPG2) interior penalty
Galerkin methods. We showed optimal a priori error estimates with respect to
two discrete norms $\|\cdot\|_{h}$ and $\|\cdot\|_{h,*}$ and presented
numerical results for closed and open surfaces that confirmed the analysis
presented. We will extend the current results to biharmonic problems with
singularities as treated for the second order elliptic PDE in [15, Chapter 4].

## References

  * [1] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces.  Pure and Applied Mathematics 140, Elsevier/Academic Press, second edition, 2008.
  * [2] Y. Bazilevs, L. Beirão da Veiga, J. A. Cottrell, T. J. R. Hughes, and G. Sangalli.  Isogeometric analysis: Approximation, stability and error estimates for $h$-refined meshes.  Comput. Methods Appl. Mech. Engrg., 194:4135–4195, 2006.
  * [3] L. Beirão da Veiga, A. Buffa, G. Sangalli, and R. Vázquez.  Mathematical analysis of variational isogeometric methods.  Acta Numerica, 23:157–287, 5 2014.
  * [4] P. G. Ciarlet.  The Finite Element Method for Elliptic Problems.  Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 2002\.
  * [5] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs.  Isogeometric analysis: Toward Integration of CAD and FEA.  John Wiley & Sons, Chichester, 2009.
  * [6] A. Dedner, P. Madhavan, and B. Stinner.  Analysis of the discontinuous Galerkin method for elliptic problems on surfaces.  IMA J. Numer. Anal., 33(3):952–973, 2013.
  * [7] G. Dziuk and C.M. Elliott.  Finite element methods for surface PDEs.  Acta Numerica, 22:289–396, 2013.
  * [8] D. Halpern and E. O. Jensen andJ. B. Grotberg.  A theoretical study of surfactant and liquid delivery into the lung.  J. Appl Physiol, 85(1):333–352, 1998.
  * [9] M. Hofer and H. Pottmann.  Energy-minimizing splines in manifolds.  ACM Trans. Graph., 23(3):284–293, 2004.
  * [10] B. Jüttler, U. Langer, A. Mantzaflaris, S. E. Moore, and W. Zulehner.  Geometry + Simulation Modules: Implementing Isogeometric Analysis.  PAMM, 14(1):961–962, 2014.
  * [11] U. Langer, A. Mantzaflaris, S. E. Moore, and I. Toulopoulos.  Multipatch discontinuous galerkin isogeometric analysis.  In Bert Jüttler and Bernd Simeon, editors, Isogeometric Analysis and Applications 2014, volume 107 of Lecture Notes in Computational Science and Engineering, pages 1–32. Springer, 2015.
  * [12] U. Langer and S. E. Moore.  Domain decomposition methods in science and engineering xxii.  In T. Dickopf, J. M. Gander, L. Halpern, R. Krause, and F. Luca Pavarino, editors, Domain Decomposition Methods in Science and Engineering XXII, chapter Discontinuous Galerkin Isogeometric Analysis of Elliptic PDEs on Surfaces, pages 319–326. Springer, Cham, 2016.
  * [13] K. Larsson and M. G. Larson.  A continuous/discontinuous galerkin method and a priori error estimates for the biharmonic problem on surfaces.  Mathematics of Computation, 8:2613–2649, 2017.
  * [14] Facundo M\̇lx@bibnewblockImplicit brain imaging.  NeuroImage, 23:S179–S188, 2004.  Mathematics in Brain Imaging.
  * [15] S. E. Moore.  Nonstandard Discretization Strategies In Isogeometric Analysis for Partial Differential Equations.  PhD thesis, Johannes Kepler University, January 2017.
  * [16] S. E. Moore.  Discontinuous galerkin isogeometric analysis for the biharmonic equation.  Computers & Mathematics with Applications, 76(4):673 – 685, 2018\.
  * [17] S. E. Moore.  Discontinuous galerkin isogeometric analysis for elliptic problems with discontinuous diffusion coefficients on surfaces.  Numerical Algorithms, 284:1075–1094, 2020.
  * [18] T. G. Myers and J. P. F. Charpin.  A mathematical model for atmospheric ice accretion and water flow on a cold surface.  International Journal of Heat and Mass Transfer, 47(25):5483 – 5500, 2004.
  * [19] O. Nemitz, M.B. Nielsen, M. Rumpf, and R. Whitaker.  Finite element methods on very large, dynamic tubular grid encoded implicit surfaces.  SIAM J. on Sci. Comput., 31(3):2258–2281, 2009.
  * [20] V. P. Nguyen, P. Kerfriden, M. Brino, S. P. A. Bordas, and E. Bonisoli.  Nitsche’s method for two and three dimensional nurbs patch coupling.  Comput. Mech., 53(6):1163–1182, Jun 2014.
  * [21] E. Süli and I. Mozolevski.  hp-version interior penalty DGFEMs for the biharmonic equation.  Comput. Methods Appl. Mech. Engrg., 196(13-16):1851 – 1863, 2007\.
  * [22] A. Tagliabue, L. Dede, and A. Quarteroni.  Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics.  Computers & Fluids, 102(0):277 – 303, 2014.
  * [23] S. Timoshenko and S. Woinowsky-Krieger.  Theory of plates and shells.  Engineering societies monographs. McGraw-Hill, 1959.
  * [24] G. Wheeler.  Fourth order geometric evolution equations.  Bulletin of the Australian Mathematical Society, 82(3):523–524, 2010.