Residual minimization for isogeometric analysis in reduced and mixed forms

TL;DR

This paper introduces a residual minimization framework for isogeometric analysis that uses highly-continuous basis functions for trial spaces and lower-regularity basis functions for test spaces, improving efficiency while maintaining accuracy.

Contribution

It proposes a mixed formulation with residual minimization that allows using less regular test functions in isogeometric analysis, enhancing computational efficiency.

Findings

Variationally-stable formulations verified numerically

Equivalent formulations demonstrated through numerical tests

Framework reduces regularity requirements for test spaces

Abstract

Most variational forms of isogeometric analysis use highly-continuous basis functions for both trial and test spaces. For a partial differential equation with a smooth solution, isogeometric analysis with highly-continuous basis functions for trial space results in excellent discrete approximations of the solution. However, we observe that high continuity for test spaces is not necessary. In this work, we present a framework which uses highly-continuous B-splines for the trial spaces and basis functions with minimal regularity and possibly lower order polynomials for the test spaces. To realize this goal, we adopt the residual minimization methodology. We pose the problem in a mixed formulation, which results in a system governing both the solution and a Riesz representation of the residual. We present various variational formulations which are variationally-stable and verify their equivalence numerically via numerical tests.

Equations

$$-\nabla \cdot (\kappa \nabla u-\mathbf{\beta }u)+\gamma u$$

$$u$$

$$\mathbf{q}-\kappa \nabla u+\mathbf{\beta }u$$

$$-\nabla \cdot \mathbf{q}+\gamma u$$

$$u$$

$$b(w,u)=\ell (w),\forall w\in W,$$

$$b((w,\mathbf{p}),(u,\mathbf{q}))=\ell ((w,\mathbf{p})),\forall (w,\mathbf{p})\in W^u\times W^{\mathbf{q}},$$

$$b(w,u)$$

$$\ell (w)$$

$$b(w,u)$$

$$\ell (w)$$

$$b((w,\mathbf{p}),(u,\mathbf{q}))$$

$$=(w,-\nabla \cdot \mathbf{q}+\gamma u{)}_{\Omega }+(\mathbf{p},\mathbf{q}-\kappa \nabla u+\mathbf{\beta }u{)}_{\Omega },$$

$$\ell ((w,\mathbf{p}))$$

$$b((w,\mathbf{p}),(u,\mathbf{q}))$$

$$=(w,-\nabla \cdot \mathbf{q}+\gamma u{)}_{\Omega }+(\mathbf{p},\mathbf{q}+\mathbf{\beta }u{)}_{\Omega }+(\nabla \cdot (\kappa \mathbf{p}),u{)}_{\Omega },$$

$$\ell ((w,\mathbf{p}))$$

$$b((w,\mathbf{p}),(u,\mathbf{q}))$$

$$=(\nabla w,\mathbf{q}{)}_{\Omega }+(w,\gamma u{)}_{\Omega }+(\mathbf{p},\mathbf{q}-\kappa \nabla u+\mathbf{\beta }u{)}_{\Omega },$$

$$\ell ((w,\mathbf{p}))$$

$$b((w,\mathbf{p}),(u,\mathbf{q}))$$

$$=(\nabla w,\mathbf{q}{)}_{\Omega }+(w,\gamma u{)}_{\Omega }+(\mathbf{p},\mathbf{q}+\mathbf{\beta }u{)}_{\Omega }+(\nabla \cdot (\kappa \mathbf{p}),u{)}_{\Omega },$$

$$\ell ((w,\mathbf{p}))$$

$$b(\mathbf{p},\mathbf{q})$$

$$\ell (\mathbf{p})$$

$${\theta }_0^j(x)$$

$${\theta }_p^j(x)$$

$$V_{p,k}^h$$

$$V_{p,k}^{\mathbf{q},h}$$

$$g(w^h,{\varphi }^h)+b(w^h,u^h)$$

$$b({\varphi }^h,v^h)$$

$$g(v,w)=\sum _{K\in {\mathcal{T}}_h}{\tau }_0(v,w{)}_K+{\tau }_1{h}^{{\iota }_1}(\nabla v,\nabla w{)}_K+{\tau }_2{h}^{{\iota }_2}(\Delta v,\Delta w{)}_K,\forall v,w\in W_{q,l}^h,$$

