k-Version of Finite Element Method for BVPs and IVPs

Abstract

The paper presents k-version of the finite element method for boundary value problems (BVPs) and initial value problems (IVPs) in which global differentiability of approximations is always the result of the union of local approximations. The higher order global differentiability approximations (HGDA/DG) are always p-version hierarchical that permit use of any desired p-level without effecting global differentiability. HGDA/DG are true $C^i$, $C^{ij}$, $C^{ijk}$, hence the dofs at the nonhierarchical nodes of the elements are transformable between natural and physical coordinate spaces using calculus. This is not the case with tensor product higher order continuity elements discussed in this paper, thus confirming that the tensor product approximations are not true $C^i$, $C^{ijk}$, $C^{ijk}$ approximations. It is shown that isogeometric analysis for a domain with more than one patch can only yield solutions of class $C^0$. This method has no concept of finite elements and local approximations, just patches. It is shown that compariso of this method with k-version of the finite element method is meaningless. Model problem studies in ${\mathbb{R}}^2$ establish accuracy and superior convergence characteristics of true $C^{ij}$ p-version hierarchical local approximations presented in this paper over tensor product approximations. Convergence characteristics of p-convergence, k-convergence and $pk$-convergence are illustrated for self adjoint, non-self adjoint and non-linear differential operators in BVPs. It is demonstrated that h, p and k are three independent parameters in all finite element computations. Tensor product local approximations and other published works on k-version and their limitations are discussed in the paper and are compared with present work.

1. Introduction

In physical sciences, the mathematical descriptions of the deformation of continuous media (solid or fluent continua) derived using conservation and balance laws of thermodynamics and associated constitutive theories lead to initial value problems (IVPs) or boundary values problems (BVPs). Initial value problems describe evolution, hence the dependent variables in their mathematical description exhibit simultaneous dependence on spatial coordinates and time. In case of boundary value problems that are stationary states of evolutions described by IVPs, the dependent variables only exhibit dependence on spatial coordinates. In the absence of permanent damage to the continuum during deformation, the solutions for the BVPs and IVPs related to deformation of the continua must be continuous and must possess continuous derivatives up to certain order which depends upon the differential operator describing the IVP or the BVP and their theoretical solutions. Physical processes do not admit being discontinuous unless there is damage to the continuum. In many cases what is viewed as discontinuous behavior is often an issue of scale. For example, the shocks in compressible flow, if viewed at a bigger scale, may appear as a discontinuous phenomenon, but on closer examination at a finer scale these are indeed continuous and differentiable. This basic assumption that solution of all BVPs and IVPs in physical sciences are generally analytic is the foundation of the work presented in this paper.

1.1. Boundary Value Problems (BVPs)

Consider the following BVP

$$\begin{array}{c}\hfill A\varphi \left(x\right)-f\left(x\right)=0\forall x\in {\mathsf{\Omega }}_x\end{array}$$

(1)

in which a A is differentiable operator, the variable $\varphi$ is unknown and depends on f, and x is an independent variable that could be x or x, y or x, y, z, but for the sake of simplicity we continue with x only and ${\mathsf{\Omega }}_x$ is the domain of definition. Let the differential operator contain derivatives of $\varphi$ up to orders $2m$. A theoretical solution $\varphi$ of (1) must be at least of class $C^{2m}\left({\overline{\mathsf{\Omega }}}_x\right)$ if derivatives of $\varphi$ of up to orders $2m$ are continuous and differentiable. This is due to the fact that the continuity and differentiability aspects of the physics used to derive BVP (1) necessitate this feature of the theoretical solution of (1). A theoretical solution $\varphi$ of (1) can be of class higher than $C^{2m}\left({\overline{\mathsf{\Omega }}}_x\right)$, this of course depends upon the non-homogeneous function $f\left(x\right)$. Thus, a solution of (1) must belong to subspace V such that

$$\begin{array}{c}\hfill \varphi \in V\subset H^k\left({\overline{\mathsf{\Omega }}}_x\right);k\ge 2m+1\end{array}$$

(2)

V is a subspace of scalar product space $H^k\left({\overline{\mathsf{\Omega }}}_x\right)$. k is the order of the scalar product space. A scalar product space of order k contains functions that have continuous derivatives up to orders $k-1$ and the derivative of order k can have isolated discontinuities for some $x\in {\overline{\mathsf{\Omega }}}_x$, but are square integrable. For the BVP (1), k must at least be $2m+1$ (functions $\varphi$ of class $C^{2m}$), but $k>2m+1$ or $k\to \infty$ (if for example $\varphi$ consists of trigonometric functions) are admissible as well. When $k=2m+1$, V is called minimally conforming space, i.e., a solution of (1) cannot be of a class lower than $C^{2m}\left({\overline{\mathsf{\Omega }}}_x\right)$.

When seeking finite element solutions for (1) we discretize the spatial domain ${\overline{\mathsf{\Omega }}}_x$ into ${\overline{\mathsf{\Omega }}}_x^T={\displaystyle \bigcup _e}{\overline{\mathsf{\Omega }}}_x^e$ in which ${\overline{\mathsf{\Omega }}}_x^e$ is a finite element. ${\overline{\mathsf{\Omega }}}_x^T$ is called discretization of ${\overline{\mathsf{\Omega }}}_x^{}$ or a finite element mesh. We consider ${\varphi }_h^e$ to be local approximation of $\varphi$ over ${\overline{\mathsf{\Omega }}}_x^e$, the domain of an element e and ${\varphi }_h$ to be the approximation of $\varphi$ over ${\overline{\mathsf{\Omega }}}_x^T$, global approximation of $\varphi$ for the discretization ${\overline{\mathsf{\Omega }}}_x^T$ such that

$$\begin{array}{cc}\hfill {\varphi }_h\left(x\right)& =\bigcup ^{}_e{\varphi }_h^e\left(x\right)\forall x\in {\overline{\mathsf{\Omega }}}_x^T\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \mathrm{and}{\varphi }_h^e\left(x\right)& =\sum _i^n{N}_i\left(x\right){\delta }_i^e\forall x\in {\overline{\mathsf{\Omega }}}_x^e\hfill \end{array}$$

(4)

In (4), $N_i\left(x\right)$ are approximation functions and ${\delta }_i^e$ are nodal dofs. The approximation functions $N_i\left(x\right)$ are generally obtained using linear combination of algebraic monomials (interpolation theory) up to degree p. This is convenient as well as prudent. Thus, ${\varphi }_h^e$ of class $C^p({\overline{\mathsf{\Omega }}}_x^e)$ is always ensured. However, the differentiability ${\varphi }_h\left(x\right)\forall x\in {\overline{\mathsf{\Omega }}}_x^T$ (global differentiability) depends upon (3), i.e., the union of local approximations (${\varphi }_h^e$’s) controls the global differentiability of ${\varphi }_h$ over the discretization ${\overline{\mathsf{\Omega }}}_x^T$. It is obvious that ${\varphi }_h$ over each ${\mathsf{\Omega }}_x^e$ is of class $C^p\left({\mathsf{\Omega }}_x^e\right)$, hence the differentiability of ${\varphi }_h\forall x\in {\overline{\mathsf{\Omega }}}_x^T$ is controlled by the differentiability of ${\varphi }_h$ over the inter-element boundaries of elements ${\overline{\mathsf{\Omega }}}_x^e$. Thus, we could design local approximations ${\varphi }_h^e$ (in (4)) of some polynomial degree p such that (3) gives us the desired global differentiability of ${\varphi }_h$ at the inter-element boundaries, always less than the differentiability of order p of ${\varphi }_h$ over ${\mathsf{\Omega }}_x^e$, interior of each element. The local approximation ${\varphi }_h^e$ over ${\overline{\mathsf{\Omega }}}_x^e$ of polynomial order p that yield global differentiability of order q at the inter-element boundaries ($q<p$) are called local approximations of class $C^q({\overline{\mathsf{\Omega }}}_x^e)$. Thus, the global differentiability of ${\varphi }_h$ is of order q due to the local approximation ${\varphi }_h^e$ yielding global differentiability of order q at the inter-element boundaries. Thus, to achieve global differentiability of ${\varphi }_h$ of certain order we must design local approximations such that (3) in fact gives us the desired global differentiability of ${\varphi }_h$ at the inter-element boundaries.

1.2. k-Version of Finite Element Method

We note that (Surana et al. [1,2,3]) order of the space k in $H^{k,p}\left({\overline{\mathsf{\Omega }}}_x^{}\right)$ giving global differentiability of order $k-1$ is an independent parameter in addition to h and p in all finite element processes, thus all finite element computations are dependent on three independent parameters h, p, and k. Just as h and p must be incorporated in local approximations, k also needs to be considered in designing local approximations ${\varphi }_h^e$ so that (3) would yield desired global differentiability of ${\varphi }_h$ at the inter-element boundaries. h, p refinement cannot change or alter k as it is independent of h and p. Thus, we have k-version of finite element method in addition to h- and p-versions. Derivation of higher order global differentiability local approximations for BVPs, their accuracy, convergence behavior and comparisons with published works is the thrust of the work presented in this paper.

1.3. Initial Value Problems (IVPs) and k-Version

Consider an initial value problem

$$\begin{array}{c}\hfill A\varphi (x,t)-f(x,t)=0\forall (x,t)\in {\mathsf{\Omega }}_{xt}={\mathsf{\Omega }}_x\times {\mathsf{\Omega }}_t\end{array}$$

(5)

in which A is a space-time differential operator. Let $2m_1$ and $2m_2$ be the highest orders of derivations of $\varphi$ in space x and time t in the (5), then for a solution of $\varphi$ we have

$$\begin{array}{c}\hfill \varphi \in V\subset H^{(k_1,k_2)}\left({\overline{\mathsf{\Omega }}}_{xt}^{}\right)\end{array}$$

(6)

The scalar product space $H^{(k_1,k_2)}({\overline{\mathsf{\Omega }}}_{xt}^{})$ contains functions that are admissible in (5). $k_1$, $k_2$ are orders of the space in space and time and we have

$$\begin{array}{c}\hfill k_1\ge 2m_1+1;k_2\ge 2m_2+1\end{array}$$

(7)

in which $k_1=2m_1+1$ and $k_2=2m_2+1$ correspond to minimally conforming space in space and time. If we consider space-time coupled finite element processes for obtaining approximation solutions for (5), then we discretize the space-time domain ${\overline{\mathsf{\Omega }}}_{xt}^{}$ into ${\overline{\mathsf{\Omega }}}_{xt}^T=\underset{e}{\stackrel{}{\cup }}{\overline{\mathsf{\Omega }}}_{xt}^e$. ${\overline{\mathsf{\Omega }}}_{xt}^T$ is the discretization of the space-time domain ${\overline{\mathsf{\Omega }}}_{xt}^{}$ in which ${\overline{\mathsf{\Omega }}}_{xt}^e$ is a space-time element. Let ${\varphi }_h^e(x,t)$ be local approximation of $\varphi (x,t)$ over ${\overline{\mathsf{\Omega }}}_{xt}^e$ and ${\varphi }_h(x,t)$ be approximation of $\varphi (x,t)$ over the discretization ${\overline{\mathsf{\Omega }}}_{xt}^T$ such that

$$\begin{array}{cc}\hfill {\varphi }_h(x,t)& =\bigcup ^{}_e{\varphi }_h^e(x,t)\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathrm{and}{\varphi }_h^e(x,t)& =\sum _{i=1}^n{N}_i(x,t){\delta }_i^e\hfill \end{array}$$

(9)

$N_i(x,t)$ are space-time local approximation functions and ${\delta }_i^e$ are nodal dofs for a space-time element e with domain ${\overline{\mathsf{\Omega }}}_{xt}^e$. For the IVP with highest orders of derivatives of orders $2m_1$ and $2m_2$ in space and time require ${\varphi }_h(x,t)$ over ${\overline{\mathsf{\Omega }}}_{xt}^T$ to be at least of class $2m_1$ in space and of class $2m_2$ in time, i.e., $k_1=2m_1+1$ and $k_2=2m_2+1$ correspond to minimally conforming space in space and time. In general we must have $k_1\ge 2m_1+1$ and $k_2\ge 2m_1+1$ and $k_1\to \infty$ and $k_2\to \infty$ are admissible. In general we can write ${\varphi }_h(x,t)$ must be of class $C^{qr}({\overline{\mathsf{\Omega }}}_{xt}^T)$, in which q and r are orders of global differentiability of ${\varphi }_h(x,t)$ in space and time. Hence ${\varphi }_h(x,t)$ over ${\overline{\mathsf{\Omega }}}_{xt}^e$ must be such that ${\varphi }_h={\displaystyle \bigcup _e}{\varphi }_h^e$ yields global differentiability of orders q and r in space and time at the inter-element boundaries. Only then ${\varphi }_h(x,t)$ of class $C^{qr}({\overline{\mathsf{\Omega }}}_{xt}^T)$ is possible. Thus, we see that regardless of BVP or IVP, the local approximation must be designed in such a way that their union yields desired global differentiability at the inter-element boundaries.

2. Literature Review

The importance of the higher order global differentiability of the finite element solutions for BVPs and IVPs have been recognized for quite some time [4,5,6,7,8,9,10]. More recently, Surana and coworkers [11,12,13,14,15,16] pointed out applications in BVPs and IVPs in which higher order global differentiability approximations were shown to be highly meritorious. First formal presentation and introduction of k-version of finite element for BVPs in which k is the order of the approximation space is due to Surana et al., in the three basic papers [1,2,3] devoted to self adjoint, non-selfadjoint and nonlinear differential operators. In these works, variationally consistent (VC) as well as variationally inconsistent (VIC) integral forms [17,18] were used in the model problem studies. In references [19,20,21] Surana et al., presented 2D HGDA/TP as well as HGDA/DG finite element formulations.

Since the publication of the works by Surana et al., on k-version, there have been many attempts with different techniques to address the higher order global differentiability issues primarily in context with the BVPs. In the following we discuss four typical publications that appeared in 2010–2020 that are typical of higher order global differentiability published works. We regret to point out that none of these works acknowledge or reference the large number of published works and two textbooks by Surana et al., that address systematic and general methodologies for deriving HGDA/DG or HGDA/TP in ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ including model problem studies, applications, error estimation and convergence rates. In reference [22], the authors claim to extend tensor product concepts of 1D approximation to derive what they call as $C1$-$Q_k$ family of finite elements. In reference [23] $C^1$ basis functions are presented, that contain same dofs at the corner nodes that arise due to tensor product of 1D $C^1$ functions in natural coordinate space. These cannot be transformed to physical coordinate space [17,18], authors in the paper do this incorrectly. Reference [24] claims the QH8-$C1$ as new innovation in higher order continuity. No reference is made to any of the published works, the work has no mathematical foundation. In reference [25], the choices of the dofs and the geometrical approach used has no mathematical foundation either. There are many other works of similar type in which some approach (mostly without mathematical basis) is used to achieve $C^1$ continuity. Isogeometric analysis [26,27] also claims the method to be k-version. In subsequent sections we clearly demonstrate that this is indeed not the case. Unfortunately we also find the motivation of k-version as presented by Surana et al. [1,2,3] and discussed in this paper is misrepresented in many of the published works on isogeometric analysis as well as other works. We point out this subsequently so that the readers of this paper have better understanding of the motivation for k-version of FEM initiated by Surana et al. First, we make some remarks.

Remark 1.

1.

In all publications only $C^1$ global differentiability is addressed for a fixed degree of polynomials of the functions in the local approximation.

2.

No published work on $C^q$, $C^{qr}$, $C^{qrs}$ (in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$) with $q,r,s>1$ higher order global differentiability local approximations is available.

3.

The published $C^1$ works are not higher degree, hierarchical, i.e., degree of local approximations cannot be increased to any arbitrary value p while maintaining $C^1$ nature of approximation, hence hierarchical nature of approximations is naturally not possible.

4.

In summary there is no unified and sound mathematical framework that exists at present for constructing $C^q$, $C^{qr}$ and $C^{qrs}$; $q,r,s\ge 1$ local approximation with global differentiability in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ that allow the use of arbitrarily higher order degree polynomials and in which the local approximations are hierarchical.

