Weakening the tight coupling between geometry and simulation in isogeometric analysis: from sub- and super- geometric analysis to Geometry Independent Field approximaTion (GIFT)

TL;DR

This paper introduces a generalized isogeometric analysis framework called GIFT that decouples geometry parameterization from solution approximation, enabling flexible, adaptive local refinement while maintaining exact CAD geometry.

Contribution

It proposes a novel approach allowing different spaces for geometry and solution approximation, enhancing flexibility and local refinement without re-parameterizing the domain.

Findings

Achieves optimal convergence rates with non-nested spaces.

Maintains exact CAD geometry during analysis.

Enables adaptive local refinement without weakening CAD link.

Abstract

This paper presents an approach to generalize the concept of isogeometric analysis (IGA) by allowing different spaces for parameterization of the computational domain and for approximation of the solution field. The method inherits the main advantage of isogeometric analysis, i.e. preserves the original, exact CAD geometry (for example, given by NURBS), but allows pairing it with an approximation space which is more suitable/flexible for analysis, for example, T-splines, LR-splines, (truncated) hierarchical B-splines, and PHT-splines. This generalization offers the advantage of adaptive local refinement without the need to re-parameterize the domain, and therefore without weakening the link with the CAD model. We demonstrate the use of the method with different choices of the geometry and field splines, and show that, despite the failure of the standard patch test, the optimum convergence rate is achieved for non-nested spaces.

Tables

Table 1. Table 1 : Control points and weights for coarsest parameterization Q 0 subscript 𝑄 0 Q_{0}

$(i,j)$	$P_{ij}^{Q_0,x}$	$P_{ij}^{Q_0,y}$	$w_{ij}^{Q_0}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	1.0000	0.7071
(1,3)	0.0000	1.0000	1.0000
(2,1)	2.0000	0.0000	1.0000
(2,2)	2.0000	2.0000	0.7071
(2,3)	0.0000	2.0000	1.0000

Table 2. Table 2 : Various parameterizations

	$p_{\xi }$	$p_{\eta }$	$s$	$t$	$P_{ij}$ & $w_{ij}$
$A_1$	1	2	2/3	1/8	Table 3
$B_1$	1	2	0.17	0.81	Table 4
$C_1$	1	2	2/3	1/8	Table 5
$D_1$	1	2	2/3	1/8	Table 6

Table 3. Table 3 : Control points and weights for parameterization A 1 subscript 𝐴 1 A_{1}

$(i,j)$	$P_{ij}^{A_1,x}$	$P_{ij}^{A_1,y}$	$w_{ij}^{A_1}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	0.0917	0.9634
(1,3)	0.8320	1.0000	0.7437
(1,4)	0.0000	1.0000	1.0000
(2,1)	1.6667	0.0000	1.0000
(2,2)	1.6667	0.1529	0.9634
(2,3)	1.3865	1.6667	0.7437
(2,4)	0.0000	1.6667	1.0000
(3,1)	2.0000	0.0000	1.0000
(3,2)	2.0000	0.1835	0.9634
(3,3)	1.6639	2.0000	0.7437
(3,4)	0.0000	2.0000	1.0000

Table 4. Table 4 : Control points and weights for parameterization B 1 subscript 𝐵 1 B_{1}

$(i,j)$	$P_{ij}^{B_1,x}$	$P_{ij}^{B_1,y}$	$w_{ij}^{B_1}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	0.7509	0.7628
(1,3)	0.1423	1.0000	0.9444
(1,4)	0.0000	1.0000	1.0000
(2,1)	1.1700	0.0000	1.0000
(2,2)	1.1700	0.8786	0.7628
(2,3)	0.1665	1.1700	0.9444
(2,4)	0.0000	1.1700	1.0000
(3,1)	2.0000	0.0000	1.0000
(3,2)	2.0000	1.5018	0.7628
(3,3)	0.2845	2.0000	0.9444
(3,4)	0.0000	2.0000	1.0000

Table 5. Table 5 : Control points and weights for parameterization C 1 subscript 𝐶 1 C_{1}

$(i,j)$	$P_{ij}^{C_1,x}$	$P_{ij}^{C_1,y}$	$w_{ij}^{C_1}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	0.0917	0.9634
(1,3)	0.8320	1.0000	0.7437
(1,4)	0.0000	1.0000	1.0000
(2,1)	1.6667	0.0000	1.0000
(2,2)	1.1000	0.3000	0.8000
(2,3)	0.7500	1.4000	0.7500
(2,4)	0.0000	1.6667	1.0000
(3,1)	2.0000	0.0000	1.0000
(3,2)	2.0000	0.1835	0.9634
(3,3)	1.6639	2.0000	0.7437
(3,4)	0.0000	2.0000	1.0000

Table 6. Table 6 : Control points and weights for parameterization D 1 subscript 𝐷 1 D_{1}

$(i,j)$	$P_{ij}^{D_1,x}$	$P_{ij}^{D_1,y}$	$w_{ij}^{D_1}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	0.0917	1.0000
(1,3)	0.8320	1.0000	1.0000
(1,4)	0.0000	1.0000	1.0000
(2,1)	1.6667	0.0000	1.0000
(2,2)	1.6667	0.1529	1.0000
(2,3)	1.3865	1.6667	1.0000
(2,4)	0.0000	1.6667	1.0000
(3,1)	2.0000	0.0000	1.0000
(3,2)	2.0000	0.1835	1.0000
(3,3)	1.6639	2.0000	1.0000
(3,4)	0.0000	2.0000	1.0000

Table 7. Table 7 : Results of various patch tests

$T_{G_i,S_j}^{\ell }$	Laplace Eq. (20)	Elasticity Eq. (22)
$T_{Q_0,A_1}^0$	1.3815e-15	3.0871e-14
$T_{Q_0,A_2}^0$	5.2147e-15	1.7986e-14
$T_{Q_0,C_1}^0$	0.0182	0.0050
$T_{Q_0,C_2}^0$	0.0023	0.0012
$T_{A_1,A_1}^1$	1.0023e-15	1.1675e-14
$T_{A_1,A_2}^1$	4.3958e-14	1.2547e-14
$T_{A_2,A_1}^1$	1.4059e-15	1.5525e-15
$T_{B_1,A_1}^2$	1.4755e-15	2.9941e-15
$T_{B_1,A_2}^2$	2.1639e-15	1.2118e-14
$T_{B_2,A_1}^2$	1.0144e-15	5.4590e-15
$T_{C_1,C_1}^3$	1.1061e-15	1.6439e-14
$T_{C_1,C_2}^3$	1.8263e-15	2.8737e-15
$T_{C_2,C_1}^3$	1.2062e-15	5.6517e-14
$T_{C_1,A_1}^4$	0.0203	0.0085
$T_{C_1,A_2}^4$	0.0016	0.0009
$T_{C_2,A_1}^4$	0.0203	0.0085
$T_{A_1,D_1}^5$	0.0188	0.0214
$T_{A_1,D_2}^5$	0.0121	0.0039
$T_{A_1,D_0}^5$	0.5418	0.1411

Table 8. Table 8 : Control points in the parametrization of Example  4.2.2 .

$(i,j)$	$P_{ij}^x$	$P_{ij}^y$	$w_{ij}$
(1,1)	1.0000	0.0000	1.0000
(1,2)	1.0000	0.4142	0.8536
(1,3)	0.4142	1.0000	0.8536
(1,4)	0.0000	1.0000	1.0000
(2,1)	2.5000	0.0000	1.0000
(2,2)	2.5000	1.5000	0.8000
(2,3)	1.5000	2.5000	0.8000
(2,4)	0.0000	2.5000	1.0000
(3,1)	4.0000	0.0000	1.0000
(3,2)	4.0000	4.0000	1.0000
(3,3)	4.0000	4.0000	1.0000
(3,4)	0.0000	4.0000	1.0000

Table 9. Table 9 : Control points for the parametrization of one-eighth of a sphere.

$(i,j,k)$	$P_{ij}^x$	$P_{ij}^y$	$P_{ij}^z$	$w_{ij}$
(1,1,1)	1.0000	0.0000	0.0000	1.0000
(1,2,1)	1.0000	1.0000	0.0000	0.7071
(1,3,1)	0.0000	1.0000	0.0000	1.0000
(2,1,1)	2.0000	0.0000	0.0000	1.0000
(2,2,1)	2.0000	2.0000	0.0000	0.7071
(2,3,1)	0.0000	2.0000	0.0000	1.0000
(1,1,2)	1.0000	0.0000	1.0000	0.7071
(1,2,2)	1.0000	1.0000	1.0000	0.5000
(1,3,2)	0.0000	1.0000	0.0000	0.7071
(2,1,2)	2.0000	0.0000	2.0000	0.7071
(2,2,2)	2.0000	2.0000	2.0000	0.5000
(2,3,2)	0.0000	2.0000	0.0000	0.7071
(1,1,3)	0.0000	0.0000	1.0000	1.0000
(1,2,3)	0.0000	1.0000	1.0000	0.7071
(1,3,3)	0.0000	1.0000	0.0000	1.0000
(2,1,3)	0.0000	0.0000	2.0000	1.0000
(2,2,3)	0.0000	2.0000	2.0000	0.7071
(2,3,3)	0.0000	2.0000	0.0000	1.0000