$$g((w^h,{\mathbf{p}}^h),({\varphi }^h,{\mathbf{\psi }}^h))+b((w^h,{\mathbf{p}}^h),(u^h,{\mathbf{q}}^h))$$

$$b(({\varphi }^h,{\mathbf{\psi }}^h),(v^h,{\mathbf{r}}^h))$$

$$g((v,\mathbf{r}),(w,\mathbf{p}))$$

$$+{\tau }_5(\mathbf{r},\mathbf{p}{)}_K+{\tau }_6{h}^{{\iota }_4}(\nabla \cdot \mathbf{r},\nabla \cdot \mathbf{p}{)}_K,$$

$$\left[\begin{array}{cc}G& B\\ B^T& 0\end{array}\right]\left[\begin{array}{c}\Phi \\ U\end{array}\right]=\left[\begin{array}{c}L\\ 0\end{array}\right],$$

$$B^T{G}^{-1}(L-BU)=B^T{G}^{-1}L-B^T{G}^{-1}BU=0.$$

$$|u-u^h{|}_{1,\Omega }\le C{h}^p,$$

$$\Vert \mathbf{q}-{\mathbf{q}}^h{\Vert }_{0,\Omega }=\Vert (\kappa \nabla u-\mathbf{\beta }u)-{\mathbf{q}}^h{\Vert }_{0,\Omega }.$$

Full text

1 Introduction

Isogeometric analysis, introduced in 2005 [26, 12], is a widely-used numerical method for solving partial differential equations (PDEs). This analysis unifies the finite element methods with computer-aided design tools. Within the framework of the classic Galerkin finite element methods, isogeometric analysis uses as basis functions the functions that describe the geometry in the computer-aided design (CAD) and engineering (CAE) technologies, namely, B-splines or non-uniform rational B-splines (NURBS) and their generalizations. These CAD/CAE basis functions may possess higher continuity. This smoothness improves the numerical approximations of PDEs which have highly regular solutions. Bazilevs et al. established the approximation, stability, and error estimates in [2]. Cottrell et al. in [13] first used isogeometric analysis to study the structural vibrations and wave propagation problems. Spectral analysis shows that isogeometric elements significantly improve the accuracy of the spectral approximation when compared with classical finite elements. In [27], the authors explored the additional advantages of isogeometric analysis on the spectral approximation over finite elements. Moreover, the work [22, 35, 6, 21] minimized the spectral approximation errors for isogeometric elements. The minimized spectral errors possess superconvergence (two extra orders) as the mesh is refined.

Galerkin isogeometric analysis uses highly-continuous B-splines or NURBS for both trial and test spaces. For PDEs with smooth solutions (high regularities), isogeometric elements based on smooth functions render solutions which have a better physical interpretation. For example, for a PDE with a solution in $H^{2}$, both the exact solution field and exact flux field are expected to be globally continuous. Isogeometric analysis with $C^{1}$ B-spline basis functions produces a globally continuous flux field while the classical finite element method does not. Therefore, for smooth problems, it makes sense to use highly-continuous B-splines for the trial spaces. However, it is not necessary to use highly-continuous B-splines for the test spaces. In the view of the work [11] which explores the cost of continuity, we expect an extra cost for solving the resulting system when we apply unnecessary high continuities for the basis functions in the test spaces. This extra cost per degree of freedom for highly-continuous discretizations was verified for direct [11, 33, 9] and iterative solvers [10]. Additionally, the work [25, 24] shows that a reduction in the continuity of the trial and test spaces may lead to faster and more accurate solutions. Thus, we seek to develop an isogeometric analysis which uses minimal regularity for the test spaces. In this work, we establish a framework which uses highly-continuous B-splines for the trial spaces and basis functions with minimal regularity for the test spaces.

To realize this goal, we adopt the residual minimization methodology sharing with the Discontinous Petrov-Galerkin (DPG) method, the idea of a Riesz representation for the residual with respect to a stabilizing norm (c.f., [17, 18, 20, 38]). The main idea of DPG is the use of discontinuous or broken basis functions for the test spaces within the Petrov-Galerkin framework. The method computes the optimal test function by using a trial-to-test operator. The goal is to automatically stabilize the discrete formulation [31, 5, 30, 32] without parameter tuning. The DPG method can be interpreted as a minimum residual method in which the residual is measured in the dual test norm [19]. By introducing an auxiliary unknown $\phi$ representing the Riesz representation of residual, the method is cast into a mixed problem. Effectively, we extend to highly continuous isogeometric discretizations the method described in [7] where standard $C^{0}$ finite element basis functions build the trial space while the test space uses broken polynomials.