3. Scope and Approach Used in the Present Work

The scope of work published in this paper is summarized in the following.

• The paper presents a unified methodology and mathematical infrastructure for deriving higher order global differentiability approximations of class $C^q$, $C^{qr}$, $C^{qrs}$; $p,q,r\ge 0$ for BVPs and IVPs in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$.
• Since the global differentiability of approximation ${\varphi }_h$ for discretization ${\overline{\mathsf{\Omega }}}_x^T$ or ${\overline{\mathsf{\Omega }}}_{xt}^T$ is due to local approximation ${\varphi }_h^e$ over ${\overline{\mathsf{\Omega }}}_x^e$ or ${\overline{\mathsf{\Omega }}}_{xt}^e$, i.e.,

  $$\begin{array}{c}\hfill {\varphi }_h=\bigcup ^{}_e{\varphi }_h^e\end{array}$$
  Thus, the local approximation ${\varphi }_h^e$ over ${\overline{\mathsf{\Omega }}}_x^e$ or ${\overline{\mathsf{\Omega }}}_{xt}^e$ must be designed in such a way that desired globally differentiability of ${\varphi }_h$ (i.e., of orders q or $qr$ or $qrs$) is achievable at the inter-element boundaries of the discretizations. This approach is essential in all finite element processes. This is due to the fact the union of element approximation must always establish the global approximation ${\varphi }_h$ over the discretization ${\overline{\mathsf{\Omega }}}_x^T$ or ${\overline{\mathsf{\Omega }}}_{xt}^T$.
• The approximations in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ of type $C^q$, $C^{qr}$ and $C^{qrs}$ must be p-version hierarchical, i.e., for chosen q, r, s, we must be able to increase p-levels to whatever value we desire without affecting order of global differentiability q, r, s. If possible, lower p-levels must be embedded (complete subset) in the higher p-levels (hierarchical property or embedding property) so that the computations performed at lower p-levels can be used when doing computations at higher p-level.
• The derivation of the local approximations must always be in natural coordinate space with transparent transformations to physical coordinate space to achieve the desired global differentiability in the physical coordinate space.
• In all derivations of the type $C^q$, $C^{qr}$, $C^{qrs}$ in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ the nodal configuration for geometry and for the dofs must remain the same regardless of the choices of q, r, s and p-levels. This permits use of a single discretization for all studies if h-refinement is not used.
• All derivations must initialize with $C^0$, $C^{00}$ and $C^{000}$p-version hierarchical local approximation of arbitrary p-level in the derivations of $C^q$, $C^{qr}$, $C^{qrs}$ HGDA/DG. The derivations must retain the hierarchical structure of increasing p-level beyond the minimum polynomial degree needed for q, $qr$, $qrs$ orders of global differentiability.

4. ${\mathit{C}}^{\mathit{q}}({\overline{\mathsf{\Omega }}}_{\mathit{x}}^{\mathit{e}})$ Local Approximation in ${\mathbb{R}}^{\mathbf{1}}$

Consider a three node 1D element e with domain ${\overline{\mathsf{\Omega }}}_x^e$ of length $h_e$ in the physical coordinate space x, mapped in the natural coordinate space $\xi$ in a two unit length (Figure 1). Let 1, 2, 3 be local node numbers then

$$\begin{array}{c}\hfill x=x\left(\xi \right)=\sum _{i=1}^3{\widehat{N}}_i\left(\xi \right)x_i\end{array}$$

(10)

describes mapping of a point in $\xi$- and x-spaces. ${\widehat{N}}_i\left(\xi \right)$ are standard quadratic Lagrange functions. Figure 2 shows a typical inter-element boundary between element $e-1$ and e.

