Manifold-based isogeometric analysis basis functions with prescribed sharp features

TL;DR

This paper develops manifold-based basis functions for isogeometric analysis that can handle surfaces with sharp features, boundaries, and creases, enabling accurate analysis of complex shell geometries.

Contribution

It introduces a novel manifold-based basis function construction combining differential geometry and conformal mappings, accommodating sharp creases and boundaries in isogeometric analysis.

Findings

Basis functions handle arbitrary sharp features and boundaries.

Suitable for Kirchhoff-Love thin shell analysis with kinks.

No normalization needed for partition of unity functions.

Abstract

We introduce manifold-based basis functions for isogeometric analysis of surfaces with arbitrary smoothness, prescribed $C^0$ continuous creases and boundaries. The utility of the manifold-based surface construction techniques in isogeometric analysis was demonstrated in Majeed and Cirak (CMAME, 2017). The respective basis functions are derived by combining differential-geometric manifold techniques with conformal parametrisations and the partition of unity method. The connectivity of a given unstructured quadrilateral control mesh in $\mathbb R^3$ is used to define a set of overlapping charts. Each vertex with its attached elements is assigned a corresponding conformally parametrised planar chart domain in $\mathbb R^2$, so that a quadrilateral element is present on four different charts. On the collection of unconnected chart domains, the partition of unity method is used for approximation. The transition functions required for navigating between the chart domains are composed out of conformal maps. The necessary smooth partition of unity, or blending, functions for the charts are assembled from tensor-product B-spline pieces and require in contrast to earlier constructions no normalisation. Creases are introduced across user tagged edges of the control mesh. Planar chart domains that include creased edges or are adjacent to the domain boundary require special local polynomial approximants. Three different types of chart domain geometries are necessary to consider boundaries and arbitrary number and arrangement of creases. The new chart domain geometries are chosen so that it becomes trivial to establish local polynomial approximants that are always $C^0$ continuous across tagged edges. The derived non-rational manifold-based basis functions are particularly well suited for isogeometric analysis of Kirchhoff-Love thin shells with kinks.

Equations

$$\begin{array}{rl}{\Psi }_1:& \eta \in \square \mapsto {\xi }_1\in [0,1]\in {\hat{\Omega }}_1,\\ {\Psi }_2:& \eta \in \square \mapsto {\xi }_2\in [-1,0]\in {\hat{\Omega }}_2.\end{array}$$

$${\Psi }_2\circ {\Psi }_1^{-1}:{\xi }_1\in {\hat{\Omega }}_1\mapsto {\xi }_2\in {\hat{\Omega }}_2.$$

$$f_j({\xi }_j)={\mathbf{p}}_j^{\mathsf{T}}({\xi }_j){\mathbf{\alpha }}_j\text{with}j\in \{1,2\},$$

$$\sum _{j=1}^2{w}_j({\xi }_j)=1\forall \eta \in \square \text{with}{\xi }_1={\Psi }_1(\eta ),{\xi }_2={\Psi }_2(\eta ).$$

$$\begin{array}{rl}{\tilde{w}}_1(\eta )& =B_1(\eta )+B_2(\eta )+B_3(\eta )+\frac{B_4(\eta )}{2},\\ {\tilde{w}}_2(\eta )& =\frac{B_4(\eta )}{2}+B_5(\eta )+B_6(\eta )+B_7(\eta ),\end{array}$$

$$w_j({\xi }_j)={\tilde{w}}_j({\Psi }_j^{-1}({\xi }_j)).$$

$$f(\eta )=\sum _{j=1}^2{w}_j({\xi }_j){\mathbf{p}}_j^{\mathsf{T}}({\xi }_j){\mathbf{\alpha }}_j\text{with}{\xi }_j={\Psi }_j(\eta ).$$

$${\mathbf{\alpha }}_j={\mathbf{P}}_j\mathbf{f},$$

$$f(\eta )=\sum _{j=1}^2\left(w_j({\xi }_j){\mathbf{p}}_j^{\mathsf{T}}({\xi }_j){\mathbf{P}}_j\right)\mathbf{f}={\mathbf{N}}^{\mathsf{T}}(\eta )\mathbf{f}\text{with}{\xi }_j={\Psi }_j(\eta ),$$

$$\frac{d^l{w}_j(-1)}{d{\xi }_j}=\frac{d^l{w}_j(1)}{d{\xi }_j}=0\forall l\le k.$$

$$\frac{d^l{\Psi }_1(0)}{d{\eta }^l}=\frac{d^l{\Psi }_2(1)}{d{\eta }^l}\forall l\le k.$$

$$\mathbf{x}(\eta )=\sum _I{N}_I(\eta ){\mathbf{x}}_I.$$

$${\Psi }_j^{\text{I}}:\mathbf{\eta }=({\eta }^1,{\eta }^2)\in \square \mapsto \mathbf{\xi }=({\xi }_j^1,{\xi }_j^2)\in {\hat{\Omega }}_j.$$

$$f(\mathbf{\eta })=\sum _{j=1}^4{w}_j({\mathbf{\xi }}_j)f_j({\mathbf{\xi }}_j)=\sum _{j=1}^4{w}_j({\mathbf{\xi }}_j){\mathbf{p}}_j^{\mathsf{T}}({\mathbf{\xi }}_j){\mathbf{\alpha }}_j\text{with}{\mathbf{\xi }}_j={\Psi }_j^{\text{I}}(\mathbf{\eta }),$$

$$f(\mathbf{\eta })=\sum _{j=1}^4\left(w_j({\mathbf{\xi }}_j){\mathbf{p}}_j^{\mathsf{T}}({\mathbf{\xi }}_j){\mathbf{A}}_j{\mathbf{P}}_j\right)\mathbf{f},$$

$${\mathbf{N}}^{\mathsf{T}}(\mathbf{\eta })=\sum _{j=1}^4{w}_j({\mathbf{\xi }}_j){\mathbf{p}}_j^{\mathsf{T}}({\mathbf{\xi }}_j){\mathbf{A}}_j{\mathbf{P}}_j.$$

$$z={\eta }^1+i{\eta }^2=|z|(\cos \varphi +i\sin \varphi )=|z|e^{i\varphi }\text{with}|z|=\sqrt{({\eta }^1{)}^2+({\eta }^2{)}^2}\text{and}\varphi =\arctan ({\eta }^2/{\eta }^1),$$

$${\left[\begin{array}{cccc}{z}_1& z_2& z_3& z_4\end{array}\right]}^{\mathsf{T}}={\left[\begin{array}{cccc}0+0i& 1+0i& 1+1i& 0+1i\end{array}\right]}^{\mathsf{T}}.$$

$${\Psi }_j^{\text{I}}={\mathbf{\Lambda }}_{\mathbf{qr}}^{\text{I}}\circ {\mathbf{\Lambda }}_{\mathbf{tr}}.$$

$${\mathbf{\Lambda }}_{\mathbf{tr}}(z;j)=(z-z_j)e^{-i\pi (j-1)/2},$$

$${\mathbf{\Lambda }}_{\mathbf{qr}}^{\text{I}}(z;v_j,n_j)=\frac{|z{|}^{\beta }}{|z{|}^{4/v_j}}{z}^{4/v_j}{e}^{i2\pi (n_j-1)/v_j}=|z{|}^{\beta }{e}^{i\varphi 4/v_j}{e}^{i2\pi (n_j-1)/v_j},$$