The proposed method has a feature which distinguishes us from DPG. We use highly-continuous B-splines for space $V^{h}$ in which we seek for the solution, while DPG uses discontinuous basis functions. Consequently, DPG introduces additional unknowns that live on the mesh skeleton; see equation (10) or (39) in [19]. We do not introduce any additional unknowns. Instead, we use basis functions ($C^{-1}$ functions in the discrete $L^{2}$ space or at most $C^{0}$ in the discrete $H^{1}$ space defined on the mesh partition) with minimal regularity in the sense that no inner products are introduced on the mesh skeleton (thus, the bilinear form only involves element integrations) while maintaining the inf-sup condition for the resulting variational formulation.

The rest of this paper is organized as follows. Section 2 describes the problem under consideration and introduces various variational formulations following the DPG framework closely. Section 3 presents isogeometric formulations at the discrete level, while Section 4 shows numerical examples to demonstrate the performance of the formulations. Concluding remarks are given in Section 5.

2 Problem statement

Let $\Omega\subset\mathbb{R}^{d},d=1,2,3,$ be an open bounded domain with Lipschitz boundary $\partial\Omega$. We consider the advection-diffusion-reaction equation: Find $u$ such that

[TABLE]

where $\kappa$ is the diffusion coefficient, ${\boldsymbol{\beta}}$ is the advective velocity vector, $\gamma\geq 0$ is the reaction coefficient, $u$ is the function to be found, and $f$ is a forcing function. This problem can be cast into an equivalent system of two first-order equations

[TABLE]

where ${\boldsymbol{q}}$ is an auxiliary variable standing for the flux.

2.1 General setting

We present the general setting for equations (1) and (2) as they have distinguishing features. Let $V$ and $W$ denote the trial and test space, resepectively. For (2) we let $(u,{\boldsymbol{q}})$ be the solution pair in its trial space $V^{u}\times V^{{\boldsymbol{q}}}$ and let its corresponding test space be $W^{u}\times W^{{\boldsymbol{q}}}$. We specify these spaces in each of the variational formulations. Let $L^{2}(\Omega)=\{v:\int_{\Omega}v^{2}\ \text{d}x<\infty\}$ and denote the induced $L^{2}$ norm as $\|\cdot\|_{0,\Omega}$. Let $C^{k}(\Omega),k=-1,0,1,\cdots,\infty,$ denote the space of functions having continuous (partial) derivatives up to the $k$-th order over the whole domain $\Omega$. In particular, for $k=-1$, we mean that the function in $C^{-1}$ is discontinuous somewhere in $\Omega$. Let $H^{1}(\Omega)=\{v:v\in L^{2}(\Omega),\int_{\Omega}(\nabla v)^{2}\ \text{d}x<\infty\}$ and $H^{1}_{0}(\Omega)$ be the space of all functions in $H^{1}(\Omega)$ vanishing at the boundary $\partial\Omega$. The $H^{1}$ seminorm is defined as $|w|_{1,\Omega}=\|\nabla w\|_{0,\Omega}$ and $H^{1}$ norm is defined as $\|w\|^{2}_{1,\Omega}=\|\nabla w\|^{2}_{0,\Omega}+|w|^{2}_{1,\Omega}$. Finally, $H(\text{div},\Omega)$ consists of all square integrable vector-valued fields on $\Omega$ whose divergence is a function that is also square integrable. Similarly, the $H(\text{div},\Omega)$ norm is defined as $\|{\boldsymbol{p}}\|^{2}_{\text{div},\Omega}=\|{\boldsymbol{p}}\|^{2}_{0,\Omega}+\|\nabla\cdot{\boldsymbol{p}}\|^{2}_{0,\Omega}$.