For the global approximation ${\varphi }_h={\displaystyle \bigcup _e}{\varphi }_h^e$ to be of class $C^q({\overline{\mathsf{\Omega }}}_x^T)$, the local approximation ${\varphi }_h^e$ over ${\overline{\mathsf{\Omega }}}_x^e$ must be of class $C^q({\overline{\mathsf{\Omega }}}_x^e)$, i.e., the union of local approximations ${\varphi }_h^e$ must yield global differentiability of order $C^q$ at the inter-element boundaries. This requires that we construct local approximation ${\varphi }_h^e$ for an element e such that $\varphi$, $\frac{d^i\varphi }{d{x}^i}$; $i=1,2,\cdots ,q$ are the dofs at the nonhierarchical nodes (boundary nodes) so that at the inter-element boundaries continuity of $\varphi$, $\frac{d^i\varphi }{d{x}^i}$; $i=1,2,\cdots ,q$ will be enforced. We shall see that in the interior $\left({\mathsf{\Omega }}_x^e\right)$ of each element ${\varphi }_h^e$, ${\varphi }_h$ is always of class higher than $C^q$. We begin with a $C^0$p-version hierarchical local approximation that is 1D local approximation in natural coordinate space $\xi$. Using the three node configuration of Figure 1b we can write the following. The p-version hierarchical local approximation (using $\varphi$ instead of ${\varphi }_h^e$) of class $C^0$ and of p-level p is given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \varphi \left(\xi \right)=\hspace{1em}\left({\displaystyle \frac{1-\xi }{2}}\right){\varphi }_1+\left({\displaystyle \frac{1+\xi }{2}}\right){\varphi }_3+\sum _{i=2}^{p_{\xi }}\left({\displaystyle \frac{{\xi }^i-a}{i!}}\right){\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_{\xi =0}\\ \hfill a=1\mathrm{when}i\mathrm{is}\mathrm{even};a=\xi \mathrm{when}i\mathrm{is}\mathrm{odd}.\end{array}\end{array}$$

(11)

Here ${\varphi }_1$ and ${\varphi }_3$ are nodal dofs at nonhierarchical nodes 1 and 3. This local approximation ensures continuity of $\varphi$ (only) at the inter-element boundaries (nodes 1 and 3), hence ${\varphi }_h^e$ is class $C^0({\overline{\mathsf{\Omega }}}_x^e)$ and therefore would yield ${\varphi }_h$ (over ${\overline{\mathsf{\Omega }}}_x^T$) of class $C^0\left({\overline{\mathsf{\Omega }}}_x^T\right)$.

4.1. $C^1({\overline{\mathsf{\Omega }}}_x^e)$ Local Approximation

A $C^1({\overline{\mathsf{\Omega }}}_x^e)$ local approximation will require $\varphi$ and $\frac{d\varphi }{dx}$ as dofs at non hierarchical nodes 1 and 3 (Figure 1). Thus, if we begin with (11), then $\varphi$ is already a nodal dof at nodes 1 and 3 and $\frac{d\varphi }{dx}$ at nodes 1 and 3 can be established by borrowing two terms from the sum in (11) and converting them to $\frac{d\varphi }{dx}$ at nodes 1 and 3. First, we write (11) as (borrowing two terms form the sum in (11)).

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\xi \right)& =\hspace{1em}\left({\displaystyle \frac{1-\xi }{2}}\right){\varphi }_1+\left({\displaystyle \frac{1+\xi }{2}}\right){\varphi }_3+\left({\displaystyle \frac{{\xi }^2-1}{2!}}\right){\left.{\displaystyle \frac{d^2\varphi }{d{\xi }^2}}\right|}_2+\left({\displaystyle \frac{{\xi }^3-\xi }{3!}}\right){\left.{\displaystyle \frac{d^3\varphi }{d{\xi }^3}}\right|}_2\hfill \\ \hfill & +\sum _{i=4}^{p_{\xi }}\left({\displaystyle \frac{{\xi }^i-a}{i!}}\right){\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_2\hfill \end{array}\end{array}$$

(12)

Using (12) we can obtain expressions for ${\left.\frac{d\varphi }{d\xi }\right|}_{\xi =-1}={\left.\frac{d\varphi }{d\xi }\right|}_1$ and ${\left.\frac{d\varphi }{d\xi }\right|}_{\xi =1}={\left.\frac{d\varphi }{d\xi }\right|}_3$ (at nodes 1 and 3). From these two equations we can solve for ${\left.\frac{d^2\varphi }{d{\xi }^2}\right|}_2$, ${\left.\frac{d^3\varphi }{d{\xi }^3}\right|}_2$ and then substitute them back in (12) to obtain

$$\begin{array}{c}\hfill \varphi \left(\xi \right)=N_0^1\left(\xi \right){\varphi }_1+{\underset{~}{N}}_1^1{\left.{\displaystyle \frac{d\varphi }{d\xi }}\right|}_1+N_3^0{\varphi }_3+{\underset{~}{N}}_3^1{\left.{\displaystyle \frac{d\varphi }{d\xi }}\right|}_3+\sum _{i=4}^{p_{\xi }}{N}_2^i\left(\xi \right){\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_{2(\mathrm{or}\xi =0)}\end{array}$$

(13)

in which

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill N_1^0\left(\xi \right)& =\hspace{1em}\left(\frac{(1-\xi )}{2}+\frac{({\xi }^3-\xi )}{4}\right)\hfill \\ \hfill N_1^1\left(\xi \right)& =\hspace{1em}\left(\frac{({\xi }^3-\xi )}{4}-\frac{({\xi }^2-1)}{4}\right)J\hfill \\ \hfill N_2^0\left(\xi \right)& =\hspace{1em}\left(\frac{(1+\xi )}{2}-\frac{({\xi }^3-\xi )}{4}\right)\hfill \\ \hfill N_2^1\left(\xi \right)& =\hspace{1em}\left(\frac{({\xi }^3-\xi )}{4}+\frac{({\xi }^2-1)}{4}\right)J\hfill \\ \hfill N_3^i\left(\xi \right)& =\hspace{1em}\left\{\begin{array}{c}\left(\frac{({\xi }^i-1)-\frac{1}{2}({\xi }^2-1)}{i!}\right)i\mathrm{is}\mathrm{even}\hfill \\ \left(\frac{({\xi }^i-\xi )-\frac{i-1}{2}({\xi }^3-\xi )}{i!}\right)i\mathrm{is}\mathrm{odd}\hfill \end{array}\right.(i=4,5,\dots ,p)\hfill \end{array}\end{array}$$

(14)

we note that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d^i\varphi }{d{\xi }^i}}& ={\left(J\right)}^i{\displaystyle \frac{d^i\varphi }{d{x}^i}};i=1,2,\cdots \hfill \\ \hfill J& ={\displaystyle \frac{dx}{d\xi }}=\frac{he}{2}\left(\mathrm{for}\mathrm{equally}\mathrm{spaced}\mathrm{nodes}\mathrm{in}x\mathrm{space}\right)\hfill \end{array}\end{array}$$

(15)

using (15), ${\left.\frac{d\varphi }{d\xi }\right|}_1$ and ${\left.\frac{d\varphi }{d\xi }\right|}_3$ can be replaced by $J{\left.\frac{d\varphi }{dx}\right|}_1$ and $J{\left.\frac{d\varphi }{dx}\right|}_3$ and we have

$$\begin{array}{c}\hfill N_1^0\left(\xi \right){\varphi }_1+N_1^1{\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_1+N_3^0{\varphi }_3+N_3^1{\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_3+\sum _{i=4}^{p_{\xi }}{N}_2^i\left(\xi \right){\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_{2(\mathrm{or}\xi =0)}\end{array}$$

(16)

in which $N_1^1=J{N}_1^1$ and $n_3^1=J{\underset{~}{N}}_3^1$. $\varphi \left(\xi \right)$ in (16) is the desired p-version local approximation of class $C^1({\overline{\mathsf{\Omega }}}_x^e)$ in which if we choose $p_{\xi }=3$ then there are no dofs at the hierarchical node 2. Each p-level increase beyond p-level of 3 adds one dof at the hierarchical node 2.

4.2. $C^2({\overline{\mathsf{\Omega }}}_x^e)$ Local Approximation

In this case we need $\varphi$, $\frac{d\varphi }{dx}$, $\frac{d^2\varphi }{d{x}^2}$ as dofs at node 1 and 3 of the element of Figure 1. Thus, compared to $C^0({\overline{\mathsf{\Omega }}}_x^e)$ local approximation of (11), we need additional two dofs at nodes 1 and 3. We borrow four terms from the sum in (11).

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\xi \right)& =\hspace{1em}\left(\frac{1-\xi }{2}\right){\varphi }_1 + \left(\frac{1+\xi }{2}\right){\varphi }_2 + \left(\frac{{\xi }^2-1}{2!}\right){\left.\frac{d^2\varphi }{d{\xi }^2}\right|}_{\xi =0} + \left(\frac{{\xi }^3-\xi }{3!}\right){\left.\frac{d^3\varphi }{d{\xi }^3}\right|}_{\xi =0}\hfill \\ \hfill & + \left(\frac{{\xi }^4-1}{4!}\right){\left.\frac{d^4\varphi }{d{\xi }^4}\right|}_{\xi =0} + \left(\frac{{\xi }^5-\xi }{5!}\right){\left.\frac{d^5\varphi }{d{\xi }^5}\right|}_{\xi =0} + \sum _{i=6}^p{N}_3^i\left(\xi \right){\left.\frac{d^i\varphi }{d{\xi }^i}\right|}_{\xi =0}\hfill \end{array}\end{array}$$

(17)

Using (17) we obtain expression for ${\left.\frac{d\varphi }{d\xi }\right|}_{\xi =-1}$, ${\left.\frac{d^2\varphi }{d{\xi }^2}\right|}_{\xi =-1}$ (at node one) and ${\left.\frac{d\varphi }{d\xi }\right|}_{\xi =1}$, ${\left.\frac{d^2\varphi }{d{\xi }^2}\right|}_{\xi =-1}$ (at node 3). Using these four equations we solve for ${\left.\frac{d^i\varphi }{d{\xi }^i}\right|}_{\xi =0}$; $i=2,3,\cdots ,5$ and substitute these back into (17) and use the relationship $\frac{d^i\varphi }{d{\xi }^i}={\left(J\right)}^i\frac{d^i\varphi }{d{\xi }^i}$ to obtain [17,18] (and noting that $\xi =-1$ and $\xi =1$ are locations of nodes 1 and 3).

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\xi \right)& ={\widehat{N}}_1^0\left(\xi \right){\varphi }_1+{\widehat{N}}_1^1\left(\xi \right){\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_1+{\widehat{N}}_1^2{\left.{\displaystyle \frac{d^2\varphi }{d{x}^2}}\right|}_1+{\widehat{N}}_3^0{\varphi }_3+{\widehat{N}}_3^1{\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_3+{\tilde{N}}_3^2\left(\xi \right){\left.{\displaystyle \frac{d^2\varphi }{d{x}^2}}\right|}_3\hfill \\ \hfill & +\sum _{i=6}^{p_{\xi }}{N}_2^i\left(\xi \right){\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_{2(\mathrm{or}\xi =0)}\hfill \end{array}\end{array}$$

(18)

This is the desired $C^2({\overline{\mathsf{\Omega }}}_x^e)$ approximation that would yield global approximation ${\overline{\mathsf{\Omega }}}_x^T$ of class $C^2({\overline{\mathsf{\Omega }}}_x^T)$.

4.3. $C^q({\overline{\mathsf{\Omega }}}_x^e)$ Local Approximation

Following the procedure presented for local approximation of types $C^1({\overline{\mathsf{\Omega }}}^e)$ and $C^2({\overline{\mathsf{\Omega }}}^e)$, we can derive the following for $C^q({\overline{\mathsf{\Omega }}}^e)$ type local approximations in one dimension that ensure continuity of $\varphi$ of order q over the discretization ${\overline{\mathsf{\Omega }}}_x^T$ of ${\overline{\mathsf{\Omega }}}_x$. For $\varphi \left(\xi \right)$ of class $C^q({\overline{\mathsf{\Omega }}}_x^e)$ and $C^q({\overline{\mathsf{\Omega }}}_x^T)$, we can write

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\xi \right)& ={\underset{~}{N}}_1^0\left(\xi \right){\varphi }_1+{\underset{~}{N}}_1^1\left(\xi \right){\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_1+{\underset{~}{N}}_1^2\left(\xi \right){\left.{\displaystyle \frac{d^2\varphi }{d{x}^2}}\right|}_1+\dots +{\underset{~}{N}}_1^q\left(\xi \right){\left.{\displaystyle \frac{d^q{\varphi }_1}{d{x}^q}}\right|}_1\hfill \\ \hfill & +{\underset{~}{N}}_3^0\left(\xi \right){\varphi }_3+{\underset{~}{N}}_3^1\left(\xi \right){\left.{\displaystyle \frac{d\varphi }{dx}}\right|}_3+{\underset{~}{N}}_3^2\left(\xi \right){\left.{\displaystyle \frac{d^2\varphi }{d{x}^2}}\right|}_3+\dots +{\underset{~}{N}}_3^q\left(\xi \right){\left.{\displaystyle \frac{d^q\varphi }{d{x}^q}}\right|}_3\hfill \\ \hfill & +\sum _{j=2q+1}^{p_{\xi }}{\underset{~}{N}}_3^j\left(\xi \right){\left.\frac{d^j\varphi }{d{\xi }^j}\right|}_{\xi =0\left(\mathrm{or}2\right)}\hfill \end{array}\end{array}$$

(19)

The expression for $\varphi \left(\xi \right)$ in (19) allows local approximations that yield global differentiability of order q and will permit increase in p-level up to any arbitrary value $p_{\xi }$ beyond the minimum p-level needed for $C^q({\overline{\mathsf{\Omega }}}_x^e)$.

5. ${\mathit{C}}^{\mathit{qr}}({\overline{\mathsf{\Omega }}}_{\mathit{xy}}^{\mathit{e}})$ Local Approximation in ${\mathbb{R}}^{\mathbf{2}}$

Regarding HGDA or HGDA/DG in ${\mathbb{R}}^2$ there is no single theory, methodology or even a unique and general mathematical procedure in the published works that can be used to derive them. The basic point of confusion begins with the very definition of what $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ local approximations are. From $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$, it is clear that we want global differentiability of orders q and r in x and y, but what is not clear is the fact in context with local approximation how do we achieve this in a unique and hierarchical manner.

Definition 1.

If n and t are normal and tangential directions at an inter-element boundary Γ (Figure 3), then global differentiability of orders q and r of ϕ in the n and t directions requires that

$$\begin{array}{c}\hfill \frac{\partial {}^{i+j}\varphi }{\partial n^i\partial t^j}i=0,1,\cdots ,q;j=0,1,\cdots ,ri+j<q+r\end{array}$$

(20)

be unique at every point Q along the boundary Γ.

Figure 3. A four element discretization: inter-element boundary Γ between elements and .

Figure 3. A four element discretization: inter-element boundary Γ between elements and .

Figure 3 shows a four element discretization in ${\mathbb{R}}^2$. If $\varphi$ is a dependent variable, then the $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ local approximation of continuity require continuity of the derivatives of $\varphi$ in normal $\left(n\right)$ and tangential (t) directions of orders q and r, respectively, along the inter-element boundary, i.e., these derivations must be unique on $\mathsf{\Gamma }$ regardless of whether we use element or . For example along the inter-element boundary $\mathsf{\Gamma }$ between elements a and b the following must hold

$$\begin{array}{c}\hfill \frac{\partial {}^{i+j}{\varphi }_h^a}{\partial n^i\partial t^j}=\frac{\partial {}^{i+j}{\varphi }_h^b}{\partial n^i\partial t^j}i=0,1,\cdots ,qj=0,1,\cdots ,ri+j<q+r\end{array}$$

(21)

in which ${\varphi }_h^a$ and ${\varphi }_h^b$ are local approximation for elements and of the four element discretization. The orthogonal normal and tangential directions n and t can be transformed to x, y orthogonal direction, thus (21) implies that the following must hold along an inter-element boundary $\mathsf{\Gamma }$ between elements a and b.

$$\begin{array}{c}\hfill \frac{\partial {}^{i+j}{\varphi }_h^a}{\partial x^i\partial y^j}=\frac{\partial {}^{i+j}{\varphi }_h^b}{\partial x^i\partial y^j}i=0,1,\cdots ,qj=0,1,\cdots ,ri+j<q+r\end{array}$$

(22)

Thus, in a discretization ${\overline{\mathsf{\Omega }}}_{xy}^T={\displaystyle \bigcup _e}{\overline{\mathsf{\Omega }}}_{xy}^e$, using (nine node p-version) elements, an element ${\overline{\mathsf{\Omega }}}_x^e$ shown in Figure 4 in $xy$ space with its map in natural coordinates space $\xi \eta$ shown in Figure 5, must have the nodal dofs at the corner nodes for (nonhierarchical nodes) $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ and $C^{qr}({\overline{\mathsf{\Omega }}}_{\xi \eta }^{})$ differentiability approximations in $xy$ and $\xi \eta$ spaces as shown in Figure 4 and Figure 5. We note that local approximations of class $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ require minimum p-levels of $2q+1$ and $2r+1$. For p-levels $p_{\xi }=2q+1$ and $p_{\eta }=2r+1$, the hierarchical nodes have no dofs. For $p_{\xi }>2q+1$ and $p_{\eta }>2r+1$ the dofs at the nonhierarchical corner nodes remain unaffected but additional dofs appear at the hierarchical nodes. It is more meaningful to represent dofs in terms of nodal variable operators (obtained by removing dependent variable). When the nodal variable operators act on a dependent variable they produce dofs. These are shown in Figure 4 and Figure 5 in $xy$ and $\xi \eta$ coordinate spaces. Nodal variable operators for the nonhierarchical nodes of a nine node p-version element in $xy$ and $\xi \eta$ spaces for $C^{11}$, $C^{22}$ and $C^{33}$ classes are listed in

Table 1. Nodal variable operators at the nonhierarchical nodes for $C^{11}$, $C^{22}$, $C^{33}$ classes in $xy$ and $\xi \eta$ spaces.

Global Differentiability		Nodal Variable Operators at Nonhierarchical Nodes (Nodes 1, 3, 5, 7)
		$\mathit{xy}$ Space	$\mathit{\xi }\mathit{\eta }$ Space
1	$C^{11}$	1, $\frac{\partial }{\partial x}$, $\frac{\partial }{\partial y}$	1, $\frac{\partial }{\partial \xi }$, $\frac{\partial }{\partial \eta }$
2	$C^{22}$ when (1) holds	$\frac{\partial {}^2}{\partial x^2}$, $\frac{\partial {}^2}{\partial y^2}$, $\frac{\partial {}^2}{\partial y\partial x}$	$\frac{\partial {}^2}{\partial {\xi }^2}$, $\frac{\partial {}^2}{\partial {\eta }^2}$, $\frac{\partial {}^2}{\partial \eta \partial \xi }$
3	$C^{33}$ when (1) and (2) hold	$\frac{\partial {}^3}{\partial x^3}$, $\frac{\partial {}^3}{\partial y^3}$, $\frac{\partial {}^3}{\partial y\partial x^2}$, $\frac{\partial {}^3}{\partial y^2\partial x}$	$\frac{\partial {}^3}{\partial {\xi }^3}$, $\frac{\partial {}^3}{\partial {\eta }^3}$, $\frac{\partial {}^3}{\partial \eta \partial {\xi }^2}$, $\frac{\partial {}^3}{\partial {\eta }^2\partial \xi }$

.

Remark 2.

1.

We note that dofs for $C^{11}$, $C^{22}$ and $C^{33}$ listed in Table 1 at the nonhierarchical nodes in $\xi \eta$ space can be transformed into those listed in $xy$ space. Transformation of the dofs at nonhierarchical nodes is always possible in general between $C^{qr}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ and $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ (see Section 6).

2.

Based on Remark 1, it is convenient to generate $C^{qr}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ approximation first (nodal dofs or nodal variable operators and associated functions) and then transform the nodal dofs or nodal variable operators at the nonhierarchical nodes into $xy$ space for ${\overline{\mathsf{\Omega }}}_{xy}^e$. We keep in mind that all derivations initiate with p-version $C^{00}$ hierarchical local approximation in $\xi \eta$ spaces, as this approach will allow us to retain the feature of increasing p-levels beyond ($2q+1$) and $(2r+1)$ in ξ and η directions without influencing the dofs at the nonhierarchical nodes responsible for $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ continuity at the inter-element boundaries and will also retain the hierarchical nature of $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ approximation.

5.1. Transformation of Dofs between $C^{qr}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ and $C^{qr}\left({\overline{\mathsf{\Omega }}}_{xy}^{}\right)$ at the Nonhierarchical Nodes

Let

$$\begin{array}{c}\hfill \begin{array}{c}\hfill x(\xi ,\eta )=\sum _i{\tilde{N}}_i(\xi ,\eta )x_i\\ \hfill y(\xi ,\eta )=\sum _i{\tilde{N}}_i(\xi ,\eta )y_i\end{array}\end{array}$$

(23)

define a mapping of points between $\xi \eta$ and $xy$ configurations of a nine node p-version hierarchical element. For suitably chosen mapping, the lengths between $\xi \eta$ and $xy$ are defined by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}dx\\ dy\end{array}\right\}=\left[\begin{array}{cc}{\displaystyle \frac{\partial x}{\partial \xi }}& {\displaystyle \frac{\partial x}{\partial \eta }}\\ {\displaystyle \frac{\partial y}{\partial \xi }}& {\displaystyle \frac{\partial y}{\partial \eta }}\end{array}\right]\left\{\begin{array}{c}d\xi \\ d\eta \end{array}\right\}=\left[\begin{array}{cc}{x}_{\xi }& x_{\eta }\\ y_{\xi }& y_{\eta }\end{array}\right]\left\{\begin{array}{c}d\xi \\ d\eta \end{array}\right\}=\left[J\right]\left\{\begin{array}{c}d\xi \\ d\eta \end{array}\right\}\end{array}$$

(24)

in which

$$\begin{array}{c}\hfill \left[J\right]=\left[\begin{array}{cc}{x}_{\xi }& x_{\eta }\\ y_{\xi }& y_{\eta }\end{array}\right]\end{array}$$

(25)

Let

$$\begin{array}{ll}{\left\{\varphi \right\}}_1^{\xi \eta }& ={\left[\begin{array}{c}{\displaystyle \frac{\partial \varphi }{\partial \xi }}{\displaystyle \frac{\partial \varphi }{\partial \eta }}\end{array}\right]}^T;{\left\{\varphi \right\}}_1^{xy}={\left[\begin{array}{c}{\displaystyle \frac{\partial \varphi }{\partial x}}{\displaystyle \frac{\partial \varphi }{\partial y}}\end{array}\right]}^T\hfill \end{array}$$

(26)

$$\begin{array}{ll}{\left\{\varphi \right\}}_2^{\xi \eta }& ={\left[\begin{array}{c}{\displaystyle \frac{{\partial }^2\varphi }{\partial {\xi }^2}}{\displaystyle \frac{{\partial }^2\varphi }{\partial \eta \partial \xi }}{\displaystyle \frac{{\partial }^2\varphi }{\partial {\eta }^2}}\end{array}\right]}^T;{\left\{\varphi \right\}}_2^{xy}={\left[\begin{array}{c}{\displaystyle \frac{{\partial }^2\varphi }{\partial x^2}}{\displaystyle \frac{{\partial }^2\varphi }{\partial y\partial x}}{\displaystyle \frac{{\partial }^2\varphi }{\partial y^2}}\end{array}\right]}^T\hfill \end{array}$$

(27)

$$\begin{array}{ll}{\left\{\varphi \right\}}_3^{\xi \eta }& ={\left[\begin{array}{c}{\displaystyle \frac{{\partial }^3\varphi }{\partial {\xi }^3}}{\displaystyle \frac{{\partial }^3\varphi }{\partial \eta \partial {\xi }^2}}{\displaystyle \frac{{\partial }^3\varphi }{\partial {\eta }^2\partial \xi }}{\displaystyle \frac{{\partial }^3\varphi }{\partial {\eta }^3}}\end{array}\right]}^T;{\left\{\varphi \right\}}_3^{xy}={\left[\begin{array}{c}{\displaystyle \frac{{\partial }^3\varphi }{\partial x^3}}{\displaystyle \frac{{\partial }^3\varphi }{\partial y\partial x^2}}{\displaystyle \frac{{\partial }^3\varphi }{\partial y^2\partial x}}{\displaystyle \frac{{\partial }^3\varphi }{\partial y^3}}\end{array}\right]}^T\hfill \end{array}$$

(28)

In the following we define the rule of transformation between ${\left\{\varphi \right\}}_i^{\xi \eta }$ and ${\left\{\varphi \right\}}_i^{xy}$; $i=1,2,\cdots$

$$\begin{array}{c}\hfill {\left\{\varphi \right\}}_i^{\xi \eta }=\left[J_i\right]{\left\{\varphi \right\}}_i^{xy}\end{array}$$

(29)

Obviously,

$$\begin{array}{c}\hfill \left[J_1\right]=\left[\begin{array}{cc}{x}_{\xi }& y_{\xi }\\ x_{\eta }& y_{\eta }\end{array}\right]={\left[J\right]}^T\end{array}$$

(30)

Using the following notation, we can transform dofs between natural and physical coordinate spaces [17,20].

$$\begin{array}{c}\hfill x_{\xi }={\displaystyle \frac{\partial x}{\partial \xi }};x_{\xi }^i={\left({\displaystyle \frac{\partial x}{\partial \xi }}\right)}^i;x_{{\xi }^i}={\displaystyle \frac{{\partial }^{i}x}{\partial {\xi }^i}};x_{{\xi }^i{\eta }^j}={\displaystyle \frac{{\partial }^{i+j}x}{\partial {\xi }^i\partial {\eta }^j}}\end{array}$$

(31)

$$\begin{array}{c}\hfill x_{\eta }={\displaystyle \frac{\partial x}{\partial \eta }};x_{\eta }^i={\left({\displaystyle \frac{\partial x}{\partial \eta }}\right)}^i;x_{{\eta }^i}={\displaystyle \frac{{\partial }^{i}x}{\partial {\eta }^i}}\end{array}$$

(32)

$$\begin{array}{c}\hfill y_{\xi }={\displaystyle \frac{\partial y}{\partial \xi }};y_{\xi }^i={\left({\displaystyle \frac{\partial y}{\partial \xi }}\right)}^i;y_{{\xi }^i}={\displaystyle \frac{{\partial }^{i}y}{\partial {\xi }^i}};y_{{\xi }^i{\eta }^j}={\displaystyle \frac{{\partial }^{i+j}y}{\partial {\xi }^i\partial {\eta }^j}}\end{array}$$

(33)

$$\begin{array}{c}\hfill y_{\eta }={\displaystyle \frac{\partial y}{\partial \eta }};y_{\eta }^i={\left({\displaystyle \frac{\partial y}{\partial \eta }}\right)}^i;y_{{\eta }^i}={\displaystyle \frac{{\partial }^{i}y}{\partial {\eta }^i}}\end{array}$$

(34)

and

$$\begin{array}{c}\hfill {\left\{\varphi \right\}}_2^{xy}={\left[J_2\right]}^{-1}\left[{\left\{\varphi \right\}}_2^{\xi \eta }-\left[J_2^1\right]{\left\{\varphi \right\}}_1^{xy}\right]\end{array}$$

(35)

$$\begin{array}{c}\hfill \left[J_2\right]=\left[\begin{array}{ccc}{x}_{\xi }^2& 2x_{\xi }{y}_{\xi }& y_{\xi }^2\\ x_{\xi }{x}_{\eta }& x_{\xi }{y}_{\eta }+x_{\eta }{y}_{\xi }& y_{\eta }{y}_{\xi }\\ x_{\eta }^2& 2x_{\eta }{y}_{\eta }& y_{\eta }^2\end{array}\right]\left[J_2^1\right]=\left[\begin{array}{cc}{x}_{{\xi }^2}& y_{{\xi }^2}\\ x_{\xi \eta }& y_{\xi \eta }\\ x_{{\eta }^2}& y_{{\eta }^2}\end{array}\right]\end{array}$$

(36)

$$\begin{array}{c}\hfill {\left\{\varphi \right\}}_3^{xy}={\left[J_3\right]}^{-1}\left[{\left\{\varphi \right\}}_3^{\xi \eta }-\left[J_3^1\right]{\left\{\varphi \right\}}_2^{xy}-\left[J_3^2\right]{\left\{\varphi \right\}}_1^{xy}\right]\end{array}$$

(37)

$$\begin{array}{c}\hfill \left[J_3\right]=\left[\begin{array}{cccc}{x}_{\xi }^3& 3x_{\xi }^2{y}_{\xi }& 3x_{\xi }{y}_{\xi }^2& y_{\xi }^3\\ x_{\xi }^2{x}_{\eta }& 2x_{\eta }{x}_{\xi }{y}_{\xi }+x_{\xi }^2{y}_{\eta }& 2x_{\xi }{y}_{\eta }{y}_{\xi }+y_{\xi }^2{x}_{\eta }& y_{\eta }{y}_{\xi }^2\\ x_{\eta }^2{x}_{\xi }& 2x_{\eta }{x}_{\xi }{y}_{\eta }+x_{\eta }^2{y}_{\xi }& 2x_{\eta }{y}_{\eta }{y}_{\xi }+y_{\eta }^2{x}_{\xi }& y_{\xi }{y}_{\eta }^2\\ x_{\eta }^3& 3x_{\eta }^2{y}_{\eta }& 3x_{\eta }{y}_{\eta }^2& y_{\eta }^3\end{array}\right]\end{array}$$

(38)

$$\begin{array}{c}\hfill \left[J_3^1\right]=\left[\begin{array}{ccc}3x_{\xi }{x}_{{\xi }^2}& 3x_{{\xi }^2}{y}_{\xi }+3y_{{\xi }^2}{x}_{\xi }& 3y_{\xi }{y}_{{\xi }^2}\\ x_{\eta }{x}_{{\xi }^2}+2x_{\xi }{x}_{\xi \eta }& x_{{\xi }^2}{y}_{\eta }+2x_{\xi \eta }{y}_{\xi }+x_{\eta }{y}_{{\xi }^2}+2x_{\xi }{y}_{\xi \eta }& y_{\eta }{y}_{{\xi }^2}+2y_{\xi }{y}_{\xi \eta }\\ x_{\xi }{x}_{{\eta }^2}+2x_{\eta }{x}_{\xi \eta }& x_{\xi }{y}_{{\eta }^2}+2x_{\xi \eta }{y}_{\eta }+x_{{\eta }^2}{y}_{\xi }+2x_{\eta }{y}_{\xi \eta }& y_{{\eta }^2}{y}_{\xi }+2y_{\eta }{y}_{\xi \eta }\\ 3x_{\eta }{x}_{{\eta }^2}& 3x_{{\eta }^2}{y}_{\eta }+3y_{{\eta }^2}{x}_{\eta }& 3y_{\eta }{y}_{{\eta }^2}\end{array}\right]\end{array}$$

(39)

$$\begin{array}{c}\hfill \left[J_3^2\right]=\left[\begin{array}{cc}{x}_{{\xi }^3}& y_{{\xi }^3}\\ x_{{\xi }^2\eta }& y_{{\xi }^2\eta }\\ x_{\xi {\eta }^2}& y_{\xi {\eta }^2}\\ x_{{\eta }^3}& y_{{\eta }^3}\end{array}\right]\end{array}$$

(40)

From (29), (35) and (37) it follows that we can write the following general expression.

$$\begin{array}{c}\hfill {\left\{\varphi \right\}}_i^{xy}={\left[J_i\right]}^{-1}\left[{\left\{\varphi \right\}}_i^{\xi \eta }-\left[J_i^1\right]{\left\{\varphi \right\}}_{i-1}^{xy}-\left[J_i^2\right]{\left\{\varphi \right\}}_{i-2}^{xy}-\dots -\left[J_i^{i-1}\right]{\left\{\varphi \right\}}_1^{xy}\right]\end{array}$$

(41)

Let us use the abbreviation HGDA/DG to denote higher order global differentiability approximation. The symbol $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ denotes a HGDA/DG of class q and r in respective coordinates x and y on ${\overline{\mathsf{\Omega }}}_{xy}^e$. To show specific details that are needed in (41), we consider $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$, $C^{22}({\overline{\mathsf{\Omega }}}_{xy}^e)$, $C^{33}({\overline{\mathsf{\Omega }}}_{xy}^e),\cdots$ local approximations.

5.2. 2D $C^{11}$ HGDA/DG Quadrilateral Elements

The dofs at the hierarchical and the nonhierarchical nodes of a 2D $C^{11}$ HGDA/DG elements are shown in Figure 6. A $C^{00}$p-version element will only have $\varphi$ as dof at the nonhierarchical nodes. The nonhierarchical nodes of 2D $C^{11}$ HGDA/DG element require additional eight dofs compared to 2D $C^{00}$ p-version element. These are generated by borrowing eight dofs from the hierarchical nodes of $C^{00}$ p-version element (two from each hierarchical edge or face node), i.e., we borrow ${\varphi }_{{\xi }^2}$, ${\varphi }_{{\xi }^3}$ and ${\varphi }_{{\eta }^2}$, ${\varphi }_{{\eta }^3}$ from nodes (2,6) and (4,8) of $C^{00}$ p-version element (see [17,18,19,20,21] for more details).

5.3. $C^{22}$ HGDA/DG for 2D Distorted Quadrilateral Elements in $xy$ Space

Figure 7 shows dofs at the nonhierarchical and the hierarchical nodes of this element. To generate ${\varphi }_x$, ${\varphi }_y$, ${\varphi }_{x^2}$, ${\varphi }_{xy}$, ${\varphi }_{y^2}$ at the nonhierarchical nodes (a total of twenty), we borrow ${\varphi }_{{\xi }^2}$, ${\varphi }_{{\xi }^3}$, ${\varphi }_{{\xi }^4}$, ${\varphi }_{{\xi }^5}$ from nodes 2 and 6 and ${\varphi }_{{\eta }^2}$, ${\varphi }_{{\eta }^3}$, ${\varphi }_{{\eta }^4}$, ${\varphi }_{{\eta }^5}$ from nodes 4 and 8, a total of sixteen. Additional four dofs are borrowed from node 9 of 2D $C^{00}$ p-version element (shown in Figure 8). See [17,18,19,20,21] for more details and for the rationale of this choice.

5.4. 2D $C^{33}$ HGDA/DG Elements

For this element ${\varphi }_x$, ${\varphi }_y$, ${\varphi }_{x^2}$, ${\varphi }_{xy}$, ${\varphi }_{y^2}$, ${\varphi }_{x^3}$, ${\varphi }_{x^{2}y}$, ${\varphi }_{x{y}^2}$, and ${\varphi }_{y^3}$ are the dofs needed at the nonheirarchical nodes in addition to $\varphi$. Compared to $C^{00}$ p-version element, we need nine additional dofs at each of the nonhierarchical nodes (a total of thirty six). We borrow ${\varphi }_{{\xi }^2}$, ${\varphi }_{{\xi }^3}$, ${\varphi }_{{\xi }^4}$, ${\varphi }_{{\xi }^5}$, ${\varphi }_{{\xi }^6}$ and ${\varphi }_{{\xi }^7}$ from each of the four mid side hierarchical nodes (a total of twenty four). The remaining twelve dofs are borrowed from node nine as shown in Figure 8. See [17,18,19,20,21] for the rationale of this choice and for more details.

5.5. 2D $C^{ij}$ HGDA/DG Elements

The dofs $\left\{{\delta }^e\right\}$ of $C^{00}$ p-version element are partitioned into $\left\{{\delta }_{co}^e\right\}$, $\left\{{\delta }_m^e\right\}$, $\left\{{\delta }_c^e\right\}$ and $\left\{{\delta }_{mc}^e\right\}$ in which $co$, m, c and $mc$ stand for corner nodes, mid side nodes, center nodes and mid side plus center node and we can write the following:

$$\begin{array}{c}\hfill \varphi (\xi ,\eta )=\left[a\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}+\left[b\right]{\left\{{\delta }_{mc}^e\right\}}_{el}+\left[c\right]{\left\{{\delta }_m^e\right\}}_{r_2}+\left[d\right]{\left\{{\delta }_c^e\right\}}_{r_3}\end{array}$$

(42)

in which $r_1,r_2,r_3$ imply retained, el implies eliminated and $\left[a\right]$, $\left[b\right]$, $\left[c\right]$, and $\left[d\right]$ are row matrices of approximation functions in $\xi ,\eta$ space for 2D $C^{00}$ p-version element. As an example, consider a 2D $C^{11}$ HGDA/DG element. In this case, from nodes 2, 4, 6 and 8, we have:

$$\begin{array}{c}\hfill {\left\{{\delta }_m^e\right\}}_e=[{\varphi }_{{\xi }^2}{|}_2,{\varphi }_{{\xi }^3}{|}_2,{\varphi }_{{\eta }^2}{|}_4,{\varphi }_{{\eta }^3}{|}_4,{\varphi }_{{\xi }^2}{|}_6,{\varphi }_{{\xi }^3}{|}_6,{\varphi }_{{\eta }^2}{|}_8,{\varphi }_{{\eta }^3}{|}_8]\end{array}$$

(43)

We note that for 2D $C^{11}$ HGDA/DG element, we do not need to borrow dofs from node 9 of 2D $C^{00}$ p-version element. Hence, ${\left\{{\delta }_{mc}^e\right\}}_{el}={\left\{{\delta }_m^e\right\}}_{el}$.

If ${\left\{{\delta }^e\right\}}_n^{xy}$ are the new dofs at the nonhierarchical nodes of 2D $C^{ij}$ HGDA/DG element, then for 2D $C^{ij}$ HGDA/DG we have (in the natural coordinate space) the following at nodes 1, 3, 5 and 7.

$$\begin{array}{c}\hfill {\left\{{\delta }^e\right\}}_n^{\xi \eta }={\left\{{\varphi }_{\xi }{|}_1,{\varphi }_{\eta }{|}_1,{\varphi }_{\xi }{|}_3,{\varphi }_{\eta }{|}_3,{\varphi }_{\xi }{|}_5,{\varphi }_{\eta }{|}_5,{\varphi }_{\xi }{{|}_7,{\varphi }_{\eta }|}_7\right\}}^T\end{array}$$

(44)

We differentiate (42) with respect to $\xi$ and $\eta$ (as needed) and evaluate these at each of the non hierarchical nodes to obtain the following:

$$\begin{array}{c}\hfill {\left\{{\delta }^e\right\}}_n^{\xi \eta }=\left[A\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}+\left[B\right]{\left\{{\delta }_{mc}^e\right\}}_{el}+\left[C\right]{\left\{{\delta }_m^e\right\}}_{r_2}+\left[D\right]{\left\{{\delta }_c^e\right\}}_{r_3}\end{array}$$

(45)

Solving for the dofs to be eliminated, i.e., ${\left\{{\delta }_{mc}^e\right\}}_{el}$ in (45), we obtain

$$\begin{array}{c}\hfill {\left\{{\delta }_{mc}^e\right\}}_{el}={\left[B\right]}^{-1}{\left\{{\delta }^e\right\}}_n^{\xi \eta }-{\left[B\right]}^{-1}\left[A\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}-{\left[B\right]}^{-1}\left[C\right]{\left\{{\delta }_m^e\right\}}_{r_2}-{\left[B\right]}^{-1}\left[D\right]{\left\{{\delta }_c^e\right\}}_{r_3}\end{array}$$

(46)

Substituting Jacobian of transformation from (41) into the above equation, we can transform the new derivative dofs from $\xi \eta$ space to $xy$ space. Equation (46) can thus be written as

$$\begin{array}{c}\hfill {\left\{{\delta }_{mc}^e\right\}}_{el}={\left[B\right]}^{-1}\left[J_i\right]{\left\{{\delta }^e\right\}}_n^{xy}-{\left[B\right]}^{-1}\left[A\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}-{\left[B\right]}^{-1}\left[C\right]{\left\{{\delta }_m^e\right\}}_{r_2}-{\left[B\right]}^{-1}\left[D\right]{\left\{{\delta }_c^e\right\}}_{r_3}\end{array}$$

(47)

Now, substituting ${\left\{{\delta }_{mc}^e\right\}}_{el}$ from (47) into (42),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi (\xi ,\eta )& =\left[a\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}+\left[b\right]({\left[B\right]}^{-1}\left[J_i\right]{\left\{{\delta }^e\right\}}_n^{xy}-{\left[B\right]}^{-1}\left[A\right]{\left\{{\delta }_{co}^e\right\}}_{r_1}\hfill \\ \hfill & -{\left[B\right]}^{-1}\left[C\right]{\left\{{\delta }_m^e\right\}}_{r_2}-{\left[B\right]}^{-1}\left[D\right]{\left\{{\delta }_c^e\right\}}_{r_3})+\left[c\right]{\left\{{\delta }_m^e\right\}}_{r_2}+\left[d\right]{\left\{{\delta }_c^e\right\}}_{r_3}\hfill \end{array}\end{array}$$

(48)

Collecting terms in the (48), we obtain the final form of the $C^{ij}$ HGDA/DG local approximations as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi (\xi ,\eta )& =\left(\left[a\right]-\left[b\right]{\left[B\right]}^{-1}\left[A\right]\right){\left\{{\delta }_{co}^e\right\}}_{r_1}+\left[b\right]{\left[B\right]}^{-1}\left[J_i\right]{\left\{{\delta }^e\right\}}_n^{xy}\hfill \\ \hfill & +\left(\left[c\right]-\left[b\right]{\left[B\right]}^{-1}\left[C\right]\right){\left\{{\delta }_m^e\right\}}_{r_2}+\left(\left[d\right]-\left[b\right]{\left[B\right]}^{-1}\left[D\right]\right){\left\{{\delta }_c^e\right\}}_{r_3}\hfill \end{array}\end{array}$$

(49)

Remark 3.

1.

The derivation is based on 2D $C^{00}$ p-version hierarchical approximation.

2.

Dofs from the hierarchical nodes of $C^{00}$ element are borrowed to generate desired degrees of freedom first in $\xi \eta$ space which are then transformed to $xy$ space.

3.

All matrices $\left[a\right],\left[b\right],\left[c\right]$ and $\left[d\right]$ contain $C^{00}p$-version hierarchical local approximations while matrices $\left[A\right],\left[B\right],\left[C\right]$ and $\left[D\right]$ contain derivatives of 2D $C^{00}$ p-version hierarchical local approximations.

4.

Dofs at the non hierarchical nodes of 2D $C^{ij}$ HGDA/DG are transformable between $\xi \eta$ and $xy$ spaces transparently (an intrinsic feature of truly $C^{ij}$ elements).

5.6. HGDA/TP: Higher Order Global Differentiability Elements Using Tensor Product

Consider a 2D distorted nine node p-version hierarchical element in $xy$ space and its map in $\xi \eta$ space in a two unit square (Figure 9). We know that tensor product requires the two orthogonal directions in which the 1D approximations are defined, hence tensor product cannot be used in $xy$ space for the distorted 2D element. However, in $\xi \eta$ space 2D functions and the dofs for $C^{qr}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ approximated can be derived using tensor product of 1D $C^q\left({\overline{\mathsf{\Omega }}}_{\xi }\right)$ and $C^r\left({\overline{\mathsf{\Omega }}}_{\eta }\right)$ approximations. For the $C^{qr}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ 2D element in $\xi \eta$ space the nodal dofs at the nonhierarchical nodes are function $\varphi$ and its derivatives of $\varphi$ with respect to $\xi$, $\eta$. A $C^{qr}\left({\overline{\mathsf{\Omega }}}_{xy}^e\right)$ element will require transformation of these dofs to the$xy$space. We consider details in the following. Consider the 1D $C^q\left({\overline{\mathsf{\Omega }}}_{\xi }\right)$ and $C^r\left({\overline{\mathsf{\Omega }}}_{\eta }\right)$ hierarchical local approximation in the $\xi$ and $\eta$ directions (following (19))

$$\begin{array}{c}\hfill \varphi \left(\xi \right)={}^{\xi }{N}_1^0{\varphi }_1+{}^{\xi }{N}_3^0{\varphi }_3+\sum _{i=1}^q{}^{\xi }{N}_1^i{\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_1+\sum _{i=1}^q{}^{\xi }{N}_3^i{\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_3+\sum _{i=2q+1}^{p_{\xi }}{}^{\xi }{N}_2^j{\left.{\displaystyle \frac{d^i\varphi }{d{\xi }^i}}\right|}_2\end{array}$$

(50)

$$\begin{array}{c}\hfill \varphi \left(\eta \right)={}^{\eta }{N}_1^0{\varphi }_1+{}^{\eta }{N}_3^0{\varphi }_3+\sum _{j=1}^r{}^{\eta }{N}_1^j{\left.{\displaystyle \frac{d^j\varphi }{d{\eta }^j}}\right|}_1+\sum _{j=1}^r{}^{\eta }{N}_3^j{\left.{\displaystyle \frac{d^j\varphi }{d{\eta }^j}}\right|}_3+\sum _{j=2r+1}^{p_{\eta }}{}^{\eta }{N}_2^j{\left.{\displaystyle \frac{d^j\varphi }{d{\eta }^j}}\right|}_2\end{array}$$

(51)

using (50) and (51) we can generate 2D approximation functions in the $\xi \eta$ space as well as nodal dofs (also in $\xi \eta$ space) by taking tensor product of 1D approximation and the tensor product of the nodal variable operators [17,18]. The dofs at the nonhierarchical nodes (1, 3, 5, 7) generated in doing so are given by

$$\begin{array}{c}\hfill \frac{\partial {}^{i+j}\varphi }{\partial {\xi }^i\partial {\eta }^j};i=0,1,\cdots ,q;j=0,1,\cdots ,r\end{array}$$

(52)

For $C^{11}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ we have the following dofs at nodes 1, 3, 5, 7

$$\begin{array}{c}\hfill \varphi ,\frac{\partial \varphi }{\partial \xi },\frac{\partial \varphi }{\partial \eta },\frac{\partial {}^2\varphi }{\partial \eta \partial \xi }\end{array}$$

(53)

For $C^{22}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ we have the following dofs at nodes 1, 3, 5, 7

$$\begin{array}{c}\hfill \varphi ,\frac{\partial \varphi }{\partial \xi },\frac{\partial \varphi }{\partial \eta },\frac{\partial {}^2\varphi }{\partial \eta \partial \xi };\frac{\partial {}^2\varphi }{\partial {\xi }^2},\frac{\partial {}^2\varphi }{\partial {\eta }^2},\frac{\partial {}^3\varphi }{\partial \eta \partial {\xi }^2},\frac{\partial {}^3\varphi }{\partial {\eta }^2\partial \xi },\frac{\partial {}^4}{\partial \eta {\xi }^2}\end{array}$$

(54)

Based on Section 6, all dofs in (53) and (54) with respect into $\xi \eta$ cannot be transformed to those with respect to $xy$.

5.7. Remarks: HGDA or HGDA/DG and HGDA/TP Elements

• The issue of not being able to transform the derivative dofs (with respect to $\xi$ and $\eta$) at nonhierarchical nodes needs further investigation. We note that (Section 5) a $C^{qr}\left({\overline{\mathsf{\Omega }}}_{xy}^e\right)$ approximation has

  $$\begin{array}{c}\hfill \frac{\partial {}^{i+j}\varphi }{\partial {\xi }^i{\eta }^j};i=0,1,\cdots ,q;j=0,1,\cdots ,r;i+j<q+r\end{array}$$

  (55)

  dofs at the nonhierarchical nodes. These can be transformed from $\xi \eta$ space to $xy$ space as shown in Section 5.
• The dofs at the nonhierarchical nodes in the tensor product element in $\xi \eta$ space are

  $$\begin{array}{c}\hfill \frac{\partial {}^{i+j}\varphi }{\partial {\xi }^i{\eta }^j};i=0,1,\cdots ,q;j=0,1,\cdots ,r\end{array}$$

  (56)

  for each order of continuity (i.e., $C^{11}$, $C^{22}$, ⋯), (56) contain additional dofs compared to (55). We note that for $C^{11}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$, sum of the orders of derivatives with respect to $\xi$ and $\eta$ should be one which is true in case of (55), but in (56) $\frac{\partial {}^2\varphi }{\partial \eta \partial \xi }$ term violates this condition confirming that HGDA/TP is not a true $C^{11}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$ approximation in $\xi$ and $\eta$. The same argument holds true for $C^{22}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$, as in (54). Thus, the tensor product in $\xi \eta$ space does not generate dofs at the nonhierarchical nodes for $C^{ij}\left({\overline{\mathsf{\Omega }}}_{\xi \eta }^{}\right)$; $i=1,2,\cdots ,q$, $j=1,2,\cdots ,r$ local approximation that conform to (55), containing only the needed dofs at the nonhierarchical nodes for $C^{ij}$ continuity.
• If we only consider rectangular elements in $xy$ space with $\xi \eta$ for each element parallel to $xy$ and pointing in the same direction, then we can transform $\frac{d^i\varphi }{d{\xi }^i}$ and $\frac{d^j\varphi }{d{\eta }^j}$ in (50) and (51) to $\frac{d^i\varphi }{d{x}^i}$ and $\frac{d^j\varphi }{d{y}^j}$ using

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d^i\varphi }{d{\xi }^i}}=\left(J_x^i\right){\displaystyle \frac{d^i\varphi }{d{x}^i}}& \mathrm{and}{\displaystyle \frac{d^j\varphi }{d{\eta }^j}}=\left(J_y^j\right){\displaystyle \frac{d^j\varphi }{d{y}^j}}\hfill \\ \hfill J_x={\displaystyle \frac{dx}{d\xi }}& \mathrm{and}{J}_y={\displaystyle \frac{dy}{d\eta }}\hfill \end{array}\end{array}$$

  (57)

  giving

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\xi \right)& ={}^{\xi }{N}_1^0\varphi +{}^{\xi }{N}_3^0{\varphi }_3\hfill \\ \hfill & +\sum _{i=1}^q\left({}^{\xi }{N}_1^i\left(J_x\right){}^i\right){\left.{\displaystyle \frac{d^i\varphi }{d{x}^i}}\right|}_1+\sum _{i=1}^q\left({}^{\xi }{N}_3^i\left(J_x\right){}^i\right){\left.{\displaystyle \frac{d^i\varphi }{d{x}^i}}\right|}_3+\sum _{i=2q+1}^{p_{\xi }}{}^{\xi }{N}_2^i{\left.{\displaystyle \frac{d^i\varphi }{d{x}^i}}\right|}_2\hfill \end{array}\end{array}$$

  (58)

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill \varphi \left(\eta \right)& ={}^{\eta }{N}_1^0\varphi +{}^{\eta }{N}_3^0{\varphi }_3\hfill \\ \hfill & +\sum _{j=1}^r\left({}^{\eta }{N}_1^j\left(J_y\right){}^j\right){\left.{\displaystyle \frac{d^j\varphi }{d{y}^j}}\right|}_1+\sum _{j=1}^r\left({}^{\eta }{N}_3^j\left(J_y\right){}^j\right){\left.{\displaystyle \frac{d^j\varphi }{d{y}^j}}\right|}_3+\sum _{j=2r+1}^{p_{\eta }}{}^{\eta }{N}_2^j{\left.{\displaystyle \frac{d^j\varphi }{d{y}^j}}\right|}_2\hfill \end{array}\end{array}$$

  (59)
  By taking tensor product of expressions (58) and (59), the resulting element will contain

  $$\begin{array}{c}\hfill \frac{\partial {}^{i+j}\varphi }{\partial {\xi }^i\partial y^j};i=0,1,\cdots ,q;j=0,1,\cdots ,r\end{array}$$

  (60)

  as dofs at the nonhierarchical corner nodes (nodes 1, 3, 5 and 7). We refer to this element as $C^{qr}({\overline{\mathsf{\Omega }}}_{xy}^e)$ HGDA/TP element. This element is not a true $C^{ij}({\overline{\mathsf{\Omega }}}_{xy}^e)$ global differentiability element.
• HGDA/TP is not very useful in applications as it requires elements to be rectangular in $xy$ space with $\xi \eta$ pointing in the same directions as $xy$. Thus, distorted element geometries required for an irregular domain cannot be supported by this higher order global differentiability 2D finite elements based on tensor product.
• If we consider a discretization of a domain ${\overline{\mathsf{\Omega }}}_{xy}^{}$ for which HGDA/TP are valid, then for this special case we may ask the following question. Are the HGDA/DG formulations meritorious over HGDA/TP or vice versa? Answering this question requires many basic and important considerations.

  (a)

  In the case of HGDA/DG elements, the dofs at the nonhierarchical nodes (corner nodes) in $\xi \eta$ and $xy$ spaces for each order of continuity $C^{ii}$; $i=1,2,\cdots$ constitute derivatives sets in the two spaces that can be transformed from $\xi \eta$ to $xy$ and $xy$ to $\xi \eta$ spaces. Thus, these dofs are in conformity with calculus.

  (b)

  HGDA/DG elements contain only those dofs at the nonhierarchical nodes for the desired orders of continuity that are essential. Thus, a $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$ element will have

  $$\begin{array}{c}\hfill \varphi ,\frac{\partial \varphi }{\partial dx},\frac{\partial \varphi }{\partial y}\end{array}$$

  as dofs at the nonhierarchical nodes. The set $\frac{\partial \varphi }{\partial x}$, $\frac{\partial \varphi }{\partial y}$ can be obtained from $\frac{\partial \varphi }{\partial \xi }$, $\frac{\partial \varphi }{\partial \eta }$ in $\xi \eta$ using calculus.

  (c)

  In the case of a $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$ HGDA/TP element we have

  $$\begin{array}{c}\hfill \varphi ,\frac{\partial \varphi }{\partial x},\frac{\partial \varphi }{\partial y},\frac{\partial {}^2\varphi }{\partial y\partial x}\end{array}$$

  as dofs at the nonhierarchical nodes. The first important point to note is that $\frac{\partial {}^2\varphi }{\partial y\partial x}$ is not part of the $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$ approximation. This is obvious from the dofs at the nonhierarchical nodes of HGDA/DG $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$ element.

  (d)

  Does the extra dof $\frac{\partial {}^2\varphi }{\partial y\partial x}$ at each of the four nonhierarchical node in case of a HGDA/TP element help in improving local approximation or just adds extra dofs without much improvement compared to the dofs in HGDA/DG $C^{11}({\overline{\mathsf{\Omega }}}_x^e)$ element? From the numerical studies presented in a later section we conclusively show that the latter is the case, i.e., $\frac{\partial {}^2\varphi }{\partial y\partial x}$ dofs do not result in an improved approximation but add additional dofs compared to HGDA/DG $C^{11}({\overline{\mathsf{\Omega }}}_{xy}^e)$ approximation, hence they result in poor accuracy for a given dofs and also result in lower convergence rates. A somewhat less rigorous but persuasive argument may be since $\frac{\partial {}^2\varphi }{\partial y\partial x}$ is only one element of the complete set needed for $C^{22}({\overline{\mathsf{\Omega }}}_{xy}^{})$, hence cannot be very effective in improving local approximation. Similar arguments hold for $C^{22}({\overline{\mathsf{\Omega }}}_{xy}^e),\cdots$ etc. HGDA/DG and HGDA/TP elements.

  (e)

  Thus, our conclusion is that in the HGDA/TP approach of deriving local approximation:

  • true $C^{iJ}({\overline{\mathsf{\Omega }}}_{xy}^e)$; $i=1,2\cdots ,q$, $j=1,2\cdots ,r$ local approximations are not possible.
  • HGDA/TP elements are inferior to HGDA/DG elements in all aspects.

6. HGDA/DG: Triangular and Hexadral Elements

Higher order global differentiability local approximations for triangular and hexadron elements can also be derived [17,18,20,21], but are not presented here for the sake of brevity.

7. Isogeometric Analysis and $\mathit{k}$-Version

In 2005 [26,27] isogeometric analysis method was presented. According to the authors Analogues of finite element h- and p-refinement schemes are presented and a new, more efficient, higher order concept, k-refinement is introduced. Thus, in our view isogeometric analysis claims to be an h, p, k method. We first present details of the isogeometric analysis, then compare it with $hpk$ finite element method, but, more importantly, we compare it with the k-version of finite element introduced by Surana et al. [1,2,3,11,12,13,14,19,20,21].

In isogeometric analysis, the complex domain is subdivided into subdomains (not finite elements). Let us refer to a subdomain as a patch. For a patch, one picks points called control points that are not necessarily on the boundary of the patch. Using the coordinates (geometric locations) of these points one constructs a geometric description (say using NURBS or other methods). This gives a relationship of the type for geometry.

$$\begin{array}{c}\hfill \left\{{\left(x\right)}_p^i\right\}=\left[N\left(\xi \right)\right]\left\{\begin{array}{c}\left\{x_p^i\right\}\\ \left\{y_p^i\right\}\\ \left\{z_p^i\right\}\end{array}\right\}\end{array}$$

(61)

Equation (61) gives position coordinates $x_p^i$, $y_p^i$, $z_p^i$ of every point in the patch in terms of the control points and the functions $\left[N\right]$. A relationship similar to (61) is established for displacements u, v, w of an arbitrary point in the patch at location x, y, z in terms of displacements of the control points and the same functions as used in (61)

$$\begin{array}{c}\hfill {\mathit{u}}_p^i=\left\{u_p^i\right\}=\left[N\right]\left\{\begin{array}{c}\left\{u_p^i\right\}\\ \left\{v_p^i\right\}\\ \left\{w_p^i\right\}\end{array}\right\}=\left[N\right]\left\{{\delta }_p^i\right\}\end{array}$$

(62)

$\left\{{\delta }_p^i\right\}$ are the total dofs for the patch i. If the approximation (62) is of degree p in x, y, z then it is also of class $C^p\left({\overline{\mathsf{\Omega }}}_p^i\right)$ where ${\overline{\mathsf{\Omega }}}_p^i$ is the domain of the ith patch. Thus, within the patch the order of the approximation space of $\mathit{u}$ is $k=p+1$. We remark that so far we do not have finite elements in the true sense in which ${\left({\overline{\mathsf{\Omega }}}_p^i\right)}^T=\bigcup _e{\overline{\mathsf{\Omega }}}^e$, in which ${\left({\overline{\mathsf{\Omega }}}_p^i\right)}^T$ is a discretization of ${\overline{\mathsf{\Omega }}}_p^i$. Instead (62) is a statement that one would use in classical methods of approximation in which ${\overline{\mathsf{\Omega }}}_p$ is not a discretized domain (patch). Now, since we have an approximation of displacements in (62), we can proceed with standard method using a method of approximation to obtain algebraic equation for the patch ${\overline{\mathsf{\Omega }}}_p^i$. If we choose linear elasticity as an example, then the balance of linear momenta gives

$$\begin{array}{c}\hfill \mathit{A}·\mathit{u}-\mathit{f}=0\end{array}$$

(63)

We construct an integral form of (63) over the domain $\overline{\mathsf{\Omega }}$ consisting of N patches, i.e., $\overline{\mathsf{\Omega }}={\displaystyle \bigcup _N^{i=1}}{\overline{\mathsf{\Omega }}}_p^i$. Then

$$\begin{array}{c}\hfill \left\{u_N\right\}=\underset{}{\bigcup ^N}\left\{u_p^i\right\}\end{array}$$

where $\left\{u_N\right\}$ is approximation of u, v, w over $\overline{\mathsf{\Omega }}$ and $\left\{u_p^i\right\}$ is approximation of $\left\{u\right\}$ over ${\overline{\mathsf{\Omega }}}_p^i$ (ith patch). Based on the fundamental lemma [17,18] we can write the following using (63)

$$\begin{array}{c}\hfill {(\mathit{A}·{\mathit{u}}_N-\mathit{f},\mathit{v})}_{\overline{\mathsf{\Omega }}}=0\end{array}$$

(64)

in which $\overline{\mathsf{\Omega }}={\displaystyle \bigcup _i}{\overline{\mathsf{\Omega }}}_p^i$. If we consider Galerkin method with weak form, then $\mathit{v}=\delta {\mathit{u}}_N$, (64) can be written as (over the patches).

$$\begin{array}{c}\hfill \sum _{i=1}^N{(\left[A\right]\left\{u_p^i\right\}-\left\{f\right\},\left\{v\right\})}_{{\overline{\mathsf{\Omega }}}_p^i}=0;v=\delta u_p^i\end{array}$$

(65)

Consider ${(\left[A\right]\left\{u_p^i\right\}-\left\{f\right\},\left\{v\right\})}_{{\overline{\mathsf{\Omega }}}_p^i}$ for ith patch using integration by parts (IBP) (since in linear elasticity adjoint of $\mathit{A}$, i.e., ${\mathit{A}}^{*}=\mathit{A}$), we obtain

$$\begin{array}{c}\hfill {(\left[A\right]\left\{u_p^i\right\}-\left\{f\right\},\left\{v\right\})}_{{\overline{\mathsf{\Omega }}}_p^i}=B^p(u_p^i,v)-l^p\left(v\right)\end{array}$$

(66)

Substituting (62) in (66) and noting that $\delta u_p^i={\left[N\right]}^T$ we obtain

$$\begin{array}{c}\hfill {(\left[A\right]\left\{u_p^i\right\}-\left\{f\right\},\left\{v\right\})}_{{\overline{\mathsf{\Omega }}}_p^i}=\left[{\int }_{{\overline{\mathsf{\Omega }}}_p^i}{\left[B\right]}^T\left[D\right]\left[B\right]d\mathsf{\Omega }\right]\left\{{\delta }_p^i\right\}-\left\{F_p^i\right\}=\left[K_p^i\right]\left\{{\delta }_p^i\right\}-\left\{F_p^i\right\}\end{array}$$

(67)

Secondary variables due to IBP are included in $\left\{F_p^i\right\}$ and strain $\left\{\epsilon \right\}=\left[B\right]\left\{\delta \right\}$, stress $\left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}$. In addition, $\left[K_p^i\right]$ is the stiffness matrix for the ith patch. Using (67) in (65), we have

$$\begin{array}{c}\hfill \left[K\right]\left\{\delta \right\}=\left\{F\right\}\end{array}$$

(68)

in which

$$\begin{array}{cc}\hfill \left[K\right]=\sum _{i=1}^N\left[K_p^i\right];& \mathrm{assembly}\mathrm{of}\mathrm{patch}\mathrm{stiffness}\mathrm{matrices}\hfill \\ \hfill \left\{\delta \right\}=\underset{i=1}{\bigcup ^N}\left\{{\delta }_p^i\right\};& \mathrm{Total}\mathrm{dofs}\mathrm{for}\overline{\mathsf{\Omega }}=\underset{i=1}{\bigcup ^N}{\overline{\mathsf{\Omega }}}_p^i\hfill \end{array}$$

Remark 4.

1.

In this process described above used in isogeometric analysis there is no concept of finite element anywhere.

2.

A patch appears like a $C^0$ p-version finite element in which p-level can be as high as desired. Naturally, if ${\overline{\mathsf{\Omega }}}^e$ is the domain of $C^0$ p-version finite element with degree of local approximation p for a field variable ϕ then ${\varphi }_h^e$ (local approximation of ϕ over ${\overline{\mathsf{\Omega }}}^e$) is naturally of class $C^p\left({\overline{\mathsf{\Omega }}}^e\right)$ and holds over the element domain ${\overline{\mathsf{\Omega }}}^e$. This is not a discovery. This has been known since the advent of the finite element method. The isogeometric analysis with patches is exactly like $C^0$ p-version finite element mesh in which each finite element is a patch. Within each $C^0$ element (i.e., a $C^0$ patch) the solution is of class $C^p$ but at the inter-element boundaries (same as inter patch boundaries) the solutions remains of class $C^0$. Thus, all isogeometric solutions with more than one patch are globally of class $C^0$.

3.

Since only displacements are dofs on the boundaries of the patches, $\left\{u_N\right\}={\displaystyle \bigcup _{i=1}^N}\left\{u_p^i\right\}$ remains of class $C^0$ at the interpatch boundaries. Thus, in all isogeometric analyses containing more than one patch, the global solutions are not of class p (order of the space $k=\raisebox{1ex}{$(p+1)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$), but are undoubtedly of class $C^0$ due to their $C^0$ nature at the interpatch boundaries. Thus, presenting isogeometric analysis as k-version is misleading and misrepresentating.

4.

Computation of the patch stiffness matrix requires integrals (2D or 3D) over an irregular domain ${\overline{\mathsf{\Omega }}}_p^i$ (in (67)). This can be done in several different maps. For example, ${\overline{\mathsf{\Omega }}}_p^i$ is subdivided and each subdivision is mapped into two unit square or two unit cube to integrate using Gauss quadrature while still maintaining (61) and (62) for each subdivision. The subdivisions are not finite elements but this is necessary to do to perform integration over ${\overline{\mathsf{\Omega }}}_p^i$. The subdivision are not finite elements as subdivisions have no concept of local approximation.

5.

The works published by Surana et al., [1,2,3,12,19,20,21] are based on true finite element methodology in which global differentiability of order $k-1$ is always ensured. There is no comparison of these works with isogeometric analysis as at present it cannot do what k-version [1,2,3,19,20,21] can, i.e., to ensure global differentiability of any desired order. Isogeometric at present always has global differentiability of class $C^0$, hence no different than traditional $C^0$ finite element method analysis.

6.

On another small note, in some published works on isogeometric analysis and in other works we find that some authors beleive that work as referenced in [1,2,3,19,20,21] is motivated due to least squares method. This is wrong. Perhaps a consequence of not reading our published works carefully enough. We always point out that minimally conforming space is dictated by the differential operator and not the integral form as the integral form may be a weak form.

7.

We also wish to point out that in all of the published works on isogeometric analysis, the influence of higher order differentiability of geometry on the accuracy of the computations has never been shown. Is it that the benefits advocated are only because of higher order global differentiability of the displacement approximation over the interior of the patches (${\overline{\mathsf{\Omega }}}_p^i$) and that the geometry has nothing to do with this? In our own works on true k-version finite elements with h and p, this is true. As a matter of scientific curiosity, this should have been addressed at the onset of the method, but continues to be ignored.

8. A Priori Error Estimation and Convergence Rates: $\mathit{h}$, $\mathit{p}$, $\mathit{k}$-Versions

Surana et al. [28] have shown that variational consistency [17,18] of the integral forms in BVP and IVP is essential in the derivation of a priori error estimates that establish the dependence of various error norms on h, p, k. This was a fundamental and essential discovery that enabled derivation of a priori error estimates for self-adjoint, non-self adjoint, and nonlinear differential operators. Current published works on a priori error estimates rely on the best approximation property in the $B(·,·)$ norm that only exists for self-adjoint operators. Thus, in the currently published works a priori error estimates are rarely considered and found for non-self-adjoint and non-linear operators. In this paper we only present final results that are useful in model problem studies. Derivation can be found in the references [17,18,28]. Surana et al., showed that for BVPs the following hold when the integral forms used in finite element processes are variationally consistent.

When the differential operator is self-adjoint, GM/WF and LSP yields a VC integral form. For non self-adjoint and non-linear differential operators only LSP based on residual functional yields a VC integral form. Following references [17,28] we have the following a priori error estimates when the integral form is VC regardless of the method used.

Let $A\varphi -f=0$ in ${\mathsf{\Omega }}_x$ be the BVP and let ${\varphi }_h$ be the finite element solution of $A\varphi -f=0$ over ${\overline{\mathsf{\Omega }}}_x^T={\displaystyle \bigcup _e}{\overline{\mathsf{\Omega }}}_x^e$, discretization of ${\overline{\mathsf{\Omega }}}_x^{}$, thus we have

$$\begin{array}{c}\hfill {\parallel e\parallel }_{L_2}=\parallel \varphi -{\varphi }_h{\parallel }_{L_2}\le c_1{h}^{p+1}{\left|\varphi \right|}_{p+1}\end{array}$$

(69)

$$\begin{array}{c}\hfill |\varphi -{\varphi }_h{|}_{H^q}\le c_2{h}^{p+1-q}{\left|\varphi \right|}_{p+1};\mathrm{seminorm}\mathrm{of}\mathrm{order}q\end{array}$$

(70)

$$\begin{array}{c}\hfill {\parallel e\parallel }_{H^q}={\parallel \varphi -\varphi \parallel }_{H^q}\le c_3{h}^{p+1-q}{\left|\varphi \right|}_{p+1}\end{array}$$

(71)

and if

$$\begin{array}{cc}\hfill E=A{\varphi }_h-f;& \mathrm{residual}\mathrm{function}\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill I={(E,E)}_{{\overline{\mathsf{\Omega }}}_x^{}};& \mathrm{residual}\mathrm{functional}\hfill \end{array}$$

(73)

then

$$\begin{array}{c}\hfill {\parallel E\parallel }_{L^2}=\sqrt{I}=c_4{h}^{p+1-2m}{\left|\varphi \right|}_{p+1}\end{array}$$

(74)

Here e is true error, $c_1$, $c_2$, $c_3$ and $c_4$ are variables that only depend on the order k of the approximation and $2m$ is the highest order of the derivatives in $A\varphi -f=0$.

Convergence Rate

Consider ${\parallel e\parallel }_{H^q}$ in (71). Taking logarithm of both sides

$$\begin{array}{cc}\hfill & {log\parallel e\parallel }_{H^q}\le log(c_3{\left|\varphi \right|}_{p+1})+(p+1-q)logh\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill \mathrm{or}& y\le c+mx\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill \mathrm{with}& y={log\parallel e\parallel }_{H^q}\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill & c=log(c_3\parallel \varphi {\parallel }_{p+1})\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill & m=p+1-q\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill & x=logh\hfill \end{array}$$

(80)

When we use the equality in (76), this is an equation of a straight line in $xy$ space in which m is the slope and $\underset{~}{c}$ is the intercept., i.e., if we plot $logh$ versus ${log\parallel e\parallel }_{H^q}$ on an $xy$ plot, then we obtain a straight line whose slope is ($p+1-q$) and intercept is $log\left(c_3\right|\varphi {|}_{p+1})$. Slope ($p+1-q$) is called the rate of convergence of ${\parallel e\parallel }_{H^q}$. Higher values of ($p+1-q$) imply faster convergence of ${\varphi }_h$ to $\varphi$ measured in ${\parallel e\parallel }_{H^q}$. Equation (77) can also be expressed in terms of total dofs for the discretization. If $h=max\left(h_e\right)$, $h_e$ being characteristic length for an element e, then

$$\begin{array}{cc}\hfill h& \propto {\displaystyle \frac{1}{\mathrm{dof}}}\mathrm{or}h=\mathrm{O}\left({\displaystyle \frac{1}{\mathrm{dof}}}\right)\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill \mathrm{using}h& ={\displaystyle \frac{1}{\mathrm{dof}}}\hfill \end{array}$$

(82)

$$\begin{array}{cc}\hfill {log\parallel e\parallel }_{H^q}& \le log\left(c_3\right|\varphi {|}_{p+1})-(p+1-q)log\left(\mathrm{dof}\right)\hfill \end{array}$$

(83)

Remark 5.

1.

We keep in mind that the a priori error estimate (75) only holds for a uniform mesh refinement or at the most for a quasi-uniform mesh refinement.

2.

The quantity ${\parallel e\parallel }_{H^q}$ requires knowledge of the theoretical solution ϕ which is generally not possible for a practical application.

3.

When the approximation spaces are minimally conforming, i.e., when $k=2m+1$, $2m$ is the order of the highest derivation in the BVP $A\varphi -f=0$, then the integrals over discretization ${\overline{\mathsf{\Omega }}}_x^T$ are Riemann and the ${\parallel E\parallel }_{L^2}$ can be computed for ${\overline{\mathsf{\Omega }}}_x^T$ as it does not require knowledge of theoretical solution and we have

$$\begin{array}{c}\hfill {log\parallel E\parallel }_{L^2}\le log(c_4{\left|\varphi \right|}_{p+1})-(p+1-2m)log\left(\mathrm{dof}\right)\end{array}$$

(84)

If the solution ϕ is sufficiently smooth, we can also use $k=2m$. In this case the integrals are Lebesgue over discretization ${\overline{\mathsf{\Omega }}}_x^T$.

4.

In the case of 1D problems (82) holds precisely. However in case of BVP in ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ this is generally not true precisely. For example, for exactly same discretization (h fixed) for HGDA/DG and HGDA/TP elements, the HGDA/TP elements have an additional dof at the nonhierarchical nodes, thus the total dof for HGDA/DG and HGDA/TP are going to be different even though h is same. Thus, even though a priori error estimates are derived using h, but it is perhaps more illustrative to exhibit the real convergence rates if we use dof instead of h. Thus, in all model problem studies we present graphs of error norms versus dof.

9. Model Problems: Numerical Studies (BVPs)

In this section we present numerical studies using higher order continuity 1D p-version hierarchical formulations for: (1) axial deformation of a rod (self adjoint differential operators), (2) 1D convection diffusion equation (non-self adjoint differential operator), (3) and 1D Burgers equation (non-linear differential operator). We also present numerical studies using HGDA/DG and HGDA/TP higher order continuity 2D element formulations for: (1) Poisson’s equation (self-adjoint operator) and (2) 2D convection equation (non-self-adjoint operator) and (3) 2D Burgers Equation (non-linear operator). Various aspects of k-version are illustrated in terms of correct description of physics in the computations, improved accuracy and convergence rates.

9.1. 1D Axial Deformation

The dimensionless form of the 1D deformation of an axial rod is described by ($Au-f=0$) (also see [1])

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -{\displaystyle \frac{d}{dx}}& \left(b\left(x\right){\displaystyle \frac{du}{dx}}\right)=f\left(x\right)0<x<L\hfill \\ \hfill u& \left(0\right)=0{\left({\displaystyle \frac{du}{dx}}\right)}_{x=L}=0;\mathrm{on}\mathrm{boundary}\mathsf{\Gamma }\hfill \end{array}\end{array}$$

(85)

For this BVP $A^{*}$, the adjoint of the differential operator is the same as $A=-\frac{d}{dx}\left(b\left(x\right)\frac{d}{dx}\right)$. A is naturally linear and the concomitant $<Au,v{>}_{\mathsf{\Gamma }}$ is zero as the BCs in (85) are homogeneous, hence A is self-adjoint. We consider finite element formulations based on Galerkin method with weak form (GM/WF) as well as based on residual functional, least squares finite element formulation. Details are straight forward and can be found in [17]. The approximation $u_h$ of u over discretization ${\overline{\mathsf{\Omega }}}_x^T={\displaystyle \bigcup _e}{\overline{\mathsf{\Omega }}}_x^e$ must be such that

$$\begin{array}{c}\hfill u_h\in V\subset H^k\left({\overline{\mathsf{\Omega }}}_x^T\right);k\ge 3\end{array}$$

(86)

in which $k=3$ is minimally conforming space (in which case $u_h$, $\frac{d{u}_h}{dx}$ and $\frac{d^2{u}_h}{d{x}^2}$ are continuous and $u_h$ and $\frac{d{u}_h}{dx}$ are differentiable for $\forall x\in {\overline{\mathsf{\Omega }}}_x^T$). This requires that the local approximation $u_h^e$ over ${\overline{\mathsf{\Omega }}}_x^e$ must also belong to the same space as $u_h$. Figure 10 shows a schematic of the problem.

9.1.1. Analytical Solution

We choose dimensionless L to be one and $b\left(x\right)$ to be constant and the following for the body force $f\left(x\right)$

$$\begin{array}{c}\hfill f\left(x\right)=\alpha \left[{(x+\sigma )}^{\theta }-{(1+\sigma )}^{\theta })\right]\end{array}$$

(87)

where $\sigma$ is a positive constant $0<\sigma <\infty$ and $\theta$ is another constant $-\infty <\theta <\infty$. Then (85) can be written as

$$\begin{array}{c}\hfill -{\displaystyle \frac{d^{2}u}{d{x}^2}}=\beta \left[{(x+\sigma )}^{\theta }-{(1+\sigma )}^{\theta }\right]\mathrm{where}\beta ={\displaystyle \frac{\alpha }{b}}\end{array}$$

(88)

The analytical solution of (88) can be obtained by integrating it twice and using the boundary conditions. We obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill u\left(x\right)& =\beta \left\{\left[{\displaystyle \frac{1}{2}}{x}^2{(1+\sigma )}^{\theta }-{\displaystyle \frac{{(x+\sigma )}^{2+\theta }}{(1+\theta )(2+\theta )}}\right]-\left[{(1+\sigma )}^{\theta }-{\displaystyle \frac{{(1+\sigma )}^{1+\theta }}{(1+\theta )}}\right]\right.\hfill \\ \hfill & +\left.\left[{\displaystyle \frac{{\sigma }^{2+\theta }}{(1+\theta )(2+\theta )}}\right]\right\}\hfill \end{array}\end{array}$$

(89)

$$\begin{array}{cc}\hfill {\displaystyle \frac{du}{dx}}& =\beta \left\{\left[x{(1+\sigma )}^{\theta }-{\displaystyle \frac{{(x+\sigma )}^{1+\theta }}{(1+\theta )}}\right]-\left[{(1+\sigma )}^{\theta }-{\displaystyle \frac{{(1+\sigma )}^{1+\theta }}{(1+\theta )}}\right]\right\}\hfill \end{array}$$

(90)

Remark 6.

1.

For $0<\sigma <\infty$, $f\left(x\right)$ is analytic and square integrable.

2.

When $\sigma =0$, $\frac{du}{dx}$ has square root singularity at $x=0$, hence this case is ruled out.

3.

For $0<\sigma <1$, $\theta =-1.5$, $f\left(x\right)\in C^{\infty }[0,1]$. Hence, $u\left(x\right)\in C^{\infty }[0,1]$.

4.

We use $\alpha =0.001,\theta =-1.5$ and $\sigma =0.0005$. As σ decreases, magnitude of $\frac{d^{i}u}{d{x}^i}$; $i=1,2,\cdots$ increase. Largest value is at $x=0$ i.e., $\frac{d^{i}u}{d{x}^i}$ becomes localized near $x=0$.

9.1.2. Numerical Results

Since the operator A is self-adjoint, the quadratic functional $\pi$ is given by

$$\begin{array}{c}\hfill \pi ={\displaystyle \frac{1}{2}}B(u_h,u_h)-l\left(u_h\right)=\sum _e\left({\displaystyle \frac{1}{2}}{B}^e(u_h^e,u_h^e)-l^e\left(u_h^e\right)\right)\end{array}$$

(91)

and the error or the residual functional is given by

$$\begin{array}{c}\hfill I={(E,E)}_{{\overline{\mathsf{\Omega }}}_x^T}=\sum _e{(E^e,E^e)}_{{\overline{\mathsf{\Omega }}}_x^e}\end{array}$$

(92)

in which $E^e=A{u}_h^e-f\forall x\in {\overline{\mathsf{\Omega }}}_x^e$. For $k\ge 3$ both $\pi$ and I can be computed for GM/WF and LSP. Using theoretical value of $\pi$ we can calculate percentage error in $\pi$, ${}^{g_{\epsilon }}\pi$ and ${}^{l_{\epsilon }}\pi$ for GM/WF and for LSP. We consider $\widehat{L}=10$ in., area of section $\widehat{a}=1{\mathrm{in}.}^2$, and Young’s modulus $\widehat{E}=30\times {10}^6$ psi. The rod is fixed at the left end ($\widehat{x}=0$), and free at the right end ($\widehat{x}=\widehat{L}$) (see Figure 10). If we choose $L=L_0=10$ and $E_0=E=30\times {10}^6={\sigma }_0$ then $F_0=E_0{L}_0^2=30\times {10}^8$ lb, and $a=\frac{\widehat{a}}{L_0^2}=0.01$.

We consider a fixed eight element graded mesh as shown in Figure 10b. Solutions are calculated using GM/WF as well as LSM. For each order of space, p-level is increased up to 19.

Plots of ${}^{g_{\epsilon }}\pi$ and ${}^{l_{\epsilon }}\pi$ versus dofs for GM/WF and LSM are shown in Figure 11a,b. Plots of residual functional I versus dof for GM/WF and LSM are shown in Figure 11c,d.

9.1.3. Discussion of Results

Since discretization is fixed, characteristic length h of the mesh is constant. Thus, from Figure 11a–d, we can observe p-convergence, k-convergence and $pk$-convergence of the solutions and the functionals.

p-convergence: For a fixed k, the sequence of computed solutions for progressively increasing p-level represent p-convergence of the solution and the functionals in Figure 11a–d. For each value of k in Figure 11a–d, we observe progressively increasing slope of the functionals versus dof graphs for progressively increasing p-level, indicating progressively increasing convergence rate (as derived in the a priori error estimation).

k-convergence: We note from the graphs of functionals versus dof, for increasing p-level, that with progressively increasing k the graphs for all k-values remain almost parallel to each other, but progressively shift downwards with progressively increasing value of k. Indicating that, increasing k does not result in higher convergence rate but results in better accuracy (lower value of the functionals) for a given degrees of freedom. Thus, k affects accuracy but not the convergence rate. Along the lines of constant p, hence constant h, p and only k varies, thus we have k-convergence. We note that increase in k results in increasing functional values indicating increasing error. This is due to the fact that for fixed p, increasing k results in reduced dofs (obvious from the graphs). We observe from the graphs that for a fixed dofs, increase in k unconditionally results in lower functional values, hence better accuracy of the solution.

pk-convergence: We can view $pk$-convergence at least in the following two ways: (1) for a fixed k we can choose p based on $p=2k-1$, minimum p-level for k. If we connect these points on graphs in Figure 11a–d, we obtain much higher convergence rate than p-convergence for any value of k. (2) The most obvious dramatic and the highest convergence rates are achieved for combination of p and k for which dofs are almost constant indicated by almost vertical lines on the graphs in Figure 11a–d. The numerical studies presented here for this model problem confirm $h,p,k$ as three independent parameters in the finite element computational processes for self-adjoint operators.

9.2. 1D Convection Diffusion Equation

Consider the following BVP [2]:

$$\begin{array}{c}\hfill {\displaystyle \frac{d\varphi }{dx}}-\frac{1}{Pe}{\displaystyle \frac{d^2\varphi }{d{x}^2}}=0\forall x\in {\mathsf{\Omega }}_x=(0,1)\end{array}$$

(93)

$$\begin{array}{c}\hfill \varphi \left(0\right)=1.0\varphi \left(1\right)=0\end{array}$$

(94)

The differential operator $A=\left(\frac{d}{dx}\right)-\left(\frac{1}{Pe}\right)\left(\frac{d^2}{d{x}^2}\right)$ in (93) is non-self-adjoint, hence GM/WF will yield a VIC integral form but the integral form based on LSP is VC. Here, we only consider LSP (see [17] for numerical studies using GM/WF). We consider $k\ge 2$ and p-levels up to 19. A theoretical solution of (93) and () is given by

$$\begin{array}{c}\hfill \varphi \left(x\right)=\frac{\left(e^{xPe}-e^{Pe}\right)}{\left(1-e^{Pe}\right)}\end{array}$$

(95)

Obviously, $\varphi$ is of class $C^{\infty }$. We consider the following meshes for $Pe=100$ and ${10}^6$.

For $Pe=100$: Element Lengths (5 element mesh)

$$\begin{array}{c}\hfill 0.826545,1.535 \times {10}^{-1},1.535 \times {10}^{-2},1.535 \times {10}^{-2},1.535 \times {10}^{-2}\end{array}$$

For $Pe={10}^6$: Element Lengths (10 element mesh)

$$\begin{array}{c}\hfill 0.722215,2.5 \times {10}^{-1},2.5 \times {10}^{-2},2.5 \times {10}^{-3},2.5 \times {10}^{-4},2.5 \times {10}^{-5},\mathrm{four}\mathrm{of}\mathrm{length}2.5 \times {10}^{-6}\end{array}$$

We consider $k\ge 2$ with p-levels up to 19 in LSFEP. For each k, p-levels are increased up to 19 uniformly for each element of the discretization.

Discussion of Results

As in previous model problem, here also characteristic discretization length h is constant as the discretization is fixed, hence we can study p-convergence, k-convergence and $pk$-convergence of the solution and functionals. Figure 12a,b present residual functional I versus dof graphs for various values of k for both $Pe$.

p-convergence: For a fixed k, the sequence of solutions for progressively increasing p-level represent p-convergence of the solution and the residual functional I in Figure 12a,b. For each value of k in Figure 12a,b we observe progressively increasing slope of residual functional I versus dof graphs for progressively increasing p-levels, indicating progressively increasing convergence rate (as shown in a priori error estimation section).

k-convergence: From the graphs of residual functional I versus dof for increasing p-level but for a fixed value of k, we note that with progressively increasing k the I verus dofs graphs for all k values remain almost parallel to each other but shift downwards with increasing k. This confirms that increase in k does not result in higher convergence rate but results in better accuracy (due to lower values of functional I) for a given degrees of freedom. Thus, k affects accuracy but not the convergence rate. In Figure 12a,b, from the curves of constant p, hence constant h and p, thus only k is varying along these curves (k-convergence). We note that increase in k results in increase in the value of residual functional I, indicating increasing error. This is due to the fact that for fixed p, increasing k results in reduced dofs (obvious from the graphs). We observe from the graphs that for a fixed dofs, increase in k unconditionally results in lower values of I, hence better accuracy.

pk-convergence: As in case of previous model problem, here also we can view $pk$-convergence in at least the following two ways: (1) for a fixed k, we can choose $p=2k-1$, minimum p-level of this choice of k. If we construct these parts on the graph in Figure 12a,b, we observe much higher convergence rate than p-convergence for any value of k. (2) The most dramatic and the highest convergence rates are achieved for combinations of p and k for which dofs are almost constant indicated by nearly vertical lines in the graphs in Figure 12a,b.

The numerical studies presented here for the model problem in which differential operator is non-self adjoint also confirm that $h,p$ and k are three independent parameters in finite element computations and the observed convergence behavior is in accordance with a priori estimates.

9.3. 1D Burgers Equation

Steady state Burgers equation ($A\varphi -f=0$) is given by (also see [3]):

$$\begin{array}{cc}\hfill \varphi & {\displaystyle \frac{d\varphi }{dx}}-\frac{1}{Re}{\displaystyle \frac{d^2\varphi }{d{x}^2}}=0\forall x\in {\mathsf{\Omega }}_x=(0,1)\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill \varphi & \left(0\right)=1.0\varphi \left(1\right)=0\hfill \end{array}$$

(97)

A theoretical solution of (96) and (97) is given by [2,17,18].

$$\begin{array}{c}\hfill \varphi =\overline{\varphi }\left({\displaystyle \frac{1-e^{\overline{\varphi }Re(x-1)}}{1+e^{\overline{\varphi }Re(x-1)}}}\right)\mathrm{where}\overline{\varphi }\mathrm{is}\mathrm{the}\mathrm{solution}\mathrm{of}{\displaystyle \frac{\overline{\varphi }-1}{\overline{\varphi }+1}}=e^{-\overline{\varphi }Re}\end{array}$$

(98)

$\varphi \in C^{\infty }\left(\mathsf{\Omega }\right)$, i.e., $\varphi \left(x\right)$ has continuous derivatives of all orders. The differential operator $A=\varphi \frac{d}{dx}-\frac{1}{Re}\frac{d^2}{d{x}^2}$ is nonlinear, therefore only LSP yields VC integral form [17,18]. LSFEP is used to present numerical studies using the following discretizations for $k\ge 2$ with p-levels up to 17.

For $Re=100$: Element lengths (5 element mesh)

$$\begin{array}{c}\hfill 0.779,0.17,0.017,0.017,0.017\end{array}$$

For $Re={10}^6$: Element lengths (10 Element mesh)

$$\begin{array}{c}\hfill 0.4944443,0.4944443,1.0\times {10}^{-2},1.0\times {10}^{-3},5.0\times {10}^{-5},5.0\times {10}^{-5}\\ \hfill 5.0\times {10}^{-6},5.0\times {10}^{-6},1.0\times {10}^{-6},1.0\times {10}^{-7},1.0\times {10}^{-7},1.0\times {10}^{-7},1.0\times {10}^{-7}\end{array}$$

Discussion of Results

In the numerical studies, the characteristic discretization length h is fixed (due to fixed mesh), hence we can study p-convergence, k-convergence and $pk$-convergence. Figure 13a,b show plots of residual functional I versus dof for $Re=100$ and ${10}^6$ for different values of k for progressively increasing p-level up to 17.

p-convergence: For a fixed value of k, the sequence of computed finite element solutions for progressively increasing p-level represent p-convergence of the solution and the residual functional I in Figure 13a,b. For each value of k in Figure 13a,b we observe progressively increasing slope of I versus dof graphs for progressively increasing p-level, indicating progressively increasing convergence rate (shown in a priori error estimation).

k-convergence: From the I versus dof graphs for increasing p-level but for a fixed value of k, we note that with progressively increasing k, the I versus dof graphs for all values of k remain almost parallel to each other but shift downward with increasing k. This confirms that increase in k does not result in higher convergence rate but results in better accuracy (lower values of I) for a given degrees of freedom. Thus, k affects accuracy but not the convergence rate. From the curves of constant p, hence constant h and p in Figure 13a,b, thus only k varying along these curves (k-convergence), we note that increase in k results in increase of the value of functional I, indicating increasing error. This is due to the fact for fixed p, increasing k results in reduced degrees of freedom (obvious from the graphs). We note from the graphs that for a fixed dof, increase in k unconditionally results in lower values of residual functional I, hence improved accuracy.

$\mathit{pk}$-convergence: We can observe $pk$-convergence in at least two ways: (1) for a fixed k we can choose $p=2k-1$, minimum p-level for this choice of k. If we connect these points in the graphs in Figure 13a,b, we observe much higher convergence rate than p-convergence for any value of k. (2) The most dramatic and the highest convergence rates are achieved for combination of p and k for which the dofs are almost constant indicated by nearly vertical lines in the graphs in Figure 13a,b. The numerical studies for the BVP in which the differential operator is nonlinear also confirm that $h,p,k$ are three independent parameters in the finite element computations and the convergence behavior is in accordance with a priori estimates.

9.4. Poisson’s Equation (BVP)

Consider the two-dimensional steady state Poisson’s equation (see [19,20,21] also):

$$\begin{array}{cc}\hfill & \frac{\partial {}^{2}u}{\partial x^2}+\frac{\partial {}^{2}u}{\partial y^2}=-f(x,y)\forall x,y\in {\mathsf{\Omega }}_{xy}=(-a,a)\times (-b,b)\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill & u(x,-b)=u(x,b)=u(-a,y)=u(a,y)=0\hfill \end{array}$$

(100)

$f(x,y)$ is chosen that the theoretical solution of (99) and (100) is given by

$$\begin{array}{c}\hfill u(x,y)=(a^n-x^n)(b^m-y^m)\end{array}$$

(101)

$u(x,y)$ in (101) satisfies boundary conditions (100). We choose $a=b=1$, $m=n=8$. The differential operator in (99) is self-adjoint, hence GM/WF and LSP both yield variationally consistent integral forms. If $u_h$ is an approximation of u over ${\overline{\mathsf{\Omega }}}_{xy}^T={\displaystyle \bigcup _e}{\overline{\mathsf{\Omega }}}_{xy}^e$, discretization of ${\overline{\mathsf{\Omega }}}_{xy}^{}$, then

$$\begin{array}{c}\hfill {\varphi }_h\subset V\subset H^k\left({\overline{\mathsf{\Omega }}}_{xy}^T\right);k\ge 3\end{array}$$

(102)

in which $k=3$ is minimally conforming space of approximation ${\varphi }_h$ and hence ${\varphi }_h^e$ for which integrals are Riemann over $\left({\overline{\mathsf{\Omega }}}_{xy}^T\right)$. We consider two types of discretization.

Mesh A

The domain ${\overline{\mathsf{\Omega }}}_{xy}^{}$ is divided into $2\times 2$, $4\times 4,\cdots$ uniform discretizations in which elements are square and $\xi \eta$ axes for each element parallel to $xy$ and pointing in the same directions as the $xy$ axes. (Figure 14). For these discretizations we can use both HGDA/DG and HGDA/TP element formulations.

Mesh B

Initial mesh for distorted elements is shown in Figure 15. Each element of this mesh is uniformly divided into $2\times 2$, $4\times 4$, ⋯ to generate progressively refined meshes.

Computations are performed for progressively refined Mesh A and Mesh B using HGDA/DG of class $C^{11}$, $C^{22}$ and $C^{33}$ as well as using HGDA/TP elements of class $C^{11}$, $C^{22}$, $C^{33}$ for p-levels of 3, 5, 7 and 9, keeping in mind that for HGDA/TP elements only Mesh A can be used whereas for HGDA/DG both Mesh A and Mesh B can be used. The p-levels are increased starting with 3 in $\xi$ and $\eta$ directions. Figure 16, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 show graphs of residual functional I versus dof for solutions for classes $C^{11}$, $C^{22}$ and $C^{33}$ for Mesh A and Mesh B using HGDA/DG and HGDA/TP formulations. We discuss details in what follows.

Discussion of Results

• In Figure 16, Figure 17 and Figure 18, residual functional I versus dof are presented for solutions for classes $C^{11}$, $C^{22}$ and $C^{33}$ for p-levels of 3, 5, 7 and 9 for both HGDA/DG and HGDA/TP formulations using Mesh A. Slopes of I versus dof graphs for HGDA/DG, $C^{11}$, $C^{22}$, $C^{33}$ solutions for all p-levels are higher compared to the corresponding graphs for HGDA/TP solutions showing higher convergence rates of HGDA/DG elements compared to HGDA/TP.
• HGDA/DG solutions yield lower values of I for given dofs. This holds true for all three classes of solutions and for each p-level considered.
• It is significant to note that even though HGDA/TP formulations are perfectly valid for Mesh A and are believed to have optimal convergence rates, even then HGDA/DG solutions yield higher convergence rates and yield better accuracy. This needs further investigation and discussion (given later).
• Figure 19, Figure 20 and Figure 21 show plots of I versus dof for solutions for class $C^{11}$, $C^{22}$ and $C^{33}$ at p-levels of 3, 5, 7 and 9 using HGDA/DG formulation with Mesh B and HGDA/TP formulation with Mesh A. We remark that HGDA/TP element formulation is believed to be meritorious formulation for Mesh A. Obviously HGDA/TP formulation cannot be used for Mesh B as the elements of the discretizations are distorted.
• HGDA/DG formulation always has higher convergence rate and lower values of I for a given dof for solutions for classes $C^{11}$, $C^{22}$ and $C^{33}$ for each p-level considered.
• Thus, regardless of whether the discretization contains square (or rectangular) elements or distorted elements in $xy$ space, HGDA/DG formulation has higher convergence rate and better accuracy (lower values of I) compared to HGDA/TP elements for Mesh A.

9.5. 2D Convection Diffusion Equation (BVP)

We consider 2D steady state convection diffusion equation

$$\begin{array}{c}\hfill \frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}-{\displaystyle \frac{1}{Pe}}\left(\frac{\partial {}^{2}u}{\partial x^2}+\frac{\partial {}^{2}u}{\partial y^2}\right)=0\forall x,y\in {\mathsf{\Omega }}_{xy}=(0,1)\times (0,1)\end{array}$$