$$\begin{split}&f_{j}(\bm{\xi}_{j})=\begin{cases}f^{s_{j}=1}_{j}(\bm{\xi}_{j})&\text{ if }\xi^{2}_{j}\geq 0\\ f^{s_{j}=2}_{j}(\bm{\xi}_{j})&\text{ if }\xi^{2}_{j}<0\\ \end{cases}&\text{with}\quad f^{s_{j}=1}_{j}(\xi^{1}_{j},\,\xi^{2}_{j}=0)=f^{s_{j}=2}_{j}(\xi^{1}_{j},\,\xi^{2}_{j}=0)\,.\end{split}$$

$$f_j^{s_j}({\mathbf{\xi }}_j)={\mathbf{p}}_j^{s_j}({\mathbf{\xi }}_j)\cdot {\mathbf{\alpha }}_j^{s_j}.$$

$${\mathbf{p}}_j^{s_j}({\mathbf{\xi }}_j)=\mathbf{b}(\mathbf{\theta })\text{with}\mathbf{\theta }={\left({{\mathbf{\Lambda }}_{\mathbf{qr}}}^{\text{I}}\right)}^{-1}({\mathbf{\xi }}_j;k_j,s_j).$$

$${\Psi }_j^{\text{II}}={{\mathbf{\Lambda }}_{\mathbf{qr}}}^{\text{II}}\circ {\mathbf{\Lambda }}_{\mathbf{tr}},$$

$${{\mathbf{\Lambda }}_{\mathbf{qr}}}^{\text{II}}(\mathbf{\eta };k_j,s_j,l_s,m_s)={\mathbf{\Lambda }}_{\mathbf{qr}}^{\text{I}}(\mathbf{\eta };k_j{l}_s,(s_j-1)l_s+m_s).$$

$${\Psi }_j^{\text{III}}={{\mathbf{\Lambda }}_{\mathbf{qr}}}^{\text{III}}\circ {\mathbf{\Lambda }}_{\mathbf{tr}}$$