Let $\mathcal{T}_{h}$ be a partition of $\Omega$ into non-overlapping mesh elements. For simplicity and the purpose of using B-splines in multiple dimensions, we assume the tensor-product structure. Let $K\in\mathcal{T}_{h}$ be a generic element and denote by $\partial K$ its boundary. Let ${\boldsymbol{n}}$ be the outward unit normal vector. Let $(\cdot,\cdot)_{S}$ denote the $L^{2}(S)$ the inner product where $S$ is a $d$ or $d-1$ dimensional domain ($S$ is typically $\Omega,K,\partial\Omega,\partial K$).

At discrete level, we define the finite spaces, namely the finite subspaces of $V,W,V^{u}\times V^{{\boldsymbol{q}}}$, and $W^{u}\times W^{{\boldsymbol{q}}}$. Since isogeometric analysis adopts highly-continuous basis functions, we specify the corresponding finite spaces with both the polynomial order $p$ and the order $k$ of global continuity (that is, $C^{k}$). Let us denote by $V^{h}_{p,k},W^{h}_{p,k}$ and $V^{u,h}_{p,k},V^{{\boldsymbol{q}},h}_{p,k},W^{u,h}_{p,k},W^{{\boldsymbol{q}},h}_{p,k}$ the finite-dimensional subspaces of $V,W$ for (1) and $V^{u},V^{{\boldsymbol{q}}},W^{u},W^{{\boldsymbol{q}}}$ for (2), respectively. Consequently, $V^{h}_{p,k}\subset V,W^{h}_{p,k}\subset W,V^{u,h}_{p,k}\subset V^{u},V^{{\boldsymbol{q}},h}_{p,k}\subset V^{{\boldsymbol{q}}},W^{u,h}_{p,k}\subset W^{u}$, and $W^{{\boldsymbol{q}},h}_{p,k}\subset W^{{\boldsymbol{q}}}$. In this work, we construct the test spaces $W^{h}_{q,l}$ and $W^{u,h}_{q,l}$ using basis functions with lower order polynomials as well as lower regularity.

2.2 Various variational formulations at continuous level

In this section, we present seven variational formulations for both (1) and (2) following closely the six formulations described within the DPG framework in [19]. In particular, we add one more formulation (to the six formulations in DPG framework) which resembles the isogeometric collocation method.

We start with the abstract variational formulations of (1) and (2). The variational weak formulations of (1) can be written as: Find $u\in V$ such that

[TABLE]

where $b(\cdot,\cdot)$ and $\ell(\cdot)$ are bilinear and linear forms, respectively. Similarly, the variational weak formulations of (2) can be written: Find $(u,{\boldsymbol{q}})\in V^{u}\times V^{{\boldsymbol{q}}}$ such that

[TABLE]

where $b(\cdot,\cdot)$ and $\ell(\cdot)$ are bilinear and linear forms defined over two fields. Herein, we refer to $(u,{\boldsymbol{q}})$ as the solution (trial) pair while to $(w,{\boldsymbol{p}})$ as the weighting (test) function pair. Equations (1) and (3) are the primal forms of the PDE and the variational formulation while equations (2) and (4) are their mixed forms. We keep this structure for the forms at discrete level in Section 3. All these forms and their associated spaces are to be specified in each particular formulation as follows.

Primal trivial formulation: Let $V=H^{1}_{0}(\Omega)\cap C^{1}(\Omega),W=L^{2}(\Omega)$. Find $u\in V$ satisfying (3) with

[TABLE] 2. 2.

Primal classical (FEM) formulation: Let $V=H^{1}_{0}(\Omega),W=H^{1}_{0}(\Omega)$. Find $u\in V$ satisfying (3) with

[TABLE] 3. 3.

Mixed trivial formulation: Let $V^{u}=H^{1}_{0}(\Omega),V^{{\boldsymbol{q}}}=H(\text{div},\Omega),W^{u}=L^{2}(\Omega),$ $W^{{\boldsymbol{q}}}=(L^{2}(\Omega))^{d}$. Find $(u,{\boldsymbol{q}})\in V^{u}\times V^{{\boldsymbol{q}}}$ satisfying (4) with

[TABLE] 4. 4.

Mixed classical formulation: Let $V^{u}=L^{2}(\Omega),V^{{\boldsymbol{q}}}=H(\text{div},\Omega),W^{u}=L^{2}(\Omega),$ $W^{{\boldsymbol{q}}}=H(\text{div},\Omega)$. Find $(u,{\boldsymbol{q}})\in V^{u}\times V^{{\boldsymbol{q}}}$ satisfying (4) with