(103)

with boundary conditions

$$\begin{array}{cc}\hfill u(x,0)& ={\displaystyle \frac{1-e^{(x-1)Pe}}{1-e^{-Pe}}};u(0,y)={\displaystyle \frac{1-e^{(y-1)Pe}}{1-e^{-Pe}}}\hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill \mathrm{and}u(x,1)& =u(1,y)=0\hfill \end{array}$$

(105)

where $Pe$ is the Peclet number. A theoretical solution is given by

$$\begin{array}{c}\hfill u(x,y)={\displaystyle \frac{\left(1-e^{(x-1)Pe}\right)\left(1-e^{(y-1)Pe}\right)}{\left(1-e^{-Pe}\right)\left(1-e^{-Pe}\right)}}\end{array}$$

(106)

From the theoretical solution, we observe that $u(x,y)$ is analytic for all value of $Pe$ and is of class $C^{jj}$, where $j\to \infty$ is admissible.

We consider the same discretizations that are described for Mesh A and Mesh B for the Poisson’s equation and the same p-levels for solutions for classes $C^{11}$, $C^{22}$, and $C^{33}$. Figure 22, Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27 show graphs of residual functional I versus dof obtained using Mesh A and Mesh B for solutions for class $C^{11}$, $C^{22}$, $C^{33}$; p-levels of 3, 5, 7, and 9 for HGDA/DG and HGDA/TP formulations. We discuss results in the following.