Table 10. Table 10 : Summary of patch tests and optimal convergence. See Section  3.2 for the notations Q 𝑄 Q , A 𝐴 A , B 𝐵 B , C 𝐶 C , and D 𝐷 D . For the degree parity, we used the naming convention with respect to the geometry, see Table  11 .

Geometry	Solution	Degree	Patch	Optimal
parameterization	basis	parity	test	convergence
$Q_0$	$A_1$	Iso-geometric	✓	✓
$Q_0$	$A_2$	Super-geometric	✓	✓
$Q_0$	$C_1$	Iso-geometric	$\times$	✓
$Q_0$	$C_2$	Super-geometric	$\times$	✓
$A_1$	$A_1$	Iso-geometric	✓	✓
$A_1$	$A_2$	Super-geometric	✓	✓
$A_2$	$A_1$	Sub-geometric	✓	✓
$B_1$	$A_1$	Iso-geometric	✓	✓
$B_1$	$A_2$	Super-geometric	✓	✓
$B_2$	$A_1$	Sub-geometric	✓	✓
$C_1$	$C_1$	Iso-geometric	✓	✓
$C_1$	$C_2$	Super-geometric	✓	✓
$C_2$	$C_1$	Sub-geometric	✓	✓
$C_1$	$A_1$	Iso-geometric	$\times$	✓
$C_1$	$A_2$	Super-geometric	$\times$	✓
$C_2$	$A_1$	Sub-geometric	$\times$	✓
$A_1$	$D_1$	Iso-geometric	$\times$	✓
$A_1$	$D_2$	Super-geometric	$\times$	✓
$A_1$	$D_0$	Sub-geometric	$\times$	✓

Table 11. Table 11 : Naming convention

	$p_u=p_g$		$p_u>p_g$
FEM	Iso-parametric	Super-parametric	Sub-parametric
GIFT	Iso-geometric	Sub-geometric	Super-geometric

Equations

$$F:\mathcal{P}\to \Omega ,\mathbf{x}=F(\mathbf{\xi }),\mathbf{x}\in \Omega ,\mathbf{\xi }\in \mathcal{P}.$$

$$F(\mathbf{\xi })=\sum _{i_1=1}^{n_1}\sum _{i_2=1}^{n_2}\dots \sum _{i_d=1}^{n_d}{\mathbf{C}}_{i_1,i_2,\dots ,i_d}{N}_{i_1,i_2,\dots ,i_d}(\mathbf{\xi }),$$

$$\displaystyle\boldsymbol{I}=\big{\{}(i_{1},i_{2},\ldots,i_{d}):i_{1}\in\{1,\ldots,n_{1}\},\ldots,i_{d}\in\{1,\ldots,n_{d}\}\big{\}}$$

$$\displaystyle\boldsymbol{J}=\big{\{}(i_{1},i_{2},\ldots,i_{d}):i_{1}\in\{1,\ldots,m_{1}\},\ldots,i_{d}\in\{1,\ldots,m_{d}\}\big{\}}.$$

$$\mathbf{x}(\mathbf{\xi })=\sum _{\mathbf{k}\in \mathbf{I}}{\mathbf{C}}_{\mathbf{k}}{N}_{\mathbf{k}}(\mathbf{\xi }).$$

$$J_{ij}(\mathbf{\xi })=\frac{\partial x_i}{\partial {\xi }_j}(\mathbf{\xi })=\sum _{\mathbf{k}\in \mathbf{I}}{\mathbf{C}}_{\mathbf{k}i}\frac{\partial N_{\mathbf{k}}(\mathbf{\xi })}{\partial {\xi }_j}.$$

$$u(\mathbf{\xi })=\sum _{\mathbf{k}\in \mathbf{J}}{u}_{\mathbf{k}}{M}_{\mathbf{k}}(\mathbf{\xi }),$$

$$\left(\begin{array}[]{c}\dfrac{\partial M_{\boldsymbol{k}}(\xi,\eta)}{\partial x}\\ \dfrac{\partial M_{\boldsymbol{k}}(\xi,\eta)}{\partial y}\end{array}\right)=\left(\begin{array}[]{cc}\dfrac{\partial\xi}{\partial x}&\dfrac{\partial\eta}{\partial x}\\ \dfrac{\partial\xi}{\partial y}&\dfrac{\partial\eta}{\partial y}\end{array}\right)\left(\begin{array}[]{c}\dfrac{\partial M_{\boldsymbol{k}}(\xi,\eta)}{\partial\xi}\\ \dfrac{\partial M_{\boldsymbol{k}}(\xi,\eta)}{\partial\eta}\end{array}\right).$$

$${\nabla }_{\mathbf{x}}{M}_{\mathbf{k}}(\mathbf{\xi })=J^{-T}(\mathbf{\xi }){\nabla }_{\mathbf{\xi }}{M}_{\mathbf{k}}(\mathbf{\xi }),$$

$$a(u,v)=\ell (v),$$

$$\mathbf{K}\mathbf{u}=\mathbf{f},$$

$${\mathbf{K}}_{ij}=a(M_i(\mathbf{x}),M_j(\mathbf{x})),{\mathbf{f}}_i=\ell (M_i(\mathbf{x})),$$

$$a(u,v)={\int }_{\Omega }\nabla u\cdot \nabla vd\Omega ,\ell (v)={\int }_{\Omega }fvd\Omega +{\int }_{\partial {\Omega }_N}gvd\Gamma ,$$

$$\begin{array}{rl}{\mathbf{K}}_{ij}& ={\int }_{\Omega }{\nabla }_{\mathbf{x}}{M}_i(\mathbf{x})\cdot {\nabla }_{\mathbf{x}}{M}_j(\mathbf{x})d\Omega ,\\ {\mathbf{f}}_i& ={\int }_{\Omega }f(\mathbf{x})M_i(\mathbf{x})d\Omega +{\int }_{\partial {\Omega }_N}g(\mathbf{x})M_i(\mathbf{x})d\Gamma .\end{array}$$

$$a(u,v)={\int }_{\Omega }\mathbf{ϵ}(u{)}^T\mathbf{D}\mathbf{ϵ}(v)d\Omega ,\ell (v)={\int }_{\Omega }fvd\Omega +{\int }_{\partial {\Omega }_N}\tilde{t}vd\Gamma ,$$

$$\boldsymbol{D}=\dfrac{E}{(1+\nu)(1-2\nu)}\left(\begin{array}[]{ccc}1-\nu&\nu&0\\ \nu&1-\nu&0\\ 0&0&(1-2\nu)/2\end{array}\right),$$

$$\boldsymbol{\epsilon}(u)=\left(\begin{array}[]{c}\partial u_{x}/\partial{x}\\ \partial u_{y}/\partial{y}\\ \partial u_{x}/\partial{y}+\partial u_{y}/\partial{x}\end{array}\right)$$

$$\begin{array}{rl}{\mathbf{K}}_{ij}& ={\int }_{\Omega }{\mathbf{B}}_i^T\mathbf{D}{\mathbf{B}}_{j}d\Omega ,\\ {\mathbf{f}}_i& ={\int }_{\Omega }f(\mathbf{x})M_i(\mathbf{x})d\Omega +{\int }_{\partial {\Omega }_N}\tilde{t}(\mathbf{x})M_i(\mathbf{x})d\Gamma ,\end{array}$$

$$\boldsymbol{B}_{i}=\left(\begin{array}[]{cc}\dfrac{\partial M_{i}(\xi,\eta)}{\partial{x}}&0\\ 0&\dfrac{\partial M_{i}(\xi,\eta)}{\partial{y}}\\ \dfrac{\partial M_{i}(\xi,\eta)}{\partial{y}}&\dfrac{\partial M_{i}(\xi,\eta)}{\partial{x}}\\ \end{array}\right),$$

$$u(x,y)=1+x+y,$$

$$\begin{array}{rl}\Delta u& =0,\text{in}\Omega \\ u{|}_{\partial \Omega }(x,y)& =1+x+y.\end{array}$$

$$\begin{array}{rl}{\sigma }_{ij,j}& =0,\text{in}\Omega \\ t_i& ={\sigma }_0{n}_i,\text{at}r=1,2,\\ u_2& =0,t_1=0,\text{at}\theta =0,\\ u_1& =0,t_2=0,\text{at}\theta =\pi /2,\end{array}$$

$$u_1(x)=\frac{(1+\nu )(1-2\nu )}{E}{\sigma }_{0}x,u_2(y)=\frac{(1+\nu )(1-2\nu )}{E}{\sigma }_{0}y.$$

$${\Sigma }_0=\{0,0,1,1\}\text{and}{\Pi }_0=\{0,0,0,1,1,1\},$$

$${\Sigma }_1=\{0,0,s,1,1\}\text{and}{\Pi }_1=\{0,0,0,t,1,1,1\}.$$

$$\begin{array}{rl}\left(\begin{array}{c}x\\ y\end{array}\right)& =\sum _{ij}\left(\begin{array}{c}{P}_{ij}^{A_1,x}\\ P_{ij}^{A_1,y}\end{array}\right)\times \{\text{Basis of }{A}_1\}\end{array}$$

$$u=\sum _{ij}{u}_{ij}\times \{\text{Basis of }{A}_1\}$$

$$u=1+x+y,$$

$$u_{ij}=1+P_{ij}^{A_1,x}+P_{ij}^{A_1,y}.$$