[TABLE] 5. 5.

Mixed classical formulation: Let $V^{u}=H^{1}_{0}(\Omega),V^{{\boldsymbol{q}}}=(L^{2}(\Omega))^{d},W^{u}=H^{1}_{0}(\Omega),$ $W^{{\boldsymbol{q}}}=(L^{2}(\Omega))^{d}$. Find $(u,{\boldsymbol{q}})\in V^{u}\times V^{{\boldsymbol{q}}}$ satisfying (4) with

[TABLE] 6. 6.

Mixed ultraweak formulation: Let $V^{u}=L^{2}(\Omega),V^{{\boldsymbol{q}}}=(L^{2}(\Omega))^{d},W^{u}=H^{1}_{0}(\Omega),W^{{\boldsymbol{q}}}=H(\text{div},\Omega)$. Find $(u,{\boldsymbol{q}})\in V^{u}\times V^{{\boldsymbol{q}}}$ satisfying (4) with

[TABLE] 7. 7.

Reduced flux formulation: Let $V=H(\text{div},\Omega),W=H(\text{div},\Omega)$. Find ${\boldsymbol{q}}\in V$ satisfying (3) with

[TABLE]

Herein, the primal trivial formulation reduces to the isogeometric collocation method when applying constant test functions.

3 Various isogeometric formulations

In this section, we present the isogeometric formulations at the discrete level for both (3) and (4). We first specify the basis functions for all the finite element spaces associated with the mesh configuration $\mathcal{T}_{h}$. For this purpose, we use the Cox-de Boor recursion formula [15, 34] on each dimension and then take tensor-product to obtain the necessary basis functions for multiple dimensions. The Cox-de Boor recursion formula generates $C^{k}$ and $p$-th order B-spline basis functions, where $p=0,1,2,\cdots,$ and $k=-1,0,1,\cdots,p-1$. $C^{-1}$ basis functions generate a finite-dimensional subspace of the $L^{2}(\Omega)$ and $(L^{2}(\Omega))^{d}$, while $C^{k},k=0,1,\cdots,p-1$ basis functions generate a finite-dimensional subspace of the $H^{1}(\Omega)$ and $H(\text{div},\Omega)$. These finite-dimensional subspaces consist of piecewise polynomials of order $p=0,1,2,\cdots$ and of continuity order $k=-1,0,1,\cdots,p-1$.

The definition of the B-spline basis functions in one dimension is as follows. Let $X=\{x_{0},x_{1},\cdots,x_{m}\}$ be a knot vector with knots $x_{j}$, that is, a nondecreasing sequence of real numbers which are called knots. The $j$-th B-spline basis function of degree $p$, denoted as $\theta^{j}_{p}(x)$, is defined as [15, 34]

[TABLE]

We then construct the finite-dimensional subspaces using the B-splines on uniform tensor-product meshes with non-repeating and repeating knots. For a $p$-th order B-spline, a repetition of $k=0,1,\cdots,p$ times of an internal node results in a function of $C^{p-1-k}$ continuity; see [34, 12]. These B-spline basis functions characterize the finite-dimensional subspaces. For example,

[TABLE]

where $p$ and $k$ specify the approximation order and continuity order in each dimension (they can be different in general), respectively. $N_{x},N_{y},N_{z}$ is the total number of basis functions in each dimension. Note that the space $V^{{\boldsymbol{q}},h}_{p,k}$ collapse to standard Raviart-Thomas mixed finite elements [36]; see [23, Section 5] or [3, Section 3] for more details on the construction of these spaces using B-splines.

We now adopt the mixed formulation for (3) at discrete level: Find $(u^{h},\phi^{h})\in V^{h}_{p,k}\times W^{h}_{q,l}$ such that

[TABLE]

where the spaces are chosen such that $\text{dim}(W^{h}_{q,l})\geq\text{dim}(V^{h}_{p,k})$. Herein, $\phi^{h}$ is the negative of the Riesz representation of the residual. The auxiliary bilinear form $g(\cdot,\cdot)$ is an inner product and it produces a Gramm matrix for the purpose of residual minimization. We define it generally as follows