• In Figure 22, Figure 23 and Figure 24, I versus dof graphs are presented for solutions for classes $C^{11}$, $C^{22}$, $C^{33}$ at p-levels of 3, 5, 7 and 9 for both HGDA/DG and HGDA/TP using Mesh A. We observe higher slopes of I versus dof graphs in case of HGDA/DG for solutions for classes $C^{11}$, $C^{22}$ and $C^{33}$ for each p-level compared to the solutions obtained from HGDA/TP.
• We note from Figure 22, Figure 23 and Figure 24 that HGDA/DG solutions yield lower values of I for a given dofs compared to HGDA/TP. This holds for solutions for classes $C^{11}$, $C^{22}$, $C^{33}$ for each p-level considered.
• Figure 25, Figure 26 and Figure 27 show graphs of I versus dof for solutions for class $C^{11}$, $C^{22}$, $C^{33}$ at p-levels of 3, 5, 7 and 9 using HGDA/DG with Mesh B and HGDA/TP with Mesh A. Obviously HGDA/TP elements cannot be used for Model B. Here also we observe a higher convergence rate of I versus dof as well as lower values of I (for given dofs) in case of HGDA/DG, hence better accuracy for given dofs. This holds for all three classes of solutions and for all p-levels considered.

9.6. 2D Burgers Equation (BVP)