$$\begin{array}{rl}\left(\begin{array}{c}x\\ y\end{array}\right)& =\sum _{ij}\left(\begin{array}{c}{P}_{ij}^{A_1,x}\\ P_{ij}^{A_1,y}\end{array}\right)\times \{\text{Basis of }{A}_1\}\\ & =\sum _{ij}\left(\begin{array}{c}{P}_{ij}^{A_2,x}\\ P_{ij}^{A_2,y}\end{array}\right)\times \{\text{Basis of }{A}_2\}.\end{array}$$

$$u=\sum _{ij}{u}_{ij}\times \{\text{Basis of }{A}_2\}$$

Full text

1 Introduction

1.1 Contribution

We present an approach which enables the use of separate approximation spaces for the field and geometry in isogeometric finite element methods (IGAFEM). For example, coarse NURBS (non-uniform rational B-splines) approximations can be used for the geometry, and locally adapted PHT splines can be used for the field variables. This endows the method with the flexibility to locally enrich the approximation space for the field variables without modifying the spline space used for parametrization of the geometry of the domain. We verify the approach with various parameterizations of the geometry and approximations of the field variables, and on two and three dimensional geometries. We give a detailed mathematical explanation on the ability of the resulting method to pass the patch test, and to converge with optimal rate. Geometry Independent Field approximaTion (GIFT), as we coined the approach, represents a simple and computationally efficient alternative to IGA, which preserves tight integration with CAD, and allows adaptive local refinement of the solution field.

1.2 Background

Isogeometric analysis (IGA) was introduced in [25] to establish a direct link between the computer-aided-design (CAD) and analysis. Over the last decade, fuelled by rapid developments in computer graphics and CAD, and due to a number of advantages offered by spline basis functions over standard (Lagrange) finite element analysis (FEA), such as higher continuity, exact representation of the geometry, simplified integration and their behaviour in dynamics, IGA has created a large amount of interest in computational engineering and science.

The method has found many applications in various areas such as structural vibrations [48], fluids-structure interaction [2], shell analysis [11], and fracture mechanics [18]. Although the finite element version of IGA suffers from significant difficulties associated with generating volumetric parameterisations directly from CAD, coupling IGA with boundary element methods enables to perform computations directly from the CAD description of the boundary of the domain [3]. Applications include stress analysis [58, 32, 51], shape optimization [54, 31, 30], fracture mechanics [44, 45, 43, 46], geomechanics [6, 33, 4, 8, 7, 34, 69, 5], acoustics [59, 27], and electromagnetics [13]. A detailed overview of the recent work in the field can be found in [39] and in [29]. Other approaches, which also aim at facilitating the transition between CAD and analysis have been developed since the late 1990’s, see e.g., [26, 21], and include (in order of appearance):

Implicit/immersed boundary definitions

with the extended finite element method and level sets [10], [35] and [36, 37], which allows an implicit definition of arbitrary solids through constructive geometry and different level set functions. A promising approach, known as the finite cell method was proopsed in 2007 by [42].

Subdivision surfaces

[15, 14], where the idea is to use subdivision surfaces as a common representation for geometric modeling and numerical simulation in a unified framework. Subdivision surfaces is a flexible and efficient tool for modeling surfaces with arbitrary topology, avoiding many problems inherent in traditional geometry modeling approach with tensor-product patches. The corresponding solvers are highly scalable, and provide an efficient computational tool for numerical simulation and design optimization.

NURBS-enhanced finite element methods

[55, 57, 56], where the idea is to construct new finite elements with at least one edge represented as a NURBS or B-spline. This, therefore, enables to use the CAD geometry directly, but to retain the standard finite element method everywhere except in the vicinity of the boundary.

1.3 NURBS: Advantages and limitations

NURBS are the most common type of splines used in the CAD industry to describe geometry. NURBS shape functions are defined by means of B-splines (piece-wise polynomial functions defined over a knot vector) and a set of control points with associated weights. The main feature of NURBS is the ability to exactly represent conic sections, e.g. ellipses, parabolas, hyperbolas, etc., which are widely used in computer-aided-design. NURBS basis functions also posses all the properties of the standard FEM basis functions, such as compact support, linear independence, and partition of unity. Moreover, the high order continuity ($C^{p-k}$, with $p$ the polynomial order and $k$ the knot repetition) of NURBS basis functions across elements facilitates the solution of partial differential equations (PDEs) of arbitrary high order. However, NURBS are based on tensor-product structure, which does not facilitate local refinement (see Fig. 2(a)), NURBS suffer from inability to produce the watertight geometries for general shapes, which imposes difficulties for mesh generation and requires coupling algorithms for multi-patch geometries [40, 16, 20, 66, 67, 64].

1.4 Motivation

The main idea of IGA is to use the same shape functions (splines) for both: parameterization of the geometry (computational domain) and approximation of the unknown solution, see Fig. 1. This is a clear advantage in engineering design for example, as any modification of the CAD geometry is directly inherited by the approximation of the unknown fields so that the mesh need not be regenerated at each iteration of the geometrical design process. Because of the nature of the approximations most commonly used in geometrical design (NURBS), which are based on tensor-product constructions, the potential power of IGA is not always realized in practice. Indeed, local refinement is not natural in NURBS approximations. Moreover, situations usually arise when the approximation used for the field variables must be refined locally, whilst keeping the same (coarser) geometrical approximation. This is the case in shape optimization [54, 31, 30], or for problems with singularities, boundary layers or steep gradients.

It appears, therefore, that it is useful to develop more general methods where the approximation space used for the geometry and that used for the field variables be decoupled. The investigation of the behaviour of such methods for linear elasticity and Laplace equation is the focus of this paper.

We base our work on the developments made by the team of Gernot Beer in two important papers for isogeometric boundary element methods: [7, 33], which were further developed in a series of papers targeting geomechanics problems [6, 4, 8, 7, 34, 69, 5].

The most important short-comings of IGA, which served as motivation to investigate the possibility to decouple the geometry and the solution fields, can be summarized as follows:

The $h$- and $p$-refinement of the original exact geometry is redundant in IGA. The computational savings on non-refining the geometry were estimated in [7], and a speed up factor of up to 2 was achieved in the context of BEM. It is expected that for large industrial problems, the savings will be even more significant. 2. 2.

NURBS, the de facto standard in CAD community, do not offer local refinement. Therefore, in order to use any other locally-refined splines with the original CAD-model given by NURBS, staying within the iso-parametric concept will require re-parameterization of the original geometry. Note that some locally-refined splines (e.g., T-splines) allow exact representation of the original geometry, but others (e.g., polynomial based splines) do not offer this advantage. Therefore, such operations (particularly, for efficient local refinement) are not only time-consuming, but may also lead to loosing the geometry exactness, and may introduce the need to communicate with the CAD model at every step of the solution refinement process.

1.5 Alternatives to NURBS

In many applications, the efficiency of a numerical method can be drastically improved by restricting the refinement to certain areas (for example, where the solution exhibits high gradients). Therefore, the local refinement issue in IGA is currently an active research topic. A recent survey on locally refinable splines is given in [28]. In the authors’ opinion, the most commonly used alternatives which allow local refinement, are THB-splines, T-splines and PHT-splines (polynomial splines over hierarchical T-meshes), which are briefly discussed below.

1.5.1 Hierarchical B-splines and Truncated Hierarchical B-splines (THB-splines)

Hierarchical B-splines were first introduced in [22] for surface fitting. They posses such properties as maximum continuity on each refinement level, and linear independence of the basis functions. Later, in [62], these shape functions were modified to reduce the local support and assure the partition of unity property. These splines were called Truncated Hierarchical B- (THB-) splines. The theoretical background, as well as applications to isogeometrical analysis, can be found in [24] and [23].

1.5.2 T-splines

T-splines, introduced by Sederberg et al. in [53, 52], are defined over so-called T-meshes, see Fig. 2(b). T-splines are piecewise rational functions, which preserve the exactness of NURBS geometry. Moreover, a multi-patch NURBS geometry parameterization can be converted into a single patch T-splines description without gaps and leaks. However, in the most general case, the linear independence of T-splines blending functions is not guaranteed, and this led to the introduction of so-called analysis suitable T-splines [9, 12]. The analysis suitable T-splines have been successfully applied to IGA, see e.g. [17]. In a recent paper [70], the degree elevation of T-splines have also been studied. However, as described in [50], the implementation, in particular the algorithm of knot insertion which preserves analysis-suitability, is complicated.

1.5.3 PHT-splines

PHT-splines were introduced by Deng et al. in 2006 [19]. In addition to main properties of B-splines, the main advantage of PHT-splines, which makes them attractive for IGA, is efficient and simple refinement algorithm [41]. However, the trade off is the reduced continuity ($C^{1}$). Nevertheless, this is enough for most applications in solid and structural mechanics. Since PHT-splines are polynomials, the CAD geometry of arbitrary topology may not be represented exactly, which led to the development of rational PHT-splines [63].

1.5.4 Geometric design and IGA considerations for unstructured quadrilateral meshes