[TABLE]

where $\tau_{i},\iota_{j}\in\mathbb{R},i=0,1,2,j=1,2$ are free parameters. The default setting is $\tau_{0}=1,\tau_{1}=1,\tau_{2}=0$, $\iota_{1}=2,\iota_{2}=0$. Once one specifies an inner product $g(\cdot,\cdot)$, then under inf-sup assumption on the discrete bilinear formulation in (14), the approximate solution $u^{h}$ has a minimal error in the energy norm induced from $g(\cdot,\cdot)$; see, for example [19, 18].

Similarly, we present the mixed formulation for (4) at discrete level: Find $(u^{h},{\boldsymbol{q}}^{h})\in V^{u,h}_{p_{1},k_{1}}\times V^{{\boldsymbol{q}},h}_{p_{2},k_{2}}$ and $(\phi^{h},{\boldsymbol{\psi}}^{h})\in\times W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$ such that

[TABLE]

where the spaces are chosen such that $\text{dim}(W^{h}_{q_{1},l_{1}})\geq\text{dim}(V^{u,h}_{p_{1},k_{1}})$ and $\text{dim}(W^{{\boldsymbol{q}},h}_{q_{2},l_{2}})\geq\text{dim}(V^{{\boldsymbol{q}},h}_{p_{2},k_{2}})$. Herein, the continuity orders corresponding to solution $u$ and flux ${\boldsymbol{q}}$ can be different. Potentially, their polynomial orders can also be different. Similarly, $g((\cdot,\cdot),(\cdot,\cdot))$ is an inner product for the purpose of residual minimization. We define the Gramm product generally as follows, for $(v,{\boldsymbol{r}}),(w,{\boldsymbol{p}})\in W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$,

[TABLE]

where $\tau_{i},\iota_{j}\in\mathbb{R},i=3,4,\cdots,6,j=3,4$ are free parameters. The default setting is $\tau_{3}=1,\tau_{4}=1,\tau_{5}=1,\tau_{6}=1,\iota_{3}=2,\iota_{4}=2$. The formulation (16) is the same as (14) if we view $(u^{h},{\boldsymbol{q}}^{h})$ as a solution pair.

Within these formulations, once we specify all free parameters, the bilinear and linear forms, and space settings, we have a different method. We present the following discrete variational formulations

Let $b(\cdot,\cdot)=b_{1}(\cdot,\cdot),\ell(\cdot)=\ell_{1}(\cdot)$ defined in (5). Let $V^{h}_{p,k}$ consist of B-spline basis functions of continuity $C^{k},k\geq 1$ and $W^{h}_{q,l}$ consist of discontinuous basis functions. The discrete primal trivial formulation is: Find $(u^{h},\phi^{h})\in V^{h}_{p,k}\times W^{h}_{q,l}$ satisfying (14). 2. 2.

Let $b(\cdot,\cdot)=b_{2}(\cdot,\cdot),\ell(\cdot)=\ell_{2}(\cdot)$ defined in (6). Let $V^{h}_{p,k}$ and $W^{h}_{q,l}$ consist of B-spline basis functions of continuity at least $C^{0}$. The discrete primal classical formulation is: Find $(u^{h},\phi^{h})\in V^{h}_{p,k}\times W^{h}_{q,l}$ satisfying (14). 3. 3.

Let $b(\cdot,\cdot)=b_{3}(\cdot,\cdot),\ell(\cdot)=\ell_{3}(\cdot)$ defined in (7). Let the solution space consist of B-spline basis functions of continuity at least $C^{0}$ while the test space consist of discontinuous basis functions. The discrete mixed trivial formulation is: Find $(u^{h},{\boldsymbol{q}}^{h})\in V^{u,h}_{p_{1},k_{1}}\times V^{{\boldsymbol{q}},h}_{p_{2},k_{2}}$ and $(\phi^{h},{\boldsymbol{\psi}}^{h})\in\times W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$ satisfying (16). 4. 4.