Consider the 2D steady state scalar Burgers equation (see [19,20,21] also):

$$\begin{array}{c}\hfill u\frac{\partial u}{\partial x}+u\frac{\partial u}{\partial y}-\frac{1}{Re}\left(\frac{\partial {}^{2}u}{\partial x^2}+\frac{\partial {}^{2}u}{\partial y^2}\right)=f(x,y)\forall (x,y)\in {\mathsf{\Omega }}_{xy}=(-a,a)\times (-b,b)\end{array}$$

(107)

with the BCs

$$\begin{array}{c}\hfill u(-a,y)=u(a,y)=u(x,-b)=u(x,b)=0\end{array}$$

(108)

We consider $Re=10$ and choose $f(x,y)$ such that the theoretical solution of (107) and (108) is given by

$$\begin{array}{c}\hfill u(x,y)=(a^n-x^n)(b^m-y^m)\end{array}$$

(109)

In the numerical studies we consider $a=b=1$ and $m=n=8$. In this case the differential operator $A=u\frac{\partial }{\partial x}+u\frac{\partial }{\partial y}-\frac{1}{Re}\left(\frac{\partial {}^2}{\partial x^2}-\frac{\partial {}^2}{\partial y^2}\right)$ is nonlinear, hence the GM/WF will yield a VIC integral form. Modified LSP [17] yields VC integral form, hence we consider finite element formulation based on this integral form in the numerical studies. We consider the same discretizations that are described for Mesh A and Mesh B for Poisson’s equation and we consider the same p-levels for solutions of classes $C^{11}$, $C^{22}$, $C^{33}$ and the formulations based on HGDA/DG and HGDA/TP. Since, in this case the differential operator is nonlinear, the system of algebraic equations resulting from LSP are nonlinear. Their solution is obtained using Newton’s linear method with line search [17,18]. A tolerance of ${10}^{-6}$ or lower is used for assuming computed quantities to be zero. We discuss results in what follows.