Recently, Toshniwal et al. [61] proposed a new framework for geometric design and IGA on unstructured quadrilateral meshes $M$. Following the construction of D-patch framework developed in [49], and the construction of $C^{1}$ smooth spline elements in the neighborhood of extra-ordinary points (albeit only for PHT-splines) in [38], the focus of [61] is on the construction of smooth, linearly independent, locally supported spline functions over $M$. To achieve the smoothness at the extra-ordinary points, affine-invariant linear transformations (called smoothing matrices) from [49] are employed. As oppose to [38], the work in [61] uses the smoothing matrices with non-negative coefficients. Although it yields higher quality surfaces for geometric modeling applications, this comes at the cost of slight (very small) variations in geometries during refinement. The design and analysis spline spaces are then simply defined as the span of these splines. The design space $S_{D}$ is contained within the analysis space $S_{A}$ at $k=0$ (the coarsest level for analysis), and under the assumption of idempotence of the smoothing matrices, the spaces $S_{A}^{k}$ are nested, i.e. $S_{D}\subseteq S_{A}=:S_{A}^{0}\subseteq S_{A}^{1}\subseteq S_{A}^{2}\ldots\subseteq S_{A}^{k}$. For isogeometric compatibility, simple transformation rules for change of basis from $S_{D}$ to $S_{A}$ are constructed.

1.6 Generalizing the isogeometric concept

In this paper, instead of aiming to construct tailored spline spaces for dealing with local refinement and/or unstructured meshes, we aim to generalize the isogeometric concept, namely, to possibly choose different spaces for design and analysis (where $S_{D}$ may not be contained in $S_{A}$). Thereby, while remaining in the realm of existing technologies for design, we can use suitable spline spaces for analysis.

In order to preserve the isogeometric concept, with the three choices of splines described above, it is required to convert the original NURBS CAD geometry to a parametrization (exact or approximate) in the same spline space. This contradicts the original objective of the IGA to bring the direct link between the CAD design and analysis. This limitation, caused by tight integration of the geometry parameterization and the approximation of the solution (Fig. 1) in IGA, motivated the authors to formulate a generalized approach: Geometry Independent Field approximaTion (GIFT ​). The main idea of GIFT ​is to retain the original CAD geometry whilst flexibly adapting the solution basis to best fit the solution field (Fig. 3). For example, a NURBS geometry can be used together with a PHT-splines approximation for the solution. The main features of GIFT ​can be summarized as follows:

• •

Preserve exact CAD geometry (provided in any form, including B-splines or NURBS) at any stage of the solution process;

• •

Allow local refinement of the solution by choosing appropriate field approximations, as independently as possible of the geometrical parameterization of the domain, including partition of unity enrichment, and

• •

Exploit computational savings by not refining the geometry during the solution refinement process, and by choosing simpler approximations for the solution, i.e. polynomial functions instead of rational functions.

This paper aims at showing a proof of concept for GIFT ​. We present the general framework of the method, followed by a number of numerical examples. In our implementation of GIFT ​, the geometry is given by NURBS because this is the most commonly used form of geometry description in CAD. For solution approximations, we use a variety of approximation schemes, including NURBS functions (which have potentially different degrees and knot vectors compared to those used for the geometry), PHT-splines and B-splines.

Remark

In the engineering community, such different choices of spline spaces for geometry and field give rise to the question of verification of patch-test. The role of the patch test in the convergence of the finite-element based method is a long standing debate. It has been shown in [60], that for the method to converge with optimal rate, the patch test is neither sufficient nor necessary. The present study also investigates the relation of the convergence properties of the method with the standard patch test. The spline spaces for the geometry and the field are chosen such that they may, or may not, pass the patch test.

The organization of the remainder of this article is as follows. We first present a simple introduction to the mathematical formulation of GIFT ​​ in Section 2. In Section 3, we present the study of various patch tests. The numerical examples on convergence studies are presented in Section 4. These patch tests and numerical examples are carefully designed to study various combinations of bases for geometry and solution approximation. Finally, some conclusions and recommendations for future work are outlined in Section 6.

2 Formulation of GIFT ​

We consider an open domain $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$, with boundary $\Gamma$ consisting of two parts $\Gamma_{D}$ and $\Gamma_{N}$, such that: $\Gamma=\overline{\Gamma_{D}\bigcup\Gamma_{N}}$, $\Gamma_{D}\bigcap\Gamma_{N}=\emptyset$. The domain $\Omega$ is parameterized on a parametric domain $\mathcal{P}$ by mapping $F$:

[TABLE]

In what follows, we will denote the NURBS basis by $\mathcal{N}_{i,j}$, and the B-splines basis by $\mathcal{B}_{i,j}$, with $i$ denoting the degree in the $\xi$ direction, and $j$ denoting the degree in the $\eta$ direction. Typically, the geometrical map $F$ is given by a set of basis functions $N_{i_{1},i_{2},\ldots,i_{d}}(\boldsymbol{\xi})$ and a set of control points $\boldsymbol{C}_{i_{1},i_{2},\ldots,i_{d}}$ as

[TABLE]

where $N_{i_{1},i_{2},\ldots,i_{d}}(\boldsymbol{\xi})$ can be a tensor product of NURBS, B-splines, T-splines, PHT-splines, etc. For brevity reasons, we introduce two sets of multi-indices $(i_{1},i_{2},\ldots,i_{d})$ of NURBS basis functions by

[TABLE]

Moreover, wherever suitable, for multi-index $(i_{1},i_{2},\ldots,i_{d})$ we will interchangeably use the collapsed notation $\boldsymbol{k}$. Thence, Eq. (1) and (2) are written as

[TABLE]

In what follows, we will refer to the set $\{N_{\boldsymbol{k}}(\boldsymbol{\xi})\}_{\boldsymbol{k}\in\boldsymbol{I}}$ as the geometry basis. For change of variables, we will also need the Jacobian matrix $J(\boldsymbol{\xi})$ of the mapping $F$, which is given by

[TABLE]

The departure from classical isogeometric analysis consists in choosing a solution basis $\{M_{\boldsymbol{k}}(\boldsymbol{\xi})\}_{\boldsymbol{k}\in\boldsymbol{J}}$, which is different from the geometry basis, and looking for the solution in the form:

[TABLE]

where $u_{\boldsymbol{k}}$ are the unknown control variables. In order to evaluate derivatives of the solution basis function $M_{\boldsymbol{k}}(\boldsymbol{\xi})$ with respect to variables $\boldsymbol{x}$, the standard chain rule is used, which in two dimensions read:

[TABLE]

In matrix notations (independent of dimensions), this can be written as follows

[TABLE]

where $J^{-T}(\boldsymbol{\xi})$ is the transpose of the inverse of the Jacobian matrix (5).

Now let the weak form of the boundary value problem be given by

[TABLE]

then substituting Eq. (4) and (6) into (9) we obtain the linear system of equations

[TABLE]

where the stiffness matrix $\boldsymbol{K}$ and the force vector $\boldsymbol{f}$ are given by

[TABLE]

and vector $\boldsymbol{u}$ consists of all unknown control variables $u_{j}$.

For our numerical results, we consider two problems, namely, a) Poisson’s equation, and b) the problem of linear elasticity. In Sections 2.1-2.2, we briefly describe the derivations of the weak form of these problems.

2.1 Poisson’s equation

For the Poisson’s problem, the terms in the weak form are defined as

[TABLE]

where $f$ is the source function, and $g$ is the flux prescribed on Neumann part of the boundary. The Eq. (11) thus results in, for $i,j=1,\ldots,M$,

[TABLE]

2.2 Linear elasticity

Consider the linear elasticity problem in two dimensions, for which the weak form is given by

[TABLE]

where $u:=(u_{x},u_{y})^{T}$ denotes the vector of unknown displacements, $f$ denotes the body forces, $\tilde{t}$ denotes the tractions prescribed on Neumann boundary, and the matrix of material parameters in plane strain is given by

[TABLE]

where $E$ is the Young’s modulus and $\nu$ is the Poisson’s coefficient. The strain operator, which is defined as

[TABLE]

applied to the solution approximation given by Eq. (6) yields the following stiffness matrix and force vector:

[TABLE]

where the strain operator $\boldsymbol{B}_{i}$ is defined as

[TABLE]

and derivatives of the solution shape functions with respect to geometry variables are evaluated according to Eq. (7).

3 Patch tests

In this Section, we will study the classical patch test for various combinations of bases for the geometry representation and solution approximation. For all the test cases, we will study two problems, namely, Laplace equation, and linear elasticity problem, which are defined in Section 3.1. We present "engineering" analysis of the patch test results (including the values of the field control variables). This can also serve as a guide for engineers to design patch tests in spline based methods.

3.1 Problem definition

We formulate two problems in a quarter annulus domain $\Omega$, as shown in Fig. 4. The first problem is for the two-dimensional Laplace equation, where we use the linear solution

[TABLE]

imposed as the Dirichlet boundary condition on the entire boundary $\partial\Omega$, i.e.

[TABLE]

The second problem is for the linear elasticity, which is given by:

[TABLE]

which admits a linear solution (plane-strain case) given by

[TABLE]

3.2 Geometry parametrization

As shown in Fig. 5, the coarsest parameterization of the quarter annulus can be given by NURBS of degrees $p_{\xi}=1$ and $p_{\eta}=2$, built upon the two knot vectors:

[TABLE]

and the 6 control points $P_{ij}$, $i=1,2$, and $j=1,2,3$, given in Table 1.

In what follows we will refer to this parametrization as $Q_{0}$.