$$\begin{array}{rl}& {{\mathbf{\Lambda }}_{\mathbf{qr}}}^{\text{III}}(\mathbf{\eta };k_j,s_j,l_s,m_s)=\{\begin{array}{ll}{\mathbf{\Lambda }}_{\mathbf{qr}}^{\text{I}}(\mathbf{\eta };4(k_j-1)l_s,(s_j-1)l_s+m_s)& \text{ if }{s}_j<k_j\\ {\mathbf{\Lambda }}_{\mathbf{qr}}^{\text{I}}(\mathbf{\eta };4l_s/3,m_s)\cdot e^{i\pi /2}& \text{ if }{s}_j=k_j\end{array},\end{array}$$

$${\mathbf{x}}_I^S=\sum _I{L}_{IJ}{\mathbf{x}}_J,$$

$$\mathbf{x}(\mathbf{\eta })=\sum _I{N}_I(\mathbf{\eta }){\mathbf{x}}_I^S,$$

Full text

1 Introduction

Smooth approximation schemes for unstructured meshes are crucial for isogeometric design and analysis of parts with arbitrary topology. Until recently isogeometric analysis was dominated by NURBS basis functions, which is the prevailing technology in present CAD systems. To represent the bounding surface of parts with arbitrary topology CAD systems usually resort to trimmed NURBS. Trimming involves the computation of the intersection between spline surfaces with other surfaces or curves. The respective nonlinear root-finding problems look deceptively simple but are extremely hard to robustly solve and lead to non-watertight geometries [1, 2, 3]. In the analysis context, the non-watertight geometries obtained from trimming pose unique challenges. Without turning to trimming, unstructured meshes with extraordinary vertices, i.e. vertices with different than four attached patches inside the domain, are required to represent parts with arbitrary topology. In computer-aided design numerous techniques have been developed to deal with extraordinary vertices, including geometrically $G^{k}$ and parametrically $C^{k}$ continuous constructions and subdivision surfaces. The application and further development of these techniques is currently a very active area of research in isogeometric analysis, see e.g. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Unfortunately, none of the techniques from computer-aided design seems to give optimal finite element convergence rates without further modifications, especially when applied to Kirchhoff-Love thin shells with arbitrary geometry. Hence, the search for easy to implement and optimally convergent schemes is still open. In a complementary line of research, there has been progress in isogeometric analysis of shells based on trimmed surfaces and weak enforcement of mechanical continuity conditions across patch boundaries, see e.g. [16, 17].

Manifold techniques for mesh-based construction of $C^{k}$ continuous surfaces were first introduced in Grimm and Hughes [18]. In contrast to most other smooth surface construction techniques, which essentially rely on glueing of surface patches along their edges, in manifold techniques a surface is created by blending of overlapping surface patches. The surface patches are first defined over unconnected planar chart domains in $\mathbb{R}^{2}$ and subsequently mapped to $\mathbb{R}^{3}$. On the collection of the unconnected planar chart domains, the partition of unity method [19] is used to construct the smooth surface patches. To that end, on each planar chart domain a local polynomial approximant and a partition of unity, or blending, function is needed. In addition, transition functions are required in order to be able to navigate between the different chart domains. In manifold techniques the planar chart domains, the transition functions, the local polynomial approximants and the blending functions can all be relatively freely chosen, which makes them extremely versatile. Based on the seminal work of Grimm and Hughes a number of complementary mesh-based manifold approaches have been proposed in geometric modelling [18, 20, 21, 22, 23]. Most of these approaches differ in the choice of the planar chart domains. For instance, in [18] each vertex, edge and face of the mesh have a corresponding planar chart domain with a suitably chosen geometry. In [21] and [22] only the vertices have a corresponding planar chart domain in the form of a conformally mapped star-shaped polygonal disk or a circular disk, respectively. In [20] the chart domains are chosen similar to the characteristic map in subdivision surfaces [24]. Each of the mentioned choices for the chart domains implies a corresponding choice for the transition functions. Following the construction proposed in Ying and Zorin [21], Majeed and Cirak [25] introduced a new set of manifold-based basis functions which can yield high convergence rates in finite element analysis. From a finite element viewpoint, the manifold-based basis functions resemble spline basis functions in the sense that each basis function has a local support and has one corresponding vertex.

In this paper, we derive new manifold-based basis functions for the isogeometric design and analysis of smooth surfaces with boundaries and with $C^{0}$ continuous sharp features, like creases and corners, by extending [25]. The crease edges are tagged as such on the control mesh by the user and during the finite element analysis different mechanical continuity conditions can be imposed across the crease edges and along the boundary edges, see Figure 1. In mesh-based construction of surfaces using manifold techniques creases have previously been considered in Della Vecchia and Jüttler [26] and boundaries in Tosun and Zorin [27]. The essential idea can be reduced to the choice of special local polynomials in the partition of unity approximation. That is, the local polynomials on each chart domain have to consist out of several polynomial pieces that are $C^{0}$ continuously connected across the crease edges. As in [25], in our present construction each vertex and its attached elements have a corresponding star-shaped planar chart domain consisting of images of conformally mapped unit squares placed around a centre vertex. There is one chart domain per vertex and each quadrilateral element is present on four different chart domains. The choice of the special $C^{0}$ continuous local polynomials can be simplified by slightly modifying the geometry of the chart domains. The new chart domains are chosen such that they are rotationally symmetric with respect to the arrangement of crease edges. Or expressed differently, the crease edges must partition the chart domain into equiangular sectors. The $C^{0}$ continuously connected polynomial pieces can subsequently be obtained by mapping a tensor-product basis, like the Lagrange or Bernstein basis, into each of the sectors of the chart domain. The new chart domains are parameterised with a quasi-conformal map so that each element on it can have a different shape depending on the arrangement of the crease edges. In contrast to the conformal map used in [25] the proposed quasi-conformal map is not angle-preserving but still provides easily computable smooth transition functions. The proposed new construction leads to three different types of chart domains. One of them is specifically designed to deal with crease arrangements which lead to creased sectors with concave corners. As in [25] we assemble the smooth partition of unity blending functions from tensor-product B-spline segments defined over a unit square. However, similar to [27], the knot-interval and the B-spline coefficients are chosen such that the assembled blending functions do not require normalisation. Hence, the new manifold-based basis functions derived in this paper are non-rational.

The outline of this paper is as follows. In Section 2 the manifold-based basis functions introduced in [25] are briefly reviewed. The treatment of sharp features is discussed in Section 3. Depending on the arrangement of crease edges we distinguish between rotationally symmetric and asymmetric chart domains. In addition, we consider the case of crease arrangements leading to concave sectors on a chart domain. Each of the three cases requires a different quasi-conformal map for parametrisation. Section 4 introduces the finite element analysis of thin shells with normal control along boundaries and across crease edges. In Section 5 the new manifold-based basis functions are applied to several Bernoulli beam and linear and nonlinear Kirchhoff-Love shell problems. The numerical convergence of $L_{2}$ and energy norm errors with decreasing element size is demonstrated.

2 Review of manifold-based basis functions

In the following we briefly review the construction of univariate and bivariate manifold-based basis functions. The discussion is focused on their application in finite element analysis, so that the underlying manifold concepts from differential geometry and the partition of unity interpolation are only mentioned in passing. For a more comprehensive discussion on manifold-based basis functions and surface construction in geometry we refer to [25, 21, 18].

2.1 Univariate basis functions

It is instructive to first review the derivation of the univariate manifold-based basis functions. We consider the dash-dotted control polygon in Figure 2 representing a part of a finite element mesh consisting of vertices $\bm{x}_{I}\in\mathbb{R}^{3}$ and elements between consecutive vertices. Our aim is to derive the basis functions for a representative element $[\bm{x}_{I},\,\bm{x}_{I+1}]\in\mathbb{R}^{3}$, as highlighted in Figure 2. As in conventional finite elements we define a reference element $\Box\coloneqq[0,\,1]\in\mathbb{R}$ that will serve as an integration domain for evaluating the finite element integrals. In the manifold-based approach the basis functions are obtained by smoothly blending local polynomials defined over several overlapping charts, or patches as they were called in [25]111We refrain in this paper from using the term patches because in computer-aided design literature patches denote what are the elements in the finite element literature.. In the following two chart domains $\hat{\Omega}_{1}\coloneqq[-1,\,1]\in\mathbb{R}$ and $\hat{\Omega}_{2}\coloneqq[-1,\,1]\in\mathbb{R}$ are introduced for the element $[\bm{x}_{I},\,\bm{x}_{I+1}]$. $\hat{\Omega}_{1}$ is associated with the vertex $\bm{x}_{I}$ and its 1-neighbourhood and $\hat{\Omega}_{2}$ is associated with the vertex $\bm{x}_{I+1}$ and its 1-neighbourhood. The 1-neighbourhood of a vertex is defined as the union of elements that contain the vertex. Hence, the basis functions in the element $[\bm{x}_{I},\,\bm{x}_{I+1}]$ will be obtained by blending local polynomials defined over $\hat{\Omega}_{1}$ and $\hat{\Omega}_{2}$.

As shown in Figure 2, the reference element $\Box$ maps to elements in both chart domains $\hat{\Omega}_{1}$ and $\hat{\Omega}_{2}$. The respective maps are defined as

[TABLE]

In this specific univariate construction both maps are simple translations. The two maps imply a transition map between the two chart domains given by

[TABLE]

We choose on each chart domain a polynomial approximant of the form

[TABLE]

where the vectors $\bm{p}_{j}(\xi_{j})$ contain a local polynomial basis, like the power, Lagrange or Bernstein basis, and the vectors $\bm{\alpha}_{j}$ contain their coefficients. In the present example a Lagrange basis is chosen, see Figure 2. The local basis on the chart $\hat{\Omega}_{1}$ is quadratic and on the chart $\hat{\Omega}_{2}$ it is piecewise linear. In addition to the polynomial approximants $f_{j}(\xi_{j})$, on each chart a smooth blending, or weight function $w_{j}(\xi_{j})$ is chosen. The two non-zero blending functions over the reference element $\Box$ must satisfy the partition of unity property, i.e.,

[TABLE]

The blending functions can be conveniently combined from suitably defined B-spline basis functions. In Figure 3 the construction of $C^{2}$-continuous blending functions from cubic B-splines is illustrated. The uniform knot interval length of the B-splines is $1/4$. The blending functions are defined as the linear combinations

[TABLE]

and yield after mapping the blending functions

[TABLE]

Evidently, both functions satisfy the partition-of-unity property (4) given that a complete B-spline basis sums up to one. It is worth emphasising that the blending functions proposed here are polynomial in contrast to the rational ones in [21, 25], see also B. The tensor products of the new univariate blending functions yield the corresponding multivariate blending functions. Note that different from the chosen knot interval of 1/4, as proposed in [27], a knot interval of 1/3 yields also blending functions which sum up to one without normalisation.

Finally, the approximant over the reference element $\Box$ is obtained by blending the approximants over the two charts

[TABLE]

Due to the choice of the Lagrange basis for $\bm{p}_{j}(\xi_{j})$ the coefficients $\bm{\alpha}_{j}$ can be interpreted as vertex coefficients. The coefficients $\bm{\alpha}_{j}$ of each chart are formally obtained with

[TABLE]

where $\bm{f}$ is the array of coefficients of all control polygon vertices and $\bm{P}_{j}$ is a gather matrix, filled with ones and zeros. This introduced in (7) yields the manifold-based basis functions

[TABLE]

where $\bm{N}(\eta)$ is the array of non-zero basis functions over the considered element. The number of non-zero basis functions in $\bm{N}(\eta)$ depends on the cardinality of the set of vertices in the two chart domains $\hat{\Omega}_{1}$ and $\hat{\Omega}_{2}$, and the specific local approximants on them.

To investigate the smoothness of the basis functions $\bm{N}(\eta)$ it is necessary to consider their derivatives. According to definition (9), their smoothness depends on the blending functions $w_{j}(\xi_{j})$, the local basis functions $\bm{p}_{j}(\xi_{j})$ and the maps $\varPsi_{j}(\eta)$. To obtain $C^{k}$-continuous basis functions each of these functions has to be $k$-times differentiable over each chart domain; and, in addition, the blending functions $w_{j}(\xi_{j})$ and their up to k-th derivatives must vanish at the chart domain boundaries, i.e.,

[TABLE]

Moreover, the differentiability of the mappings $\varPsi_{1}(\eta)$ and $\varPsi_{2}(\eta)$ requires that they satisfy at the centre of the chart domain

[TABLE]

The smoothness of the basis functions $\bm{N}(\eta)$ can be pointwise reduced by selecting a suitable local polynomial basis $\bm{p}_{j}(\xi_{j})$, as in Figure 2 with a piecewise linear basis on chart $\hat{\Omega}_{2}$. The resulting basis functions $\bm{N}(\eta)$ and their first derivatives are plotted in Figure 4. The smoothness of the basis is reduced to $C^{0}$ at the vertex $\bm{x}_{I+1}$. This has also an effect on the number of non-zero basis functions in an element. For instance, in the two elements $[\bm{x}_{I}\,,\bm{x}_{I+1}]$ and $[\bm{x}_{I+1},\,\bm{x}_{I+2}]$ with the linear piecewise $C^{0}$ continuous local basis functions $\bm{p}_{j}(\xi_{j})$ there are three non-zero basis functions. In elements with a quadratic local basis on each chart there are four non-zero basis functions. Finally, with the generated manifold-based basis functions points on the reference element $\eta\in\Box$ are mapped onto the manifold according to

[TABLE]

2.2 Bivariate basis functions

The bivariate manifold-based basis functions provide smooth approximants even on unstructured surface meshes with extraordinary vertices. We consider the construction of the manifold-based basis functions for a representative element in the quadrilateral finite element mesh shown in Figure 6. The reference element is now defined as the unit square $\Box:=[0,\,1]\times[0,\,1]\in\mathbb{R}^{2}$. Similar to the univariate case each of the four vertices of the considered element and their 1-neighbourhoods have an associated chart domain $\hat{\Omega}_{j}\in\mathbb{R}^{2}$ with $j\in\{1,\,2,\,3,\,4\}$. The number of elements in a chart $\hat{\Omega}_{j}$ depends on the number of elements $v_{j}$ connected to the respective vertex, which is called the valence of the vertex. In quadrilateral meshes, the interior vertices with $v_{j}\neq 4$ are the extraordinary vertices. The basis functions will be obtained by smoothly blending local polynomials defined over each of the four overlapping charts $\hat{\Omega}_{j}$.

As shown in Figure 6, the reference element $\Box$ maps to elements in the four chart domains $\hat{\Omega}_{j}$ according to

[TABLE]

The construction of these maps requires special care and will be detailed further below. The bivariate approximant over the reference element $\Box$ is obtained by blending the approximants over the four charts

[TABLE]

where $w_{j}(\bm{\xi}_{j})$ and $\bm{p}_{j}(\bm{\xi}_{j})$ are the blending functions and the vector of local basis functions on the chart domain $\hat{\Omega}_{j}$, respectively. As mentioned in Section 2.1, the blending functions $w_{j}(\bm{\xi}_{j})$ are obtained by taking the tensor-product of univariate blending functions (5). Next, the coefficients $\bm{\alpha}_{j}$ are to be expressed in dependence of the vertex coefficients collected in the array $\bm{f}$. The number of unique vertices in the union of the four charts $\cup_{j=1}^{4}\hat{\Omega}_{j}$ is not fixed and depends on the valences $v_{j}$ of the four vertices of the element. Hence, a least-squares projection is applied to express (14) in dependence of the vertex coefficients

[TABLE]

where we used $\bm{\alpha}_{j}=\bm{A}_{j}\bm{P}_{j}\bm{f}$ on each chart $\hat{\Omega}_{j}$. Here, $\bm{P}_{j}$ is the gather matrix and $\bm{A}_{j}$ the least-squares projection matrix. It is evident that the number of polynomial coefficients in $\bm{\alpha}_{j}$ must not be greater than the number of vertices in the chart domain $\hat{\Omega}_{j}$. It is worth mentioning that the projection matrix $\bm{A}_{j}$ depends only on the valence $v_{j}$ of the vertex and the chosen local basis function $\bm{p}_{j}$ so that it can be precomputed and tabulated. Finally, the array of basis functions is defined with

[TABLE]

As in the univariate case, the smoothness of the basis functions $\bm{N}(\bm{\eta})$ depends on the blending functions $w_{j}(\bm{\xi}_{j})$, the local basis functions $\bm{p}_{j}(\bm{\xi}_{j})$ and the mappings $\bm{\varPsi}^{\textrm{I}}_{j}(\bm{\eta})$. In the presented construction conformal maps are critical for the smooth parametrisation of the chart domains $\hat{\Omega}_{j}$ and the composition of the blending functions $w_{j}(\bm{\xi}_{j})$. Specifically, the map $\bm{\varPsi}^{\textrm{I}}_{j}$ is composed of several conformal maps. To write conformal maps more succinctly, the coordinates $\bm{\eta}=(\eta^{1},\,\eta^{2})$ of points in the reference element $\Box$ are expressed as a complex number

[TABLE]

where $|z|$ is the radius and $\phi$ the angle in the complex plane. In the complex plane the coordinates of the four corners of the reference element are given by

[TABLE]

As illustrated in Figure 7, the map $\bm{\varPsi}^{\textrm{I}}_{j}\colon\bm{\eta}\in\Box\mapsto\bm{\xi}_{j}\in\hat{\Omega}_{j}$ is composed of two maps, i.e.,

[TABLE]

The auxiliary linear map $\bm{\Lambda_{tr}}$ is responsible for translating and rotating the reference element and is given by

[TABLE]

where $(z-z_{j})$ is a translation and $e^{-i\pi(j-1)/2}$ is a (rigid body) rotation. The quasi-conformal map $\bm{\Lambda^{\textrm{I}}_{qr}}$ maps the reference element into a wedge-shaped domain and is chosen as

[TABLE]

where $\beta$ is a free parameter, $v_{j}$ denotes the valence of the centre vertex and $n_{j}$ is the number of the sector onto which the reference element is mapped. According to the first expression in (21), the map $\bm{\Lambda^{\textrm{I}}_{qr}}$ is composed of a standard conformal map $z^{4/v_{j}}$, a scaling of the radius with $|z|^{\beta-4/v_{j}}$ and a rotation $e^{i2\pi(n_{j}-1)/v_{j}}$. The parameter $\beta$ can be chosen relatively freely [21]. For instance, when $\beta=4/v_{j}$ the scaling term drops and the map $\bm{\Lambda^{\textrm{I}}_{qr}}$ is composed of a conformal map and a rotation as in [25]. The second expression in (21) shows that the radius $|z|$ will remain constant when $\beta=1$. Hence, in this paper we choose $\beta=1$ to obtain a more uniform mapping $\bm{\varPsi}^{\textrm{I}}_{j}$ with smaller entries in its Jacobian matrix. As it will become clear, the choice $\beta=1$ also simplifies the parametrisation of the asymmetric chart domains and chart domains with concave corners in Sections 3.2 and 3.3. See A for a review on conformal maps and a discussion on the choice of $\beta$.

With the mapping $\bm{\varPsi}^{\textrm{I}}_{j}$ at hand, we can now evaluate the manifold-based basis functions defined in (16). For instance, for a given integration point $\bm{\eta}$ in the reference element, first its image $\bm{\xi}_{j}=\bm{\varPsi}^{\textrm{I}}_{j}(\bm{\eta})$ in each of the four overlapping charts is found and after that the blending functions $w_{j}(\bm{\xi}_{j})$ and local polynomials $\bm{p}_{j}(\bm{\xi}_{j})$ are evaluated.

3 Creased bivariate manifold-based basis functions

We introduce domain boundaries and $C^{0}$ continuous creases by modifying the local polynomials $\bm{p}(\bm{\xi}_{j})$, and, as necessary, the chart domains $\hat{\Omega}_{j}$. Creases are introduced along the edges of the control mesh prescribed by the user. The local polynomials $\bm{p}(\bm{\xi}_{j})$ are chosen to be $C^{0}$ continuous across the creases. The new chart domains require new maps, i.e. $\bm{\varPsi}^{\textrm{II}}_{j}\colon\bm{\eta}\in\Box\mapsto\bm{\xi}_{j}\in\hat{\Omega}_{j}$, from the reference element to the chart domains. In turn, the blending functions $w_{j}(\bm{\xi}_{j})$ are combined from tensor-product B-splines using the new maps $\bm{\varPsi}^{\textrm{II}}_{j}$. The remaining steps in manifold-based basis construction are identical to the non-creased case.

Depending on the number and arrangement of creases meeting at a vertex three different crease types are possible. As discussed in the following, each of them requires a different treatment. In Figure 8 the three different crease types are illustrated with the help of a sample geometry.

3.1 Crease Type 1: Rotationally symmetric chart domains

Type 1 creases have, as shown in Figure 9, a rotationally symmetric chart domain with respect to the arrangement of the creased edges. The centre vertex of a chart $\hat{\Omega}_{j}$ can have $k_{j}\geq 2$ attached edges tagged as crease222In this section, $k$ is used for the number of creased spoke edges, which is different from its meaning in $G^{k}$ and $C^{k}$ continuity.. The $k_{j}$ crease edges split the chart domain into $k_{j}$ equiangular sectors, and the number of elements in each of the sectors must be the same.

To introduce the basic idea in defining a suitable local approximant $f_{j}(\bm{\xi}_{j})$, at first the case with two creased edges, $k_{j}=2$ , as depicted in Figure 9a is considered. The crease splits the chart domain into the two sectors $s_{j}=1$ and $s_{j}=2$ along the $\xi_{j}^{1}$ coordinate axis. The local approximant is chosen as the piecewise continuous function

[TABLE]

As in the non-creased case both approximants are given by

[TABLE]

The coefficients $\bm{\alpha}_{j}^{s_{j}=1}$ and $\bm{\alpha}_{j}^{s_{j}=2}$ must be matched so that $f_{j}(\bm{\xi}_{j})$ is piecewise continuous. Especially for a tensor-product Lagrange basis $\bm{p}_{j}^{s_{j}}(\bm{\xi}_{j})$, it is straightforward to determine a new set of coefficients $\bm{\alpha}_{j}$ from the coefficients $\bm{\alpha}_{j}^{s_{j}=1}$ and $\bm{\alpha}_{j}^{s_{j}=2}$ which ensure piecewise continuity. The two Lagrange approximants must share the same coefficients along the common crease edge.

In chart domains with three or more creases, $k_{j}\geq 3$, the bases $\bm{p}^{s_{j}}_{j}(\bm{\xi}_{j})$ are obtained by mapping a tensor-product Lagrange basis $\bm{b}(\bm{\theta})$ to each sector. The respective local approximants $f^{s_{j}}_{j}(\bm{\xi}_{j})$ in two adjacent sectors meet $C^{0}$ continuously because they share the same coefficients along the common crease edge. The Lagrange basis $\bm{b}(\bm{\theta})$ is mapped onto a sector $s_{j}$ using the quasi-conformal map (21) introduced earlier

[TABLE]

This mapping is best understood in conjunction with Figure 9b. Each of the three sectors in a different colour represents the image of a square mapped with $\bm{\Lambda^{\textrm{I}}_{qr}}(\bm{\theta};\,k_{j},\,s_{j})\colon\bm{\theta}\in\Box\mapsto\bm{\xi}_{j}\in\hat{\Omega}_{j}$ with $s\in\{1,\,2,\,3\}$. Obviously, the polynomial degree of $\bm{b}(\bm{\theta})$ must be chosen such that it is less than or equal to the number of vertices in each sector. More vertices have to be added, e.g. by refining the elements in each chart domain, if a higher degree basis is requested. Although we chose in our implementation a Lagrange basis for $\bm{b}(\bm{\theta})$ it is possible to use other boundary interpolating tensor-product basis, like the Bernstein basis.

3.2 Crease Type 2: Asymmetric chart domains

Type 2 creases have a rotationally asymmetric chart domain with respect to the arrangement of crease edges. Again, the centre vertex of a chart $\hat{\Omega}_{j}$ can have $k_{j}\geq 2$ attached edges tagged as crease. However, the $k_{j}$ crease edges split the chart domain into $k_{j}$ sectors, which are not equiangular and have different number of elements. This makes it difficult to establish $C^{0}$ continuous local approximants $f_{j}(\bm{\xi}_{j})$ following the approach introduced for Type 1 creases.

Therefore, Type 2 creases are first mapped onto a chart domain, which is equiangular with respect to the arrangement of crease edges, see Figures 10 and 11. This is possible in a manifold-based approach because the chart domains $\hat{\Omega}_{j}$ can be freely chosen, subject to some smoothness and invertibility constraints. The respective maps $\bm{\varPsi}^{\textrm{II}}_{j}\colon\bm{\eta}\in\Box\mapsto\bm{\xi}_{j}\in\hat{\Omega}_{j}$ have the same structure like for non-creased charts (19), i.e.,

[TABLE]

where $\bm{\Lambda_{tr}}$ consists as defined in (20) of a translation and rotation and $\bm{\Lambda_{qr}}^{\textrm{II}}$ is a quasi-conformal map with the arguments of $\bm{\Lambda^{\textrm{I}}_{qr}}$ in (21) replaced according to

[TABLE]

Here, $k_{j}$ is the number of creases, $s_{j}\in\{1,\,2,\,\dotsc,\,k_{j}\}$ is the sector number, $l_{s}$ is the number of elements in the sector $s_{j}$ and $m_{s}\in\{1,\,\dotsc,\,l_{s}\}$ is the local element number in the sector $s_{j}$. The meaning of the introduced variables is further clarified in Figure 12. After the equiangular chart domain for Type 2 creases is established, as in Figures 10 and 11, the local approximants $f_{j}^{s_{j}}(\bm{\xi}_{j})$ in each sector are obtained following the approach for a Type 1 crease chart. That is, they are defined as in (22) in case of two creases or otherwise mapped according to (24) from a tensor-product Lagrange basis $\bm{b}(\bm{\theta})$. The local approximant $f_{j}(\bm{\xi}_{j})$, composed from the local approximants in each of the sectors $f_{j}^{s_{j}}(\bm{\xi}_{j})$, is $C^{0}$ continuous across the crease edges when the approximants share the same coefficients along the edges.

3.3 Crease Type 3: Chart domains with concave corners

The arrangement of creases in Type 2 chart domains can lead sometimes to non-convex sectors on the chart domain. However, there is no regular mapping from a convex to a non-convex domain with a non-zero determinant of the Jacobian. Therefore, when the map $\bm{\varPsi}^{\textrm{II}}_{j}$ defined in (25) is used in a non-convex sector it leads to an unwanted folding-over of the surface as visible in Figure 13. This behaviour is similar to the one observed with isoparametrically mapped non-convex Lagrange finite elements.

The folding-over can be avoided by composing the non-convex sector from several smoothly connected convex pieces. To this end, the chart domain depicted in Figure 14 is introduced. The respective map $\bm{\varPsi}^{\textrm{III}}_{j}\colon\bm{\eta}\in\Box\mapsto\bm{\xi}_{j}\in\hat{\Omega}_{j}$ has the same structure like the previously introduced maps, c.f. (19) and (25),

[TABLE]

with the quasi-conformal map $\bm{\Lambda_{qr}}^{\textrm{III}}$ obtained by replacing the arguments of $\bm{\Lambda^{\textrm{I}}_{qr}}$ in (21) in the following way

[TABLE]

where $k_{j}$ is the number of creases, $s_{j}\in\{1,\,2,\,\dotsc,\,k_{j}\}$ is the sector number, $l_{s}$ is the number of elements in a sector and $m_{s}$ is the local element number in a sector. All the mentioned variables have the same meaning as in Section 3.2 and are clarified in Figure 15. As shown in Figure 14, with the introduced map $\bm{\varPsi}^{\textrm{III}}_{j}$ the concave sector with $s_{j}=k_{j}$ is mapped to an L-shaped domain and all other sectors are mapped to the first quadrant. For obtaining the local approximant $f_{j}(\bm{\xi}_{j})$ a similar approach as discussed for the Type 1 and 2 crease charts is followed with the exception of the non-convex sector. The local approximants $f_{j}^{s_{j}<k_{j}}(\bm{\xi}_{j})$ on each convex sector are mapped according to (24) from a tensor-product Lagrange basis $\bm{b}(\bm{\theta})$ using $\bm{\theta}=(\Lambda^{I}_{qr})^{-1}(\bm{\xi}_{j};\,4(k_{j}-1),\,s_{j})$. The local approximant on the non-convex sector $f_{j}^{s_{j}=k_{j}}(\bm{\xi}_{j})$ is composed out of three (cubic) Bézier patches that are smoothly connected within the sector while they are matched $C^{0}$ continuously with polynomial pieces in other sectors across the crease edges.

3.4 Examples of manifold surfaces with creases

We consider the construction of the manifold-based basis functions for one of the elements in the unstructured quadrilateral control mesh shown in Figure 17. The element is adjacent to a prescribed crease and belongs to four overlapping chart domains as illustrated in Figure 17. Two of the chart domains, $\hat{\Omega}_{1}$ and $\hat{\Omega}_{2}$, contain each two crease edges and the other two, $\hat{\Omega}_{3}$ and $\hat{\Omega}_{4}$, contain no crease edges. In the chart domain $\hat{\Omega}_{2}$ the arrangement of the crease edges is rotationally symmetric so that it is Type 1 and in the chart domain $\hat{\Omega}_{1}$ it is asymmetric so that it is Type 2. For the chart $\hat{\Omega}_{1}$ the mapping (25) and for the chart $\hat{\Omega}_{2}$ the mapping (19) is to be used. Hence, in comparison to the smooth case illustrated in Figure 6 only the mapping for the chart $\hat{\Omega}_{1}$ is different while the other mappings are the same. The presence of creased edges only changes the chart parametrisation for the construction of local basis $\bm{p}_{j}(\bm{\xi}_{j})$. The construction of blending functions $w_{j}(\bm{\xi}_{j})$ is the same for all charts regardless of the presence of crease edges. The tensor-product B-spline blending functions are defined in the reference element and mapped to each element in the chart domains.

In Figure 18 the application of the manifold-based basis functions in representing a smooth surface with sharp features is demonstrated. The blending functions are assembled from tensor-product cubic B-splines and the local polynomial basis is a cubic tensor-product Lagrange polynomial. The relatively coarse unstructured quadrilateral mesh has a number of tagged crease edges. The location and arrangement of creases have been chosen in order to obtain as many as possible distinct chart configurations. The eight charts with valence $v=6$ include a smooth chart, four Type 1 crease charts and three Type 2 crease charts. The rendered manifold surface $\bm{x}(\bm{\eta})$ is obtained by first projecting the control vertices $\bm{x}_{I}$ on a Catmull-Clark subdivision surface

[TABLE]

where $L_{IJ}$ is the limit matrix (or, mask) of the subdivision surface. The multiplication of the control vertices with the limit matrix involves only the one-neighbourhoods of the vertices $\bm{x}_{I}$ and projects them onto a Catmull-Clark limit surface. The entries of the limit matrix depend on the valence of the vertex and can be found, for instance, in [28]. There are also so-called tuned subdivision masks available which can provide better shapes [15, 29, 24]. Moreover, note that the Catmull-Clark subdivision surface reduces to cubic B-splines on structured meshes. With the projected control vertices the image of a reference element $\Box$ on the manifold surface is obtained as

[TABLE]

where $N_{I}(\bm{\eta})$ are the introduced manifold-based basis functions. The manifold construction ensures that the images of the reference elements on the surface are connected with the required smoothness. The projection of the control vertices onto a Catmull-Clark surface is only used to obtain a visually pleasing surface. It is inconsequential when manifold-based basis functions are used as finite element basis functions.

Finally, evidently, the geometry shown in Figure 18 can be modelled using subdivision surfaces, or even using a collection of B-spline patches, without manifold constructions. The obtained smooth manifold-based basis functions are however essential for the finite element discretisation, in particular, of higher-order partial differential equations, such as the Kirchhoff-Love thin shell equations.

4 Thin-shell formulation and discretisation

This section provides a brief review of the Kirchhoff-Love thin-shell equations and their discretisation with manifold-based basis functions. Also, the boundary and crease constraints used in the presented numerical examples are introduced. Further details of our specific Kirchhoff-Love implementation may be found in [30, 31, 32].

We denote the reference and deformed mid-surfaces of the shell with $\Omega_{0}$ and $\Omega$, respectively, and the corresponding boundaries with $\Gamma_{0}=\partial\Omega_{0}$ and $\Gamma=\partial\Omega$. Any configuration of the shell is assumed to be defined as

[TABLE]

where $\bm{\varphi}(\eta^{1},\,\eta^{2},\,\eta^{3})$ is the position vector of a material point with the convective coordinates $(\eta^{1},\,\eta^{2},\,\eta^{3})$. Similarly, $\bm{x}(\eta^{1},\,\eta^{2})$ is the position vector of a material point with the convective coordinates $(\eta^{1},\,\eta^{2})$ on the shell mid-surface. Furthermore, $\bm{a}_{3}(\eta^{1},\,\eta^{2})$ is the unit normal to the mid-surface $\bm{x}(\eta^{1},\,\eta^{2})$ and $t$ is the shell thickness. The corresponding vectors in the reference configuration are denoted with uppercase letters $\bm{\varPhi}(\eta^{1},\,\eta^{2},\,\eta^{3})$, $\bm{X}(\eta^{1},\,\eta^{2})$ and $\bm{A}_{3}(\eta^{1},\,\eta^{2})$ so that the reference configuration reads

[TABLE]

The covariant basis vectors of the tangent space of the mid-surface are defined as333Throughout Section 4, the Greek indices take the values $\{1,\,2\}$, the lowercase Latin indices take the values $\{1,\,2,\,3\}$, a comma denotes differentiation and the summation convention over repeated indices applies.

[TABLE]

with the corresponding normal vectors

[TABLE]

The contravariant basis vectors $\bm{A}^{\alpha}$ and $\bm{a}^{\alpha}$ follow from the relations

[TABLE]

where $\delta_{\alpha}^{\beta}$ is the Kronecker delta. The contravariant metric of the reference configuration needed in the following is defined as

[TABLE]

The deformation gradient can now be expressed as

[TABLE]

and a straightforward calculation yields the Green-Lagrange strain tensor

[TABLE]

with the membrane and bending strain tensors

[TABLE]

As usual, the quadratic terms in the Green-Lagrange strain tensor $\bm{E}$ are neglected for thin shells.

The potential energy of a hyper-elastic shell is given by

[TABLE]

where $W(\bm{\alpha},\,\bm{\beta})$ is the internal energy density and $\Pi^{\text{ext}}(\bm{x})$ is the potential of the external forces. The internal energy density depends, in addition to the two strain tensors $\bm{\alpha}$ and $\bm{\beta}$, on the Young’s modulus $E$, the Poisson’s ratio $\nu$ and a mainly geometry related fourth-order tensor with the components

[TABLE]

When rigid joints or clamped boundaries are present, it is necessary to constrain the surface normals. The introduced manifold-based basis functions are interpolating at the boundaries so that it is straightforward to impose displacement boundary conditions. Along the clamped boundaries $\Gamma_{0,c}$ the change of the mid-surface normal $(\bm{a}_{3}-\bm{A}_{3})$ is required to be zero. Along the joints $\Gamma_{0,j}$, the $C^{0}$ continuity of the basis functions ensures that the displacements to the left $(l)$ and to the right $(r)$ are the same. When the joints are structurally rigid the relative change of the mid-surface normals $\bm{a}_{3}^{(l)}\cdot\bm{a}_{3}^{(r)}-\bm{A}_{3}^{(l)}\cdot\bm{A}_{3}^{(r)}$ is constrained to be zero. In our current implementation the rigid joint and clamped boundary constraints are enforced with the penalty method

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are penalty parameters. It is straightforward to consider alternative constrain enforcement techniques, like the Nitsche method [33, 34] or Lagrange multipliers [31, 35].

As usual, the discrete finite element equilibrium equations are derived by first writing the potential (41) as a sum over the set of reference elements $\{\Box_{e}\}$ using the Jacobian $|\partial\bm{X}/\partial\bm{\eta}|$. Here, the index $e$ denotes the number of an element. It is emphasised that each element corresponds to a quadrilateral in the control mesh, as in isogeometric analysis with B-splines. The charts considered in the manifold construction are only used for obtaining the basis functions. With the determined basis functions the reference and deformed mid-surface of a reference element $\Box_{e}$ are approximated by

[TABLE]

where  $\bm{X}_{I}$ and $\bm{x}_{I}$ are the coordinates of the control vertex coordinates. After introducing $\bm{X}^{h}(\eta^{1},\,\eta^{2})$ and $\bm{x}^{h}(\eta^{1},\,\eta^{2})$ into the element-wise expressed potential energy (41) and numerical integration the discrete equilibrium equations follow from the stationarity principle

[TABLE]

We solve this set of nonlinear equations iteratively with the Newton-Raphson method. That is, at each iteration step $(n)$ the linear equation

[TABLE]

is solved to determine $\bm{x}_{J}^{(n)}$ for a given $\bm{x}_{J}^{(n-1)}$. The expression $\partial^{2}\Pi/\partial\bm{x}_{I}\partial\bm{x}_{J}$ on the left is the tangent stiffness matrix whereas the expression $\partial\Pi/\partial\bm{x}_{I}$ on the right is the residual vector. The full expression for (44) can be found in [4, 36].

5 Examples

We consider beam, plate and shell examples to demonstrate the utility and convergence of the proposed creased manifold-based basis functions in finite element analysis. In all the examples the Kirchhoff-Love model reviewed in Section 4 is discretised with manifold-based basis functions. The blending functions are assembled from tensor-product cubic B-splines and the local polynomial basis functions are quadratic tensor-product polynomials. In convergence studies all the finite element integrals are integrated with $9\times 9$ Gauss quadrature points. The used standard Gauss integration rules do not take into account that the introduced polynomial manifold-based basis functions consist in each element out of $4\times 4$ piecewise polynomials. For the geometrically nonlinear simulation of a pinched tube we use $3\times 3$ quadrature points to save computing time.

In constructing the basis functions all element edges on the domain boundaries are tagged as crease. This eliminates the need for ghost cells as originally used in [25]. As discussed, the number of basis functions in a chart domain has to be equal or less than the number of vertices. Hence, to be able to use a quadratic tensor-product basis in all charts, mid-face and mid-edge vertices are introduced. This increases the number of vertices in a chart domain from $2v_{j}+1$ to $6v_{j}+1$, where $v_{j}$ is the valence of the centre vertex. In all examples, the rigid joint and clamped boundary constraints are enforced with the penalty method. The penalty parameters are chosen to be $\gamma=\mathcal{O}(E)$, where $E$ is the Young’s modulus.

5.1 Propped cantilever beam

As an introductory example, we compute the deflection of a beam subjected to a uniformly distributed transversal loading, see Figure 19. To be able to use bivariate basis functions the beam is modelled as a plate with uniform width and is discretised with a structured quadrilateral mesh. The chosen boundary conditions, loading and the Poisson ratio of $\nu=0$ ensure that the structure deflects like a beam. Both the left and right ends of the beam are tagged as crease, eliminating the need for ghost cells. The left end is clamped so that the displacement and the change of the beam axis normal are required to be zero. The right end is simply supported so that only the displacement is required to be zero.

In the computations the beam midspan is modelled either as continuous, hinged or rigidly-jointed hinged. The considered meshes have always a vertex at the midspan. For the continuous beam (with no hinge) the chart domain corresponding to the midspan vertex is treated in the same way like the other chart domains. In the two cases with a hinge the midspan vertex is tagged as a crease. In geometric terms, across a hinge the displacement is continuous but its derivative may be discontinuous. And, in structural terms, across a hinge only shear forces may be transferred while the bending moment at the hinge is zero. In the case of the rigidly-jointed hinged beam, the two beam axis normals across the hinge are enforced to be weakly continuous (with the penalty method) so that it becomes possible to transfer moments.

As can be obtained by a straightforward integration the analytical deflection of the continuous propped cantilever is

[TABLE]

and the deflection of the propped cantilever with a hinge is

[TABLE]

In Figure 20 the convergence of the relative $L_{2}$ norm and energy norm errors with decreasing mesh size are plotted. As can be seen, both $L_{2}$ and energy norm errors converge in all cases with the optimal rate. The errors for the beam with the hinge are slightly smaller. As to be expected, the continuous beam and rigidly-jointed hinged beam give very similar results. This confirms the soundness of the enforcement of rigid joint constraints with the penalty method.

5.2 Square plate

Next, we consider the deflection of a square plate with different boundary conditions subjected to uniformly distributed transversal loading, see Figure 21. Two different types of boundary conditions are considered. In the first case all boundary edges are simply supported and in the second case two opposite boundary edges are simply supported and the other two are clamped. At the simply supported boundary edges the displacements are constrained to be zero and at the clamped edges both the displacements and their derivatives are constrained to be zero. Similar to the propped cantilever example also a hinge line is introduced along the centre of the plate, which is considered to be present in only some of the computations. If a hinge line is considered to be present, it is always modelled as rigidly-jointed.

The structured and unstructured coarse meshes used in the computations are depicted in Figure 22. The unstructured mesh has eight extraordinary vertices, with four vertices of valence $v_{j}=3$ and the other four of valence $v_{j}=5$. In the convergence studies the refined meshes are obtained by subdividing the shown coarse meshes with the Catmull-Clark subdivision. The analytical deflections for both boundary conditions can be found in Timoshenko and Woinowsky-Krieger [37, Chapters 5 and 6]. Two representative finite element solutions for the considered two boundary conditions are shown in Figure 23. As to be expected the deflection of the plate with two clamped edges is much smaller than the deflection of the fully simply supported plate.

In Figure 24 the convergence of the relative $L_{2}$ norm error with decreasing mesh size is plotted. In all cases the errors converge with the optimal rate of two. Figure 24a depicts the convergence for the structured mesh and Figure 24b for the unstructured mesh. Both plots contain the convergence of the fully simply supported plate as well as the partially clamped plate. In addition, the convergence of the plates with rigidly-jointed hinges are included. The continuous and rigidly-jointed plates yield the same results providing an evidence for the soundness of the proposed weak enforcement of rigid joint constraints. The errors for the less constrained fully simply supported plate are smaller than the errors for the partially clamped plate. It is noteworthy that the manifold-based basis functions can achieve the optimal convergence order even on an unstructured mesh. Moreover, the magnitude of the relative errors for the structured and unstructured meshes is almost the same.

5.3 Pinched square tube

This last example involves the geometrically nonlinear analysis of a pinched square tube subjected to two diametrically opposite concentrated forces, see Figure 25. The tube consists out of two horizontal and two vertical plates that are rigidly connected along their edges. It is discretised either with a uniform structured mesh with an element size $h=L/16$ or with an unstructured mesh shown in Figure 26. The structured and unstructured meshes have 4352 and 6400 nodes, respectively. To simulate the rigid joints between the plates the relevant edges in the mesh are tagged as a crease and the angle between the two normals across a crease are constrained to remain constant during deformation. This example is an adaptation of the pinched cylinder benchmark example widely used for comparing the performance of shell elements. The pinched square tube has also been previously considered in [31]. The deflected shapes at three different load values are shown in Figure 26. As can be inferred from the deflected shapes, the tube exhibits large membrane deformations and localised bending deformations around the rigid creases and two load application points. Furthermore, it is evident that the right angle between the plates is preserved during deformation. In the load-displacement curve depicted in Figure 27 the results for the structured and unstructured meshes are visually indistinguishable. Both are in close agreement with the reference solution obtained with the commercial finite element software Abaqus. The Abaqus result is for a structured fine mesh with 10612 nodes and a six-parameter Reissner-Mindlin type shell theory.

6 Conclusions

We have presented, by extending Majeed and Cirak [25], new manifold-based basis functions for isogeometric design and analysis of surfaces with arbitrary smoothness, prescribed sharp features and boundaries. The surface is described with a quadrilateral control mesh and $C^{0}$ continuous creases are introduced along element edges tagged by the user. Manifold techniques are extremely versatile in the sense that the chart domains, their respective transition functions, local polynomial approximants, and blending functions can all be chosen to fit the needs of the specific application at hand. We introduce creases by choosing polynomial approximants that are piecewise $C^{0}$ continuous along the element edges and by modifying the geometry of the chart domains. Due to the similarities between the $C^{0}$ continuous creases and boundaries, the new basis functions also simplify the treatment of boundaries. That is, different from [25], the boundaries can be described without introducing an additional layer of ghost elements outside the domain. The chart domains have the shape of polygonal disks with curved boundaries and are composed out of quasi-conformally mapped unit squares representing the reference finite elements. The respective transition functions are easily computed by expressing the parametric coordinates of a unit square as a complex number. Furthermore, a new type of blending function is used which is assembled from tensor-product cubic B-spline pieces. In contrast to the ones in [25], the new blending functions do not require normalisation so that the partition of unity approximation leads to non-rational basis functions. Finally, in order to obtain a mesh-based approximation scheme, the coefficients of the local polynomials are expressed as vertex coefficients using a least-squares procedure. The degree of the local polynomials has to be chosen such that there are no more coefficients than vertices in a chart domain. The number of vertices in a chart domain can be increased by introducing more vertices, i.e. refining the control mesh, while keeping the chart domain size constant. The obtained basis functions have a closed form analytic description, are locally supported and are polynomial in regular regions of the mesh.

In closing, we note that the introduced manifold-based basis function construction may be interpreted as the extension of the conventional partition of unity method [19] to manifold surfaces. As detailed in [38], depending on the choice of the local polynomials and blending functions univariate manifold-based basis functions can be made to either reproduce B-splines or are identical to B-splines. Their convergence properties can be deduced from standard partition of unity results [19, 39]. Hence, manifold-based constructions may serve as a bridge between isogeometric analysis using splines [40] and the many partition of unity method inspired discretisation techniques, such as the generalised finite element method [41], hp-clouds [42] or the extended finite element method [43]; see also the recent review [44]. Specifically, manifold-based constructions can facilitate the consideration of industrial CAD geometries, in form of NURBS and other spline representations, in the partition of unity methods. In turn, the very many enrichment techniques developed for partition of unity methods can be applied to isogeometric analysis by utilising them as local polynomial approximants in the manifold construction. The exploration of these links suggests itself as a promising direction for future research.

Appendix A Conformal maps

In this section we briefly motivate the quasi-conformal map used throughout this paper; see also [45] for a general visually orientated introduction to complex analysis. In contrast to the conformal map used in Majeed and Cirak [25] the quasi-conformal map is not angle-preserving, but it is infinitely smooth except at the extraordinary vertex. With a conformal map an infinitesimal circle is mapped to a circle whereas with a quasi-conformal map it is mapped to an ellipse. The quasi-conformal map is essential for the introduced new chart domains, which are necessary for the construction of the creased basis functions. As discussed, three different quasi-conformal maps $\bm{\varPsi}^{\textrm{I}}_{j}$, $\bm{\varPsi}^{\textrm{II}}_{j}$ and $\bm{\varPsi}^{\textrm{III}}_{j}$ are needed to map the reference element $\Box$ onto elements in the chart domains $\hat{\Omega}_{j}$ according to

[TABLE]

In the complex plane the reference element coordinates $\bm{\eta}=(\eta^{1},\,\eta^{2})$ can be, as depicted in Figure LABEL:fig:complexPlane1, expressed either in Cartesian or polar form

[TABLE]

with the radius and the phase

[TABLE]

With the Euler’s formula the polar form can be equivalently expressed by

[TABLE]

Considering a complex number as a vector from the origin to the point with the coordinates $z$, a multiplication by a scalar represents a scaling and a multiplication by $e^{i\theta}$ represents a rotation. For instance, Figure LABEL:fig:complexPlane2 illustrates that $ze^{i\theta}$ is obtained by rotating $z$ anticlockwise by an angle $\theta$ and that $1/2z$ is obtained by scaling the radius of $z$ by one half.

In the quasi-conformal map introduced in (21), repeated here for convenience,

[TABLE]

the parameter $\beta$ controls how the radius of $z$ is scaled. For $\beta=4/v_{j}$ the mapping reduces to the conformal map used in [25], i.e.,

[TABLE]

which scales the radius of $z$ to $|z|^{4/v_{j}}\,$. In the present paper we choose $\beta=1$ so that the radius is not scaled. The choice of $\beta$ on the conformal maps is illustrated in Figure 29. As visible, for $\beta=1$ the parameter lines are not orthogonal to each other within the elements and the map is not angle-preserving. Therefore, it is called a quasi-conformal map. For the introduced maps $\bm{\varPsi}^{\textrm{II}}_{j}$ and $\bm{\varPsi}^{\textrm{III}}_{j}$, i.e.  (26) and (28), it is important that the radius of $z$ is not scaled. This is necessary because each creased sector may have different number of elements. Not scaling the radius ensures that the quasi-conformal map and its derivatives across crease edges are continuous c.f. Figures 10, 11 and 14.

Appendix B Blending functions

The blending functions proposed in this paper are polynomial in contrast to the rational blending functions used in Majeed and Cirak [25]. The difference between the two constructions is best understood by comparing the two Figures 3 and 30. In Figure 30 the construction process of the rational blending functions is illustrated. Notice that the cubic B-spline basis used within the reference element is not complete so that it does not add up to one. Therefore, normalisation is necessary to satisfy the partition of unity property, resulting in rational blending functions. With the numbering introduced in Figure 30 the rational blending functions are given by

[TABLE]

It is straightforward to extend this construction to the bivariate case. As a final remark, the rational blending functions have one knot and the polynomial blending functions used in this paper have three knots within the reference element.