(1)

Figure 28, Figure 29 and Figure 30 show plots of residual functional I versus dof for solutions for class $C^{11}$, $C^{22}$, $C^{33}$ at p-levels of 3,5,7 and 9 obtained using HGDA/DG and HGDA/TP for Mesh A. We observe higher slopes of I versus dof graph for HGDA/DG compared to HGDA/TP almost in all cases.

(2)

We also note from Figure 28, Figure 29 and Figure 30 that HGDA/DG solutions yield lower values of I compared to HGDA/TP solutions.

(3)

Figure 31, Figure 32 and Figure 33 show plots of I versus dof for solutions for class $C^{11}$, $C^{22}$ and $C^{33}$ at p-levels of 3,5,7 and 9 using HGDA/DG for Mesh B and HGDA/TP for Mesh A. Obviously, HGDA/TP formulation cannot be used for mesh B. Here also we observe higher convergence rate of I versus dof as well as lower values of I (for a given dof). This holds true for all three classes of solutions and for all three p-levels.

10. Model Problem: IVPs

In the case of IVPs, space-time coupled finite element processes with a space-time variationally consistent integral form with local approximation in $H^{\left(k\right)}\left({\overline{\mathsf{\Omega }}}_{xt}^{}\right)=H^{(k_1,k_2)}\left({\overline{\mathsf{\Omega }}}_{xt}^{}\right)$ scalar product spaces containing space-time functions of global differentiability $k_1-1$ in space and $k_2-1$ in time are most meritorious over all others [18]. Since the space-time differential operators are either non-self adjoint or non-linear, only space-time finite element processes based on space-time residual functional (space-time least squares processes) are STVC [18]. In case of initial value problems, rather than discretizing the entire space-time domain, we could consider an increment of time and discretize the space-time strip or slab for the increment of time, with only one element in time, but adequate discretization in space. Need for more elements in time can be addressed by reducing the size of the time increment. We compute converged solution (residual functional $I\to 0$) for the first increment of time. The solution for the second increment of time is calculated using second space-time strip for which initial conditions are determined using the converged solution from the first increment of time, i.e., first space-time strip.