In this section, together with $Q_{0}$, we will consider four 4-element parameterizations of the quarter annulus. Three will be called "uniform", denoted by letters $A$, $B$, and $D$, and one will be called "non-uniform", denoted by letter $C$. Note that the parameterization $D$ is obtained from the parameterization $A$ by setting all the weights to $1$. It does not exactly represent the quarter annulus, and it will only be used for solution approximation. The "uniform" parameterizations are obtained by performing the standard operation of knot-insertion on knot vectors $\Sigma_{0}$ and $\Pi_{0}$ to calculate the positions $P_{ij}$ and the weights $w_{ij}$ of the corresponding control points. We add knot value $s$ in $\xi$-direction and knot value $t$ in $\eta$-direction. Then the geometry is given by:

[TABLE]

The "non-uniform" parametrization is obtained by adding two knot values $s$ and $t$ to the knot vectors $\Sigma_{0}$ and $\Pi_{0}$, analogously to Eq. (24), but arbitrarily choosing the weights and the position of two inner control points: $P_{2,2}$ and $P_{2,3}$.

The control points $P_{ij}$ and the corresponding weights $w_{ij}$ for parameterizations $A_{1}$, $B_{1}$, $C_{1}$, and $D_{1}$ are listed in Tables 3, 4, 5, and 6, respectively. The higher order parameterizations $A_{2}$, $B_{2}$, $C_{2}$, and $D_{2}$ are obtained by standard degree elevation (by one in both directions) of parameterizations $A_{1}$, $B_{1}$, $C_{1}$, and $D_{1}$, respectively. Thereafter, the degrees of the parameterizations $A_{2}$, $B_{2}$, $C_{2}$, and $D_{2}$ will be given by $p_{\xi}=2$ and $p_{\eta}=3$. Moreover, for a special case, we will also consider a parameterization $D_{0}$, where $p_{\xi}=1$ and $p_{\eta}=1$, and which is obtained by degree reduction of parameterization $D_{1}$.

As discussed in the following section, a total of 19 patch tests for each of the problems of Section 3.1 (in two-dimensions) are studied. Therefore, for better readability and classification, we introduce the following notation. We will denote a patch test by $T^{\ell}_{G_{i},S_{j}}$, where the superscript $\ell=0,1,\ldots,5$, denotes the patch test case, $G_{i}$ denotes the basis for geometry parameterization, and $S_{j}$ denotes the basis for solution approximation. The basis $G_{i}$ and $S_{j}$ are chosen from the bases set given in Table 2, and their degree elevated/reduced versions.

3.3 Results and discussion of patch tests

We start the series of patch tests by choosing the coarsest parametrization of the quarter annulus, i.e. $Q_{0}$, paired with different bases for the solution.

The results of the first patch test are denoted by $T^{0}$ in Table 7. Note that all the combinations of basis functions fall within the isoparametric or sub-parametric finite element concept. However, the first two combinations pass the test, while the other two fail. Therefore, in order to understand the connection between the geometry and the solution bases, that leads to either passing or failing the patch test, we introduce the following examples, investigating other choices of the geometry parameterization and solution approximation.

In the next patch test, denoted by $T^{1}$ in Table 7, we demonstrate the performance of uniform parameterizations $A_{1}$ and $A_{2}$. All three combinations pass the patch test. Note that the test $T^{1}_{A_{1},A_{1}}$ is a standard IGA patch test within isoparametric concept, which can exactly represent a constant/linear solution on any basis which preserves partition of unity. To understand the relation between geometry representation and solution approximation, we substitute the geometry parametrization

[TABLE]

and the solution approximation

[TABLE]

into the solution

[TABLE]

and obtain the solution control variables as

[TABLE]

Now recall the following fundamental property of the degree elevation algorithm [47, Alg. 5.9, P.206] (which preserves the exactness of the geometry representation)

[TABLE]

The test $T^{1}_{A_{1},A_{2}}$ can be viewed analogous to the FEA sub-parametric approach (the degree of the solution approximation is higher than the degree of the geometry). If the solution is approximated by

[TABLE]

then, using (29), the corresponding control variables can be found as

[TABLE]

The test $T^{1}_{A_{2},A_{1}}$ can be viewed analogous to the FEA super-parametric approach (the degree of the solution approximation is lower than the degree of the geometry). However, the essential difference is that the geometry degree was artificially lifted from $A_{1}$ to $A_{2}$. Since using (29) in (27) yields the solution given by (28), this means that both the bases $A_{1}$ and $A_{2}$, as well as any two parameterizations with the same property, yield identical solutions for the solution control variables (within the numerical tolerance).

For the second patch test, denoted by $T^{2}$ in Table 7, we consider the situation when we have two uniform, but different, parameterizations of the domain, and one of them is used to parametrize the geometry and the second one is used as a basis to approximate the unknown solution. All the three combinations pass the test. Note that in this test, we used the same basis for the solution approximation as in the test $1$, but a different geometry parametrization $B_{1}$ (instead of $A_{1}$). However, both parameterizations, $A_{1}$ and $B_{1}$, were obtained from the coarsest parametrization $Q_{0}$ of the quarter annulus. According to the property of the knot insertion algorithm [47, Alg. 5.3, P.155], this gives

[TABLE]