Let $b(\cdot,\cdot)=b_{4}(\cdot,\cdot),\ell(\cdot)=\ell_{4}(\cdot)$ defined in (8). Let the solution space and test space for flux consist of B-spline basis functions of continuity at least $C^{0}$ while the test space for $u$ consist of discontinuous basis functions. The discrete mixed classical formulation I is: Find $(u^{h},{\boldsymbol{q}}^{h})\in V^{u,h}_{p_{1},k_{1}}\times V^{{\boldsymbol{q}},h}_{p_{2},k_{2}}$ and $(\phi^{h},{\boldsymbol{\psi}}^{h})\in\times W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$ satisfying (16). 5. 5.

Let $b(\cdot,\cdot)=b_{5}(\cdot,\cdot),\ell(\cdot)=\ell_{5}(\cdot)$ defined in (9). Let the solution space and test space for $u$ consist of B-spline basis functions of continuity at least $C^{0}$ while the test space for flux consist of discontinuous basis functions. The discrete mixed classical formulation II is: Find $(u^{h},{\boldsymbol{q}}^{h})\in V^{u,h}_{p_{1},k_{1}}\times V^{{\boldsymbol{q}},h}_{p_{2},k_{2}}$ and $(\phi^{h},{\boldsymbol{\psi}}^{h})\in\times W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$ satisfying (16). 6. 6.

Let $b(\cdot,\cdot)=b_{6}(\cdot,\cdot),\ell(\cdot)=\ell_{6}(\cdot)$ defined in (10). Let the solution and test spaces consist of B-spline basis functions of continuity at least $C^{0}$. The discrete mixed ultraweak formulation II is: Find $(u^{h},{\boldsymbol{q}}^{h})\in V^{u,h}_{p_{1},k_{1}}\times V^{{\boldsymbol{q}},h}_{p_{2},k_{2}}$ and $(\phi^{h},{\boldsymbol{\psi}}^{h})\in\times W^{w,h}_{q_{1},l_{1}}\times W^{{\boldsymbol{q}},h}_{q_{2},l_{2}}$ satisfying (16).

The primal classical formulation can be reduced to the standard finite element and isogeometric element methods. If we set $k=0,l=0$ for $p\geq 1$, then we have $\phi^{h}=0$ due to the orthogonal condition (second equation) in (14). Thus, the method reduces to finite element method. Similarly, if we set $k=l=1,\cdots,p-1$ for $p\geq 2$, then we have $\phi^{h}=0$ in (14) as well and the method reduces to isogeometric analysis. For all other discrete variational formulations, we may constrain the solution and test spaces in such a way that the variational forms we discuss render standard discretization techniques when the Riesz representation of the residual is identically zero.

For the scenarios where we use different trial and test spaces, we obtain a non-zero discrete representation of the residual $\Phi$, which we use as an error estimator to guide the refinements of the meshes accordingly. The error estimators are defined in the sense of $G$ which is a result of the bilinear form $g(\cdot,\cdot)$. We refer to the DPG work [17, 18, 20, 38] for more details in this direction.

These discrete variational formulations lead to a linear matrix problem

[TABLE]

where $L$ is the corresponding forcing term arising from the linear form $\ell(\cdot)$, $\Phi,U$ are the solution pairs for the mixed formulations, $G$ represents the Gramm product in which we minimize the residual that arises from the bilinear form $g(\cdot,\cdot)$ and $B$ is the matrix arising from the bilinear form $b(\cdot,\cdot)$. We solve the first equation in (18) and substitute to the second equation in (18) to obtain

[TABLE]

Herein, $BU-L$ is the residual and (19) is a least-square type of problem.

For all the discrete variational formulations, we verify the following optimal error convergence rates numerically in the Section 4:

[TABLE]

where $C$ is a constant independent of $h$. For the approximate fluxes, we define the following errors in $L^{2}$ norm

[TABLE]

The optimal convergence rate is

[TABLE]

where $C$ is a constant independent of $h$.

4 Numerical experiments

The main result of these numerical tests is that the various formulations we discuss above are equivalent in the sense of resulting the same (optimal) error convergence rates. We focus on 2D and consider the problem (1) with $\kappa=1,\gamma=1,{\boldsymbol{\beta}}=(1,1)^{T}$ and a manufactured solution

[TABLE]

$f$ is the corresponding forcing satisfying (1). The true flux ${\boldsymbol{q}}$ is calculated from (2). We apply all six variational formulations, namely, two primal and four mixed formulations, to solve this problem.