On obtaining converged solution for second space-time strip we move on to third space-time strip and so on until entire evolution is calculated. This process is referred to as space-time strip or slab with time marching. This is the most efficient and meritorious approach of computing evolution described by the IVP. We make some remarks in the following.

Remark 7.

1.

If $2m_1$ and $2m_2$ are the highest orders of the derivatives of the dependent variables in space and time in the initial value problem, then $k_1=2m_1+1$ and $k_2=2m_2+1$ correspond to the minimally conforming scalar product space $H^{(k_1,k2)}(·)$ containing functions of space and time for which all space-time integrals will remain Riemann type over the space-time discretization.

2.

In the space-time strip with time marching method, we must ensure that the computed solution for the first space-time strip is converged before we move on to the second space-time strip as the ICs for the second space-time strip are obtained using the solution for the first space-time strip.

3.

Based on Remark 2, the convergence studies in IVPs can be done only for the first space-time strip. In reference [18], it has been shown using many model problems that the convergence characteristics and various aspects of space-time HGDA/DG and HGDA/TP element formulations as well as their performance in terms of accuracy remain the same as shown for 2D BVPs discussed in this paper. For this reason and also for the sake of brevity, numerical studies are not presented for IVPs.

4.

We emphasize that the space-time finite elements of classes $C^{qr}$, $C^{qrs}$. have precisely same formulations as the $C^{qr},C^{qrs}$ formulations for BVPs. Just like BVPs, here also HGDA/DG elements only have the required dof at the nonhierarchical nodes for desired orders of continuity and the HGDA/TP space-time elements are inferior to HGDA/DG in all aspects for the same reasons that have been discussed earlier in the context with BVPs.

11. Summary and Conclusions

In the k-version of the finite element method initiated by Surana et al. [1,2,3], the authors showed that k, the order of the approximation space defining global differentiability of approximation over a discretization is an independent parameter in all finite element processes in addition to h and p, hence the terminology k-version of the finite element method in addition to h- and p-versions is quite natural. In this paper, we considered the works of Surana et al. [1,2,3,11,12,13,14,15,16,17,18,19,20,21] on k-version finite element method and have presented a unified computational methodology incorporating the k-version that addresses all higher order global differentiability issues in BVPs and IVPs. The k-version of the finite element method presented in this paper is compared with other published works on higher order global differentiability approaches including isogeometric analysis. We present a short summary and draw some conclusions.

• The global differentiability of an approximation is due to union of local approximations, i.e., in the k-version presented in this paper, the element local approximation over finite elements are designed such that their union automatically gives the desired globally differentiability everywhere in the discretized domain. This must be intrinsic and fundamental aspect of all finite element processes considering k-version.
• The HGDA/DG are p-version and are hierarchical, i.e., in $C^q$, $C^{qr}$ and $C^{qrs}$ approximations, p-levels can be increased beyond those needed for q, $qr$, $qrs$ orders of continuity. This increase in p-level does not increase dofs at the nonhierarchical nodes that are responsible for the desired degree of global smoothness.
• It is shown and established that HGDA/DG are truly higher order global differentiability local approximations. HGDA/TP elements contain additional dofs at the nonhierarchical nodes compared to HGDA/DG elements. The consequence of this is: (a) the HGDA/TP elements have dofs at the hierarchical nodes that cannot be transformed between $\xi \eta$ and $xy$ ($\xi \eta \zeta$ and $xyz$) spaces (b) poor convergence rate compared to HGDA/DG (c) poor accuracy compared to HGDA/DG (d) HGDA/TP elements do not have true $C^q$, $C^{qr}$, $C^{qrs}$ global differentiability.
• HGDA/DG elements work perfectly well in discretizations containing distorted element geometries while HGDA/TP elements are restricted to rectangular shapes with additional restrictions on $\xi \eta$ for the elements to be parallel to $xy$ and pointing in the same directions.
• Graphs of residual functional and the error in quadratic functional are reported in the paper with respect to dof instead of using characteristic discretization length h. Since HGDA/TP elements have additional dofs at the nonhierarchical nodes compared to HGDA/DG elements, use of dof is a true measure of their performance in terms of convergence rate and the accuracy as h is not effected by additional dof in HGDA/TP elements.
• We find some published works on k-version that has been cited in the literature review section of this paper. Unfortunately, the approaches used in those publications are specific to some order of differentiability and are not p-version hierarchical and cannot be extended to orders higher than those presented in the paper. The work presented here allows any desired order of global differentiability local approximations in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ with p-version hierarchical structure. To the best of our knowledge, there is no other published work on the k-version that has these features.
• We have considered and presented considerable details of isogeometric analysis. (1) We have clearly shown that isogeometric analysis is nothing more than $C^0$ finite element discretization with higher p-level, each element being a patch in isogeometric analysis. We have clearly demonstrated that in the isogeometric approach, the global solution for the entire domain is always of class $C^0$ when the domain contains more than one patch regardless of the solution class within the patch. This is due to the fact that in isogeometric analysis, inter-patch continuity is $C^0$. (2) Isogeometric analysis is not a finite element method. A patch is like a complex domain, for which displacement approximation is constructed (no finite elements yet). The patch is subdivided into smaller patches that are mapped in two unit squares or cubes for the sake of integration over the patch. Smaller patches or subdomains of the patch have no local approximation. Each smaller patch or subdomain has to inherit the displacement approximation constructed for the whole patch (no finite element yet either).
• In isogeometric analysis, benefits of higher order smoothness of geometry have never been shown. Our work on the k-version shows that the benefit of accuracy and convergence rates are due to local approximation and not due to higher order geometry.
• A comparison of isogeometric method with k-version [1,2,3,19,20,21] is meaningless for two reasons (i) isogeometric analysis is not a finite element method (ii) the global approximation in isogeometric is only $C^0$, so where is the k-version?

In conclusion this paper presents a true k-version finite element method (HGDA/DG) in which local approximation can be of class $C^q$, $C^{qr}$, $C^{qrs}$p-version hierarchical that always yield global differentiability of orders q, $qr$ and $qrs$ in ${\mathbb{R}}^1$, ${\mathbb{R}}^2$, ${\mathbb{R}}^3$. There is no other published works on the k-version of finite element method or higher order global differentiability approximation that has these features.