Due to the property (32), geometry parameterizations $A_{1}$ and $B_{1}$ (as well as their higher degree versions $A_{2}$ and $B_{2}$, or any other two parameterizations with the same property) analytically yield the results for the solution control variables in $T^{2}$ to be identical (within the numerical tolerance) to the solution control variables in $T^{1}$, which is given by Eq. (28) (or for higher degree by Eq. (31). Note that the same property holds between bases $Q_{0}$ and $A_{1}$ in the test $T^{0}$, which implies that the solution in $T^{0}_{Q_{0},A_{1}}$ is expressed by Eq. (28), and by adding the property of degree elevation Eq. (29), we can conclude that the solution in $T^{0}_{Q_{0},A_{2}}$ is given by Eq. (31).

In the third patch test, denoted by $T^{3}$ in Table 7, we investigate the performance of non-uniform parametrization. The results of the third patch test show that all the three combinations pass the test. In fact, in this test only degree elevation algorithm (from $C_{1}$ to $C_{2}$) is used, which makes it fully analogous to $T^{1}$, albeit with a non-uniform parameterization.

In tests $T^{1}$, $T^{2}$ and $T^{3}$, the bases of solution and geometry were related, either both the bases were obtained from the geometry parameterization $Q^{0}$ (as in tests $T^{1}$ and $T^{2}$), or the solution basis was obtained from the geometry basis (as in test $T^{3}$). The idea of the fourth patch test, denoted by $T^{4}$ in Table 7, is to show the combinations of the geometry parameterizations and solution approximations which will fail to produce the linear solution. For this purpose we use the same approximation bases $A_{1}$ and $A_{2}$ for the solution, as in $T^{1}$ and $T^{2}$, but combined with the non-uniform geometry parameterization $C_{1}$ and $C_{2}$. As it is seen from Table 7, this test fails for all the combinations. Due to the fact that $C_{1}$ was not derived from the coarsest parameterization $Q_{0}$ of the quarter annulus, the property analogous to (32) between $A_{1}$ and $C_{1}$ does not hold.

In the fifth patch test, denoted by $T^{5}$ in Table 7, the conditions (29) and (32) between the geometry parameterization and the solution approximation do not hold, and therefore as expected, all combinations of $T^{5}$ fail. However, it is important for practical applications, and the suggested bases combinations will be used later in the numerical examples. We use the NURBS basis $A_{1}$ for geometry parameterization, and the B-splines basis $D_{i}$, $i,=\{0,1,2\}$ for solution approximation. Two important features of $T^{5}$ are the following:

• •

The sub-parametric case is different from what is considered in the previous patch tests, because instead of artificially lifting the degree of the geometry parametrization, we kept the minimum required order of the geometry basis, and reduced the degree of the solution approximation.

• •

In all previous patch tests, the basis used for the solution approximation, could also be used to exactly represent the geometry. However, in $T^{5}$, the B-splines basis, which is used for the solution, cannot exactly represent the geometry.

Conclusion: The sufficient condition for the patch test to pass consists in the requirement for the geometry and field bases to be equivalent, up to operations of degree elevation and/or knot insertion.

4 Numerical examples

In this section we present convergence results of various bases choices for the geometry and the numerical solution. We consider two classes of problems, namely, Laplace and linear elasticity.

4.1 Laplace equation in two-dimensions

For Laplace equation, we choose two geometries, the quarter annulus (see Fig. 4) which can be exactly represented by $\mathcal{N}_{1,2}$, and a wedge geometry (see Fig. 11) which is represented by $\mathcal{N}_{1,4}$.

4.1.1 Example 1

As a first numerical example, we consider the problem of Laplace equation in the quarter annulus domain of Fig. 4. We consider the exact solution

[TABLE]

which is prescribed as Dirichlet boundary condition on all boundaries. The convergence study cases are based on the choices of geometry and solution bases which are used in the patch tests in Section 3 (see Table 7). The convergence rates for all study cases are shown in Figures 8, 9, and 10. The results are combined as follows.

In Fig. 8, we have collected three choices of geometry parametrization, namely the coarsest $Q_{0}$, and its variations $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$. Since $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$ are equivalent to $Q_{0}$ up to knot-insertion and degree-elevation, the convergence results are identical, and the slope of the graph depends only on the lowest degree of the solution approximation ($1$ in the case of $A_{1}$, and $2$ in the case of $A_{2}$).

In Fig. 9, we show the results for the non-uniform geometry parameterizations $C_{1}$ and $C_{2}$, and the solution bases $C_{1}$, $C_{2}$, $A_{1}$, and $A_{2}$, in comparison with the coarsest geometry parameterization $Q_{0}$ paired with the solution bases $C_{1}$ and $C_{2}$. Despite the slight difference in the results associated with the choice of the geometry representation, all graphs exhibit the expected convergence rates, which is governed by the lowest degree of the solution approximation.

Finally, in Fig. 10, we included the results for the fifth case (represented by $T^{5}$ in Table 7), characterized by the fact that the solution B-splines basis can not represent the geometry exactly. The pairs $A_{1}$-$D_{0}$, $A_{1}$-$D_{1}$, $A_{1}$-$D_{2}$ are shown in comparison with $A_{1}$-$A_{1}$ and $A_{1}$-$A_{2}$. Despite the difference in weights of the basis functions between $A_{1}$-$A_{1}$ and $A_{1}$-$D_{1}$ (as well as between $A_{1}$-$D_{2}$ and $A_{1}$-$A_{2}$) the difference in the results is very minor. This makes us conclude that though the B-splines solution basis fails the patch test, it is nevertheless a suitable basis for the analysis. Moreover, the case $A_{1}$-$D_{0}$, which can be regarded as truly super-parametric in FEA context, also exhibits the expected convergence rate. This surprising result counters the established practice of FEM, where super-parametric approach is not recommended [71].

4.1.2 Example 2

In the second example, we study the Laplace problem in a wedge geometry, see Fig. 11. The choice of this geometry is motivated by the observation in the super-parametric case of Example 4.1.1. As opposed to the patch test cases where the degree for geometry representation was artificially lifted, the geometry in this case is given by a tensor product of NURBS of degrees $p_{\xi}=1$ and $p_{\eta}=4$. We consider the exact solution

[TABLE]

which is prescribed as Dirichlet boundary condition on all boundaries. In this example we focus on sub- and super- parametric solution approximations. As the geometry parameterization is fixed (we denote the corresponding NURBS basis by $\mathcal{N}_{1,4}$), to ease with the naming, we use the notation of paired bases as $(\mathcal{N}_{1,4},S_{k,l})$, where $S_{k,l}$ denotes the basis for solution approximation ($S=\mathcal{N}$ for NURBS and $S=\mathcal{B}$ for B-splines), $k$ denotes the degree in $\xi$ direction, and $l$ denotes the degree in $\eta$ direction. For sub-parametric solution approximations, we only elevate the degree in $\xi$ direction, i.e. $k=2$, and for super-parametric solution approximations, we only consider degree reduction in $\eta$ direction, i.e. $1\leq l\leq 3$. In Fig. 12, we present the results for the following $7$ cases: $(\mathcal{N}_{1,4},\mathcal{N}_{1,4})$, $(\mathcal{N}_{1,4},\mathcal{B}_{1,4})$, $(\mathcal{N}_{1,4},\mathcal{N}_{2,4})$, $(\mathcal{N}_{1,4},\mathcal{B}_{2,4})$, $(\mathcal{N}_{1,4},\mathcal{B}_{1,1})$, $(\mathcal{N}_{1,4},\mathcal{B}_{1,2})$, and $(\mathcal{N}_{1,4},\mathcal{B}_{1,3})$. As we can observe from the numerical studies presented in Fig. 12, together with the exact representation of the geometry, the convergence rate depends only on the lowest degree in the solution approximation, and the results for NURBS and B-splines solution bases are almost identical.

4.2 Linear elasticity in two-dimensions

For linear elasticity problem, we choose two geometries, thick walled cylinder (quarter annulus for symmetry boundary conditions) and a plate with a circular hole (see Fig. 17(a)).

4.2.1 Example 1 (pressurized cylinder)

In the first example of linear elasticity, we consider a problem of a thick walled cylinder under a uniform internal and external pressure, as shown in Fig. 13. Instead of the full problem, a quarter annulus is considered with the symmetry boundary conditions along $x=0$ and $y=0$, and the following boundary conditions on the inner and outer boundary of the cylinder are imposed:

[TABLE]

The solution to this problem in polar coordinates, under plane strain assumption, is given by

[TABLE]

where $s_{1}=\dfrac{r_{1}^{2}r_{2}^{2}(p_{2}-p_{1})}{r_{2}^{2}-r_{1}^{2}}$, $s_{2}=(1-2\nu)\dfrac{r_{1}^{2}p_{1}-r_{2}^{2}p_{2}}{r_{2}^{2}-r_{1}^{2}}$, and $s_{3}=\dfrac{1+\nu}{E}$. The numerical tests were performed on the same choices of the combination of geometry-solution bases, as in the first example for Laplace equation (see Table 7). The results of numerical simulations are organized in three plots: Fig. 14, 15 and 16. The obtained numerical solutions exhibit a similar pattern as that observed in Example 4.1.1, i.e. the graphs show that the convergence rate of the solution is defined by the lowest degree in the solution basis, independent of the (exact) geometry parameterization.

4.2.2 Example 2 (plate with a circular hole)

In the second example of linear elasticity, we consider a typical problem of a plate weakened by a circular hole of radius $a$, and subject to remote tension $T$ in $x$-direction. Only a quarter of a plate (of a finite size $L=4a$), as shown in Fig. 17(a), is modeled with the symmetry boundary conditions along $x=0$ and $y=0$, and analytical tractions prescribed on the rest of the boundary according to the solution below:

[TABLE]

where $\alpha_{0}=(1+\nu)Ta/(4E)$. In all study cases for this problem, the geometry is parametrized by a basis of second order with the following knot vectors:

[TABLE]

The corresponding control points are listed in Table 8. The parametrization consists of two elements (see Fig. 17(b)), and remains unchanged during the solution refinement process (denoted by $\mathcal{N}_{2,2}$). Following the notation of paired bases introduced in Section 4.1.2, we consider the following choices for the solution approximation basis: $(\mathcal{N}_{2,2},\mathcal{N}_{2,2})$, $(\mathcal{N}_{2,2},\mathcal{N}_{3,3})$, $(\mathcal{N}_{2,2},\mathcal{B}_{2,2})$, and $(\mathcal{N}_{2,2},\mathcal{B}_{3,3})$. Note, that B-Spline basis $\mathcal{B}_{2,2}$ is built on knot vectors (38). We also consider another geometry parameterization which is built on the knot vectors (38), with the boundary control points as listed in Table 8, but the weights of two inner points being changed to $w_{2,2}=w_{2,3}=0.9$. This parameterization will be denoted by $\tilde{\mathcal{N}}_{2,2}$. Thereafter, we consider two more choices for the solution approximation basis: $(\mathcal{N}_{2,2},\tilde{\mathcal{N}}_{2,2})$, and $(\mathcal{N}_{2,2},\tilde{\mathcal{N}}_{3,3})$. The results of the these study cases are shown in Fig. 18, where it can be seen that the convergence rate for all choices of the solution basis depends only on the order of the solution basis, and for the bases of the same order the results for different basis functions are almost identical.

4.3 Linear elasticity in three-dimensions

For three-dimensional linear elasticity problem, we choose a thick walled sphere (one-eighth of a hollow sphere for symmetry boundary conditions) geometry with inner radius $r_{1}$ and outer radius $r_{2}$, see Fig. 19. The coarsest NURBS parametrization of this geometry can be given by degrees $p_{\xi}=1$, $p_{\eta}=2$ and $p_{\zeta}=2$, built upon the knot vectors:

[TABLE]

with the control points listed in Table 9. We will refer to this parametrization as $Q_{1}$. Using the notations introduced earlier, the basis for this parameterization will be denoted by $\mathcal{N}_{1,2,2}$.

We also consider the following choices of the bases for solution approximation: $\mathcal{N}_{1,2,2}$, $\mathcal{N}_{2,2,2}$, $\mathcal{B}_{1,2,2}$, $\mathcal{B}_{2,2,2}$, and $\mathcal{B}_{1,1,1}$. Note that, we use uniform refinement of $Q_{1}$ for these bases (similar to the case $A_{1}$ of Section 3.2). And, the basis $\mathcal{B}_{1,1,1}$ is built on the knot vectors:

[TABLE]

The problem in consideration is a well known thick-walled pressurized sphere, i.e. the hollow sphere subjected to the following boundary conditions (in spherical coordinates $r,\theta,\phi$):

[TABLE]

On the rest of the boundary the symmetry conditions are prescribed. The analytical solution is given by

[TABLE]

where $\alpha_{1}=2(p_{1}r_{1}^{3}-p_{2}r_{2}^{3})(1-2\nu)$, $\alpha_{2}=(p_{1}-p_{2})(1+\nu)r_{1}^{3}r_{2}^{3}$, and $\alpha_{3}=2E(r_{2}^{3}-r_{1}^{3})$. The numerical results are shown in Fig. 20, where it can again be seen that the convergence rate in all the five cases depends only on the approximation basis for the numerical solution.

4.4 Numerical solution with PHT splines

To further demonstrate the effectiveness of the proposed method, and the use of weak coupling of bases for geometry and simulation, in this section we present numerical results using PHT-splines (to avoid the disruption of this section’s readability, the details of PHT-splines based construction is provided in Appendix A).

We again consider an annulus region in two-dimensions, see Fig. 4, described by a quadratic $C^{1}$ NURBS surface with $6\times 6$ control points, and following knot vectors on the parametric domain

[TABLE]

We consider the Poisson problem, and choose the source function such that the exact solution of the problem has the following form

[TABLE]

where

[TABLE]

In Fig. 21, we show the results during local refinement operations. From the first row to the fourth row, we show the T-mesh on the parametric domain (left), the numerical solution with $C^{1}$ PHT-splines representation (middle), and the corresponding exact error color-map (right) on the parametric domain. In Fig. 22, we present the convergence behavior of four choices of geometry and simulation bases, which are as follows:

• •

IGA with cubic NURBS (for the geometry as well as numerical solution). Note that, a quadratic NURBS is sufficient for this geometry, however, to have a fair comparison with the remaining studies, we elevate the degree while maintaining the exact geometry representation.

• •

IGA with cubic PHT-splines (for the geometry as well as the numerical solution). Note that, in this case, the computational geometry is only approximate (not exact as in IGA with cubic NURBS).

• •

GIFT with cubic B-splines for the numerical solution, and quadratic NURBS for exact geometry representation.

• •

GIFT with cubic PHT-splines for the numerical solution, and quadratic NURBS for exact geometry representation.

From the convergence plots, we conclude the following:

As observed in earlier studies, IGA with cubic NURBS and GIFT with cubic B-splines exhibit same convergence rate. This is because, with the exact geometry representation, both the solution bases are of the same degree. 2. 2.

Owing to the local adaptive refinement of the solution basis, both the cases of PHT-splines (IGA as well as the proposed GIFT) exhibit higher convergence rate than the tensor-product based solution bases (B-splines as well as NURBS). 3. 3.

The comparison of IGA with PHT-splines and GIFT with PHT-splines solution highlights an important difference. The advantage of the exact geometry representation in the latter case over an approximate geometry in the former case is very minor. This is due to the fact that the geometry of the computational domain is simple, and can be accurately approximated with the PHT-splines of third degree (geometry approximation error is below the discretization error). However, in realistic industrial problems with complex domains, this advantage will become more pronounced. In comparison to standard IGA with PHT splines, employing the exact coarse NURBS geometry parametrization in GIFT (together with PHT-splines solution) brings two distinct advantages:

• •

It eliminates the need to communicate with the original CAD model at each step of the solution refinement process, and the approximation of the boundaries.

• •

It also eliminates the need to refine the original coarse geometry, as well as to store and process the refined data, which can lead to significant computational savings for big problems.

5 Discussion on numerical results

The numerical results of Section 3, indicate that a sufficient condition for the combination of the geometry and the solution bases to pass the classical patch test is the requirement that these two bases are equivalent up to the degree elevation (29) and/or knot insertion operations (32), which ensure the preservation of the geometry exactness. However, we see from the numerical results in Sections 4.1-4.4, that despite failing the classical patch test, all cases presented in Table 7 exhibit optimal order convergence. This includes the bases derived from the geometry parametrization $A$ and $B$, the non-uniform parametrization $C$, and the B-splines basis $D$ (which cannot represent the geometry exactly). This is not surprising, however. Thanks to [1, Lemma 3.4-3.5, Theorem 3.1-3.2], with a suitable projector $\Pi_{\mathcal{V}_{h}}$, the optimal order global error estimates hold for every function $v\in H^{l}(\Omega)$. The effect of geometry mapping is included in the definition of the element size $h_{K}=\|\nabla F\|_{L^{\infty}(Q)}h_{Q}$, where $F$ is the geometry mapping, $Q$ is the element in the parameteric domain, and $h_{Q}$ is the element size in the parameteric domain. Moreover, the constant $C_{\text{shape}}$ appearing in [1, Lemma 3.4-3.5, Theorem 3.1-3.2] depends only on the weight function $w$ and its reciprocal $1/w$ on $\tilde{Q}$, where $\tilde{Q}$ is the support extension of $Q$, and is uniformly bounded with respect to the mesh size. This observation is in alignment with [60] that the classical patch test is neither necessary nor sufficient for convergence. Therefore, together with exact geometric representation, a suitable basis (derived-from or related-to geometry parameterization) can be used for optimal order convergence. For a quick comparison, in Table 10 we recall all the cases studied in this paper.

Nevertheless, it is important to exercise caution while devising a basis for the numerical solution. In all our test cases of Sections 4.1-4.3, the solution basis was constructed based on the same knot vectors as the geometry parametrization. This assures the continuity of the geometry parametrization within solution elements. For Example 4.2.2, we now consider the coarsest basis consisting of B-splines of degree $2\times 2$ defined on the knot vectors (denoted by $\tilde{\mathcal{B}}_{2,2}$):

[TABLE]

In this case, the singular point $(x=L,~{}y=L)$ is inside an element of any mesh for the numerical solution. As it can be seen from Fig. 23, the solution does not exhibit the expected convergence rate.

A thorough mathematical derivation of the presented approach is a subject of future research, but it is recommended to avoid those combination of geometry-solution bases where the continuity of the geometry parametrization is violated in the solution elements.

5.1 Naming convention

Let $p_{u}$ and $p_{g}$ denote the degrees of the basis functions for the solution, and the geometry, respectively. In the standard FE context, since the primary quantity of interest is the numerical solution, the naming convention is with respect to the parameter $p_{u}$. However, in IGA, the primary quantity is the geometry. Accordingly, for better readability, we propose the naming convention with respect to the parameter $p_{g}$, as presented in Table 11.

Note that, in standard FEM, the use of $p_{u}<p_{g}$ (super-parameteric) case is not recommended [71, P.172]. However, the results presented in Section 4.1.2 show that, with exact geometry representation, super-parametric approximations can also deliver optimal orders of convergence.

6 Conclusions

We presented a method which relaxes the requirement for a tight coupling between the spaces for the representation of the geometry and the approximation of the field variables, but retains both geometrical exactness and the ability to operate directly from CAD files.

The increased independence of the choice of the boundary parameterization and field approximations enables local refinement to capture sharp gradients, without modifying the geometry parameterization. We formulate recommendations regarding the relative choice of discretizations which, as our numerical results indicate, yield optimal convergence.

Important future work includes a detailed mathematical analysis of the method to prove a priori error estimates, and the investigation of this approach for other types of partial differential equations such as wave propagation. A detailed numerical analysis of the suitability of GIFT in a boundary element approach will follow the work of [6, 33, 4, 8, 7, 34, 69, 5].

Appendix A PHT-splines

In this section we present some details on PHT-splines, which were used for the numerical results in Section 4.4. PHT-splines, proposed by Deng et al. [19], are piecewise bicubic polynomials over a hierarchical T-mesh, which inherit the advantageous properties of T-splines. Unlike T-splines, PHT-splines are non-rational polynomial splines, and the refinement algorithm of PHT-splines is local and simple. The blending functions of PHT splines are linearly independent, an important property needed for finite element approximations.

A T-mesh is a rectangular partition of a planar domain with grid lines parallel to the boundary of the domain which allows T-junctions. In T-meshes, the end points of each grid line must lie on two other grid lines, and each cell or facet in the grid must be a rectangle. If a vertex is inside of the domain, it is called an interior vertex, otherwise, it is called a boundary vertex. There are two types of interior vertices, namely crossing vertices (i.e., it possesses valency 4) and T-vertices with valency 3, respectively.

Let $\Omega\in\mathbb{R}^{m}$ be a rectangular domain with boundary $\partial\Omega$. Denote by $\mathbb{T}=\cup\mathcal{K}$ a hierarchical T-mesh over domain $\Omega$, where $\mathcal{K}$ is a cell of the mesh. We further define a spline space

[TABLE]

where $\mathbb{P}_{pq}$ is the space of all the bivariate polynomials with degree $(p,q)$, and the space $\mathcal{C}^{\alpha,\beta}(\Omega)$ consists of all continuous bivariate spline functions up to order $\alpha$ in the $\xi$-direction and order $\beta$ in the $\eta$-direction. The dimension formula for the spline space $S(p,q,\alpha,\beta,\mathbb{T})$, with $p\geq 2\alpha+1$ and $q\geq 2\beta+1$, has already been provided in [19]. For the cubic PHT-splines space, the dimension formula can be written as

[TABLE]

where $V^{b}$ stands for boundary vertices and $V^{+}$ stands for interior crossing vertices. From the dimension formula, four basis functions are associated with each basis vertex (boundary vertex or crossing vertex), and they can be built with a hierarchical approach.

For the initial level, i.e., level [math], denoted as $\mathbb{T}_{0}$, the standard bi-cubic tensor-product B-splines are used as basis functions. For simplicity, we set the initial mesh to be a uniform rectangular grid. Suppose that the grid is $[\xi_{1},\xi_{2},\xi_{3},...,\xi_{s}]\times[\eta_{1},\eta_{2},\eta_{3},...,\eta_{t}]$. Since all the vertices are either crossing vertices or boundary vertices, there are four basis functions to be defined on any vertex $(\xi_{i},\eta_{j})$. Each basis function at $(\xi_{i},\eta_{j})$ has support $[\xi_{i-1},\xi_{i+1}]\times[\eta_{j-1},\eta_{j+1}]$. These four basis functions are defined to be the B-splines basis functions with knots

[TABLE]

respectively, such that their function values and derivatives vanish outside $[\xi_{i-1},\xi_{i+1}]\times[\eta_{j-1},\eta_{j+1}]$.

In the GIFT framework with PHT-splines, the computational domain is in the NURBS space whereas the numerical solution is in the PHT-splines form. Firstly, we construct the parametric domain of the PHT-splines model to represent the numerical solution. We can then get the initial numerical solution, in which the unknown control variables can be solved by the method presented in Section 2. By using an a posteriori error estimation technique, the supporting cell with large errors in the numerical solution (in the parametric domain) can be marked, and local $h$-refinement is performed only on the numerical solution. Several local refinement steps can be performed until the desired error level is achieved.

A.1 Construction of common parametric domain for NURBS geometry and PHT-splines solution

As a first step, we need to construct a parametric domain of the numerical solution. This can be seen as a preprocessing stage. Depending on the parametric domain of the NURBS patches under consideration, two kinds of operations are proposed in this step:

• •

If the computational domain is made of a single NURBS patch (Fig. 24(a)), or multi-patches in which the parametric domain of each patch forms a quad-mesh as presented in Fig. 24(b), then the parametric domain of the numerical solution is constructed as the partition-mesh formed by the knot lines on the parametric domain of the NURBS patches.

• •

For other cases, such as the examples shown in Fig. 24(c), we need to construct the initial parametric mesh $\mathbb{T}$ according to the topological connection information of the NURBS patches. A reparameterization operation should be performed in this case as described below.

Here we present an example to show the reparameterization method. Suppose that the patch $\textbf{\emph{S}}_{1}(u,v)$ in Fig. 24(c) has parametric domain $\mathcal{P}_{1}=[a_{1},b_{1}]\times[c_{1},d_{1}]$, the patch $\textbf{\emph{S}}_{2}(u,v)$ has parametric domain $\mathcal{P}_{2}=[a_{1},b_{1}]\times[c_{2},d_{2}]$, and the patch $\textbf{\emph{S}}_{3}(u,v)$ has parametric domain $\mathcal{P}_{3}=[a_{3},b_{3}]\times[c_{1},d_{1}]$. Then according to the topological information of all three patches, the parametric domain for the PHT-splines numerical solution should be constructed as shown in Fig. 24(c). In other words, the parametric domain of $\textbf{\emph{S}}_{1}(u,v)$ does not change, while the parametric domain of $\textbf{\emph{S}}_{2}(u,v)$ changes to $[a_{1},b_{1}]\times[d_{1},d_{1}+(d_{2}-c_{2})]$, and the parametric domain of $\textbf{\emph{S}}_{3}(u,v)$ changes to $[b_{1},b_{1}+(b_{3}-a_{3})]\times[c_{1},d_{1}]$.

In order to maintain the geometry of the patches during the transformation of parametric domain, a reparameterization technique should be adopted to obtain the new parametric representation of each patch. Suppose that the initial parametric domain of $\textbf{\emph{S}}(u,v)$ is $[a,b]\times[c,d]$, we can then use the following parameter transformation to achieve a new parameterization $\boldsymbol{F}(\xi,\eta)$ with parametric domain $[e,f]\times[g,h]$

[TABLE]

Remark 1

There is actually no need to derive the explicit parametric representation of the reparameterized surface $\boldsymbol{F}(\xi,\eta)$. The derivative information required in the solving stage, as shown in (13), can be evaluated from the original parameterization $\textbf{\emph{S}}(u,v)$ through the Jacobian transformation matrix.

After the planar T-mesh $\mathbb{T}$ is constructed as the parametric domain of PHT-splines, we can write the initial formula of the numerical solution in PHT-splines form as follows

[TABLE]

where $n$ is the number of basis vertices on $\mathbb{T}$, and $T_{i}$, $i=1,2,\cdots,4n$, are the control variables to be solved.

A.2 Solution process and local refinement

After constructing PHT-splines representation of the numerical solution, the control variables in (47) is obtained by solving the problem as introduced in Section 2.

In order to obtain a solution with desired accuracy, refinement operations are often needed to increase the approximation power of the solution space locally, and thereby optimize the computational expenses required for a given accuracy level. Contrary to the classical IGA, the refinement operation in GIFT is only performed on the numerical solution, while the parameterization of the computational domain remains the same during the refinement process. Choosing a PHT-splines approximation offers a natural local refinement. Local refinement is performed on the sub-patches of the numerical solution, on which the local error indicator is larger than a given marking threshold. The corresponding working flow is summarized in Algorithm 1.

A.2.1 Residual-based error indicator

Suppose that $u_{h}$ is the PHT-splines solution of the problem (9) and (12) by using the GIFT method, and $u$ is the exact solution. Let $e_{h}=u-u_{h}$ be the error of the GIFT approximation $u_{h}$. As the refinement operation in GIFT is only performed on the numerical solution, in order to determine the parametric cell to be refined, it is required to give an error indicator on each cell of the T-mesh in the parametric domain of the numerical solution, rather than an error indicator on the sub-patch with respect to the knot span on the NURBS parameterization $\boldsymbol{F}(\xi,\eta)$.

Suppose that ${\mathcal{K}}$ is the cell on the T-mesh $\mathbb{T}$ of the parametric domain $\mathcal{P}$ for the PHT-splines solution $u_{h}$. The residual-based a posteriori error estimate $\|e_{h}\|^{2}_{\mathcal{P}}$ over the parametric domain $\mathcal{P}$ proposed in [65, 68] can be rewritten as follows:

[TABLE]

where $\boldsymbol{x}={\boldsymbol{F}}(\xi,\eta)=(x(\xi,\eta),y(\xi,\eta))$, $C$ is a positive constant, and $h_{\mathcal{K}}$ is the circumference of the sub-patch in the NURBS parameterization $\boldsymbol{F}(\xi,\eta)$ of $\Omega$ with respect to the cell $\mathcal{K}$ in the parametric domain $\mathcal{P}$.

From (48), we can obtain the local error indicator $e_{\mathcal{K}}$ on each parametric cell $\mathcal{K}$ as follows

[TABLE]

If the parametric cell ${\mathcal{K}}$ is written as $[\xi_{0},\xi_{1}]\times[\eta_{0},\eta_{1}]$, we have then

[TABLE]

A.2.2 Marking strategies

The local refinement in GIFT requires a marking strategy to decide which elements should be refined. That is, we should determine a subset $\tilde{\mathbb{T}}$ of the parametric T-mesh $\mathbb{T}$ consisting of all those cells $\mathcal{K}$ that must be refined when the value of $e_{\mathcal{K}}$ is bigger than certain threshold. To determine the set $\tilde{\mathbb{T}}$, we use the mean-value strategy shown in Algorithm 2.

For some physical problems, highly varying distribution of estimated errors may be observed in the following cases (I) very few cells having extremely large estimated errors, (II) some cells having extremely small estimated errors, and (III) other cells having estimated errors that are much smaller than that of (I) while much larger than that of (II). We group the cells into three categories (I), (II) and (III), respectively. In this case, if only the parametric cells in the first group are refined by using the above marking algorithm, it would deteriorate the performance of the local $h$-refinement method. In this paper, a simple modification is proposed. For a given small percentage $\epsilon$, the $\epsilon\%$ parametric cells are firstly marked with largest estimated error for refinement, and then the mean-value marking approach is further applied to the remaining parametric cells.

A.2.3 Local refinement

In order to achieve accurate and efficient simulation results, local refinement should be performed on the marked parametric cells. In this part, the refinement rules of PHT-splines will be introduced.

Suppose that the T-mesh on the parametric domain at level $k$ is denoted by $\mathbb{T}_{k}$, and the PHT-splines basis functions on $\mathbb{T}_{k}$ are constructed as ${M_{j}^{k}}$, $j=1,...,d_{k}$. Then the basis functions on $\mathbb{T}_{k+1}$ can be constructed as follows: some basis functions are from the modifications of the former basis functions on $\mathbb{T}_{k}$, and others are from the new basis functions associated with the new basis vertices of $\mathbb{T}_{k+1}$. We represent a PHT-splines basis function by specifying 16 Bézier coefficients in each cell within the compact support of the basis function. When a cross vertex is added in a cell, the cell can be refined into four subcells. Each subcell supports the original basis function, and also has 16 Bézier coefficients. Five new vertices are obtained by adding a cross vertex, and some new basis vertices are also introduced. Hence, for the old basis functions, all the Bézier coefficients associated with the new basis vertices should be reset to zero. The new basis vertices are introduced from two sources, i.e., some are crossing vertices while others are T-vertices from the previous level. The latter become basis vertices as the addition of cross points to the neighboring cells. The new basis functions can be further constructed over their supporting cells as in the initial level. For further details, the readers can refer to [19].

Note that in the proposed GIFT framework, the local refinement is only performed on the numerical solution, hence we only need to update the PHT-splines basis functions on the locally-refined T-mesh for the new numerical solution. It is not necessary to derive the updated control variables from the old control variables after local refinement, and the new control variables for the refined numerical solution can be obtained by re-solving the PDE in the new PHT-splines space over the refined T-mesh.

Acknowledgements

G. Xu is supported by the National Nature Science Foundation of China under Grant Nos. 61472111, Zhejiang Provincial Natural Science Foundation of China under Grant Nos. LR16F020003, LQ16F020005, and the Open Project Program of the State Key Lab of CAD&CG (A1703), Zhejiang University.

S.P.A. Bordas and S. Tomar thank the financial support of the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578) entitled “Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery”, and the support of the Luxembourg National Research Funds INTER/FWO/15/10318764 and INTER/MOBILITY/14/8813215/CBM/Bordas. E. Atroshchenko and S.P.A. Bordas acknowledge the financial support of University of Luxembourg through its Computational Sciences research priority.