Both the primal and mixed trivial formulations do not involve any integration by parts in their variational formulations and their test functions are in $L^{2}$. The difference is that the primal formulation results from (1) while the mixed formulation results from (2). We use $C^{-1}$ basis functions for the test spaces for these formulations. Consequently, we obtain a matrix $G$ in (18) which is a block-diagonal matrix. The independence of each elemental block from the rest allows these methods to be computable efficiently (cheap elemental inversion) and thus relevant for practical purposes. All other formulations involve integration by parts to pass derivatives to the test functions, which in return results in a matrix $G$ in (18) which is not block-diagonal. This makes these formulations interesting from the theoretical point of view, but untractable in most general meshes, even though splitting schemes make these methods viable on tensor-product meshes [29]. Thus, once we show that all these formulations deliver equivalent results at the discrete level, we focus on these strong/trivial formulations as they are computationally advantageous. We chose not to compare against well established DPG technologies to brake test spaces, as the goal is simply to show the equivalence of the different variational forms rather than derive alternative computable methods. Lastly, to compare the primal trivial formulation with the first-order system least-square (FOSLS) method [4], the difference is that the primal trivial formulation does not lead to first order system and it introduces a Gramm matrix to solve for the residual errors simultaneously.

Figure 1 shows the $H^{1}$ semi-norm errors when using all six formulations for $p=2,3,4,5$. For the mixed formulations, we approximate the flux ${\boldsymbol{q}}^{h}$ using basis functions of the same polynomial order as for the solution $u^{h}$. The parameters of the bilinear form $g(\cdot,\cdot)$ are set to be the default values. The mesh configurations are $5\times 5,10\times 10,20\times 20,40\times 40$. Herein, we plot the errors in natural logarithmic scale. As predicted, in all scenarios, the $H^{1}$ semi-norm errors converge in optimal rates, that is, order $p$ for all the variational formulations. This confirms the theoretical result (20). Therefore, all formulations are equivalent in the quality of the approximation they deliver. Interestingly, all primal trivial forms for even orders have optimal convergence rates, but the constants seem to be worse for this discretization for even orders than all the other weak forms we compared against. Nevertheless, for odd order polynomials both the rate of convergence and constant are comparable to those observed for all the other variational forms.

The mixed formulations introduce auxiliary variables for the fluxes. Thus, for the same mesh configuration, the resulting matrix systems of the mixed formulations have dimensions which are three times larger than the primal ones. However, the mixed formulations deliver more accurate fluxes approximations.

Figure 2 shows the flux errors in the $L^{2}$ norm when using the mixed forms for $p=2,3,4,5.$ The fluxes are of optimal convergence rates $p+1$, which confirms the error estimate (22). For $p=2,4$, we observe super-convergence rates approximately $p+2$. Once again, these results numerically verify that these four formulations are equivalent. Therefore, for the following numerical results, we only show the results for the strong primal and mixed formulations.

We also study the formulations with different choices of the auxiliary bilinear forms $g(\cdot,\cdot)$ in (15) and (17). Figure 3 shows the numerical errors $u^{h}$ in $H^{1}$ semi-norm when using the inner product $g(v,w)=(v,w)_{\Omega}$ in (15) for the trivial primal form and $g((v,{\boldsymbol{r}}),(w,{\boldsymbol{p}}))=(v,w)_{\Omega}+({\boldsymbol{r}},{\boldsymbol{p}})_{\Omega}$ in (17) for the trival mixed form. Other inner products are also possible, which is in agreement with the DPG methodology. The convergence rates in the $H^{1}$ semi-norm of all these scenarios are optimal.

5 Concluding remarks

We introduce a residual minimization based mixed formulation for solving partial differential equations. A key feature of this method is the framework which uses highly-continuous B-splines for the trial spaces and basis functions with minimal regularity and lower order polynomials for the test spaces. The method shares with the discontinuous Petrov-Galerkin methodology the idea of stabilizing the formulation considering an adequate norm for the test space and unifies the interpretation of several methods such as the classical finite element method, isogeometric analysis, and isogeometric collocation methods. Under the standard assumption, the proposed variational formulations are stable and result in optimal approximation properties.

Acknowledgement

This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth Scientific Industrial Research Organisation, CSIRO, of Australia. The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Sciences (ICES) of the University of Texas at Austin has partially supported the visits of VMC to ICES.