Patchwise Mapping Method for Solving Elliptic Boundary Value Problems Containing Multiple Singularities
Author(s) Hyunju Kim
ABSTRACT
In the paper , the geometrical mapping techniques based on Non-Uniform Rational B-Spline (NURBS) were introduced to solve an elliptic boundary value problem containing a singularity. In the mapping techniques, the inverse function of the NURBS geometrical mapping generates singular functions as well as smooth functions by an unconventional choice of control points. It means that the push-forward of the NURBS geometrical mapping that generates singular functions, becomes a piecewise smooth function. However, the mapping method proposed is not able to catch singularities emerging at multiple locations in a domain. Thus, we design the geometrical mapping that generates singular functions for each singular zone in the physical domain. In the design of the geometrical mapping, we should consider the design of control points on the interface between/among patches so that global basis functions are in Cspace. Also, we modify the B-spline functions whose supports include the interface between/among them. We put the idea in practice by solving elliptic boundary value problems containing multiple singularities.

1. Introduction

It has been introduced to solve multiple crack problems by using various numerical methods. First, converting the multiple crack problems into Fredholm integral equation using two elementary solutions was introduced in  . A numerical method by using both Fredholm integral equation method and the weighted residual method was introduced in  . Numerical methods based on Galerkin approximation such as the finite element methods, boundary element methods, and meshless method were also applied to solve them  -  .

In this paper, we solve the elliptic boundary value problems with multiple singularities based on the mapping method  . But, the mapping technique proposed is not able to catch singularities emerging at multiple locations in a domain. In order to resolve the drawback, we introduced the enriched Isogeometric Analysis (IGA) in  . In the paper  , we approximate the solution on the small circular zone centered at the crack tip or point singularity by enriching the finite approximation space generated by the singular mapping introduced in the mapping method. However, it is hard to evaluate the inverse functions of the singular mapping, and the NURBS mapping so that tracking the domains of integrals whose integrand is involved both B-spline function from the singular mapping and NURBS function from the NURBS geometrical mapping, is an additional work. Also, the conditional number of the stiffness matrix could be an issue for the enriched IGA. In order to alleviate these problems, we design the geometrical map having multiple singularities by using the singular mappings only. To do that, we divide the physical domain into multiple patches which may have a singularity for each, and then design the geometrical maps by the mapping methods for each patch having a singularity. Here, we consider the design of control points on the interface between/among patches. Because this is related to the smoothness of the global basis functions. Also, we modify the B-spline functions whose supports include the interface between/among them due to the compatibility condition. In this paper, the potential of the mapping method proposed with multiple patches regarding to handling the multiple fatigue-cracks propagation in various types of plate will be shown by solving the elliptic boundary value problems with multiple singularities or cracks.

In Section 1 and 2, we briefly review definitions and terminologies that are used throughout this paper. We follow those in the book  , and we thus refer to these texts for details. And then we introduce the mapping method that generates singular functions by using B-spline or NURBS in Section 3. In Section 4, we introduced the patchwise mapping method by solving elliptic boundary value problems containing multiple singularities. Finally, the conclusions is in Section 5.

2. Nomenclature

In this section, we introduce B-spline, NURBS, and design geometrical mappings referring to  .

2.1. B-Splines

A knot vector $U=\left\{{u}_{1},{u}_{2},\cdots ,{u}_{m}\right\}$ is a nondecreasing sequence of real numbers in the parameter space $\left[a,b\right]$, and the components ${u}_{i}$ are called knots. An open knot vector of order $p+1$ is a knot vector that satisfies

${u}_{1}=\cdots ={u}_{p+1}<{u}_{p+2}\le \cdots \le {u}_{m-p-1}<{u}_{m-p}=\cdots ={u}_{m},$

in which the first and the last $p+1$ knots are repeated.

The B-spline functions ${B}_{i,k}\left(\xi \right)$ of order $k=p+1$ corresponding to the knot vector $U=\left\{{u}_{1},{u}_{2},\cdots ,{u}_{m}\right\}$ are piecewise polynomials of degree p which are constructed recursively by the Cox-de Boor recursion formula:

$\begin{array}{l}{B}_{i,1}\left(\xi \right)=\left(\begin{array}{ll}1\hfill & \text{if}\text{\hspace{0.17em}}{u}_{i}\le \xi <{u}_{i+1},\hfill \\ 0\hfill & \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{otherwise},\hfill \end{array}\\ {B}_{i,k}\left(\xi \right)=\frac{\xi -{u}_{i}}{{u}_{i+k-1}-{u}_{i}}{B}_{i,k-1}\left(\xi \right)+\frac{{u}_{i+k}-\xi }{{u}_{i+k}-{u}_{i+1}}{B}_{i+1,k-1}\left(\xi \right),\end{array}$

for $1\le i\le \left(m-k\right)$ For example, the piecewise quadratic polynomial B-spline functions ${B}_{i,5}\left(\xi \right)$ corresponding to the open knot vector

$U=\left\{0,0,0,0,0,0.15,0.5,0.75,0.9,1,1,1,1,1\right\}$

are depicted in Figure 1.

The B-spline functions are useful in design as well as finite element analysis because they have the following properties: variation diminishing, convex hull, non-negativity, piecewise polynomial, compact support, and partition of unity.

A B-spline curve is defined as follows:

$C\left(\xi \right)=\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){P}_{i},$

where $n=m-k$ and $\left\{{P}_{i}\right\}$ are control points that make B-spline functions draw a desired curve as shown in Figure 2(a).

Let ${U}_{\eta }=\left\{{v}_{1},\cdots ,{v}_{{m}^{\prime }}\right\}$ be an open knot vector and let ${p}_{\eta }$ and ${k}^{\prime }={p}_{\eta }+1$, respectively, be the polynomial degree and order of B-spline functions ${B}_{j,{k}^{\prime }}\left(\eta \right)$. Then a B-spline surface is defined by

$S\left(\xi ,\eta \right)=\underset{i=1}{\overset{n}{\sum }}\underset{j=1}{\overset{{n}^{\prime }}{\sum }}{B}_{i,k}\left(\xi \right){B}_{j,{k}^{\prime }}\left(\eta \right){P}_{i,j},$

Figure 1. B-spline functions ${B}_{i,5}\left(\xi \right)$, $i=1,2,\cdots ,9$ of order $k=5$ corresponding to the knot vector $U=\left\{0,0,0,0,0,0.15,0.5,0.75,0.9,1,1,1,1,1\right\}$.

Figure 2. (a) B-spline curve and control points on the open knot vector $\left\{0,0,0,0.15,0.42,0.76,0.76,0.91,0.91,1,1,1\right\}$. (b) B-spline basis functions corresponding to the B-spline curve shown in (a). (a) B-spline curve and control points; (b) B-spline basis functions.

where ${n}^{\prime }={m}^{\prime }-{k}^{\prime }$ and ${P}_{i,j}$ are control points that make a bidirectional control net as shown in Figure 2(b).

2.2. Nonuniform Rational B-Spline (NURBS)

In this section, we review the non-uniform rational B-splines (NURBS), NURBS curve and surface, and primary properties of NURBS.

2.2.1. NURBS Curve

A pth-degree NURBS curve is defined by

$C\left(\xi \right)=\frac{\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){w}_{i}{P}_{i}}{\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){w}_{i}},\text{ }a\le \xi \le b$ (1)

where the $\left\{{P}_{i}\right\}$ are the control points, the $\left\{{w}_{i}\right\}$ are the weights, $k=p+1$, and the $\left\{{B}_{i,k}\left(\xi \right)\right\}$ are the pth-degree B-spline basis functions defined on the nonperiodic (and non-uniform) knot vector

$U=\left\{\underset{p+1}{\underset{︸}{a,\cdots ,a}},{u}_{p+2},\cdots ,{u}_{m-k},\underset{p+1}{\underset{︸}{b,\cdots ,b}}\right\}.$

We assume that $a=0,b=1$, and ${w}_{i}>0$ for all i. Setting

${R}_{i,k}\left(\xi \right)=\frac{{B}_{i,k}\left(\xi \right){w}_{i}}{\underset{j=1}{\overset{n}{\sum }}\text{ }{B}_{j,k}\left(\xi \right){w}_{j}}$ (2)

allows us to rewrite Equation (1) in the form

$C\left(\xi \right)=\underset{i=1}{\overset{n}{\sum }}\text{ }{R}_{i,k}\left(\xi \right){P}_{i}$ (3)

The $\left\{{R}_{i,k}\left(\xi \right)\right\}$ are the rational basis functions; they are piecewise rational functions on $\xi \in \left[0,1\right]$.

2.2.2. NURBS Surface

A NURBS surface of degree ${p}_{\xi }$ in the $\xi$ direction and degree ${p}_{\eta }$ in the $\eta$ direction is a bivariate vector-valued piecewise rational function of the form

$S\left(\xi ,\eta \right)=\frac{\underset{i=1}{\overset{n}{\sum }}\underset{j=1}{\overset{{n}^{\prime }}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){B}_{j,{k}^{\prime }}\left(\eta \right){w}_{i,j}{P}_{i,j}}{\underset{i=1}{\overset{n}{\sum }}\underset{j=1}{\overset{{n}^{\prime }}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){B}_{j,{k}^{\prime }}\left(\eta \right){w}_{i,j}},\text{ }0\le \xi ,\eta \le 1$ (4)

The $\left\{{P}_{i,j}\right\}$ form a bidirectional control net, the $\left\{{w}_{i,j}\right\}$ are the weights, and the $\left\{{B}_{i,k}\left(\xi \right)\right\}$ and $\left\{{B}_{j,{k}^{\prime }}\right\}$ are the nonrational B-spline basis functions defined on the knot vectors

$U=\left\{\underset{{p}_{\xi }+1}{\underset{︸}{0,\cdots ,0}},{u}_{{p}_{\xi }+1},\cdots ,{u}_{m-\left({p}_{\xi }+1\right)},\underset{{p}_{\xi }+1}{\underset{︸}{1,\cdots ,1}}\right\},$

$V=\left\{\underset{{p}_{\eta }+1}{\underset{︸}{0,\cdots ,0}},{v}_{{p}_{\xi }+1},\cdots ,{v}_{{m}^{\prime }-\left({p}_{\eta }+1\right)},\underset{{p}_{\eta }+1}{\underset{︸}{1,\cdots ,1}}\right\}.$

Introducing the piecewise rational basis functions

${R}_{i,j}\left(\xi ,\eta \right)=\frac{{B}_{i,k}\left(\xi \right){B}_{j,{k}^{\prime }}\left(\eta \right){w}_{i,j}}{\underset{s=1}{\overset{n}{\sum }}\underset{t=1}{\overset{{n}^{\prime }}{\sum }}\text{ }{B}_{s,k}\left(\xi \right){N}_{t,{k}^{\prime }}\left(\eta \right){w}_{s,t}}$

the surface Equation (4) can be written as

$S\left(\xi ,\eta \right)=\underset{i=1}{\overset{n}{\sum }}\underset{j=0}{\overset{{n}^{\prime }}{\sum }}\text{ }{R}_{i,j}\left(\xi ,\eta \right){P}_{i,j}.$

An example of the NURBS surface is shown in Figure 3.

2.3. Weak Solution in Sobolev Space

Let $\Omega$ be a connected open subset of ${ℝ}^{d}$. We define the vector space ${\mathcal{C}}^{m}\left(\Omega \right)$

Figure 3. An example of B-spline surface with control net in three dimensional space.

to consist of all those functions $\varphi$ which, together with all their partial derivatives ${\partial }^{\alpha }\varphi \left(={\partial }_{1}^{{\alpha }_{1}}\cdots {\partial }_{d}^{{\alpha }_{d}}\varphi \right)$ of orders $|\alpha |={\alpha }_{1}+\cdots +{\alpha }_{d}\le m$, are continuous on $\Omega$. A function $\varphi \in {\mathcal{C}}^{m}\left(\Omega \right)$ is said to be a ${\mathcal{C}}^{m}$ -function. If $\Psi$ is a function defined on $\Omega$, we define the support of $\Psi$ as

$\text{supp}\Psi =\stackrel{¯}{\left\{x\in \Omega |\Psi \left(x\right)\ne 0\right\}}.$

For an integer $k\ge 0$, we also use the usual Sobolev space denoted by ${H}^{k}\left(\Omega \right)$. For $u\in {H}^{k}\left(\Omega \right)$, the norm and the semi-norm, respectively, are

$\begin{array}{l}{‖u‖}_{k,\Omega }={\left(\underset{|\alpha |\le k}{\sum }{\int }_{\Omega }{|{\partial }^{\alpha }u|}^{2}\text{d}x\right)}^{1/2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{‖u‖}_{k,\infty ,\Omega }={\mathrm{max}}_{|\alpha |\le k}\left\{\text{ess}.\mathrm{sup}|{\partial }^{\alpha }u\left(x\right)|:x\in \Omega \right\};\\ {|u|}_{k,\Omega }={\left(\underset{|\alpha |=k}{\sum }{\int }_{\Omega }{|{\partial }^{\alpha }u|}^{2}\text{d}x\right)}^{1/2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{|u|}_{k,\infty ,\Omega }={\mathrm{max}}_{|\alpha |=k}\left\{\text{ess}.\mathrm{sup}|{\partial }^{\alpha }u\left(x\right)|:x\in \Omega \right\}.\end{array}$

Suppose we are concerned with an elliptic boundary value problem on a domain $\Omega$ with Dirichlet boundary condition $g\left(x,y\right)$ along the boundary $\partial \Omega$. Let

$\mathcal{W}=\left\{w\in {H}^{1}\left(\Omega \right):{w|}_{\partial \Omega }=g\right\}\text{ }\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }\text{ }\mathcal{V}=\left\{w\in {H}^{1}\left(\Omega \right):{w|}_{\partial \Omega }=0\right\}.$

The variational formulation of the Dirichlet boundary value problem can be written as: Find $u\in \mathcal{W}$ such that

$\mathcal{B}\left(u,v\right)=\mathcal{L}\left(v\right),\text{ }\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{all}\text{\hspace{0.17em}}\text{ }v\in \mathcal{V},$ (5)

where $\mathcal{B}$ is a continuous bilinear form that is $\mathcal{V}$ -elliptic (  ) and $\mathcal{L}$ is a linear functional. The solution to (5) is called a weak solution which is equivalent to the strong (classical) solution corresponding elliptic PDE whenever u is smooth enough. The energy norm of the trial function u is defined by

${‖u‖}_{\text{eng}}={\left[\frac{1}{2}\mathcal{B}\left(u,u\right)\right]}^{1/2}.$

Let ${\mathcal{W}}^{h}\subset \mathcal{W}$, ${\mathcal{V}}^{h}\subset \mathcal{V}$ be finite dimensional subspaces. Since the NURBS basis functions do not satisfy the Kronecker delta property, in this paper we approximate the nonhomogenuous Dirichlet boundary condition by the least squares method as follows: ${g}^{h}\in {\mathcal{W}}^{h}$ such that

${\int }_{\partial \Omega }{|g-{g}^{h}|}^{2}\text{d}\gamma =\text{minimum}.$

We can write the Galerkin form (a discrete variational equation) of (5) as follows: Given ${g}^{h}$, find ${u}^{h}={w}^{h}+{g}^{h}$, where ${w}^{h}\in {\mathcal{V}}^{h}$, such that

$\mathcal{B}\left({u}^{h},{v}^{h}\right)=\mathcal{L}\left({v}^{h}\right),\text{ }\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{all}\text{\hspace{0.17em}}\text{ }{v}^{h}\in {\mathcal{V}}^{h},$

which can be rewritten as: Find the trial function ${w}^{h}\in {\mathcal{V}}^{h}$ such that

$\mathcal{B}\left({w}^{h},{v}^{h}\right)=\mathcal{L}\left({v}^{h}\right)-\mathcal{B}\left({g}^{h},{v}^{h}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{all}\text{\hspace{0.17em}}\text{test}\text{\hspace{0.17em}}\text{functions}\text{\hspace{0.17em}}{v}^{h}\in {\mathcal{V}}^{h}.$ (6)

2.4. Variational Formulation of Equilibrium Equations of Elasticity

In elasticity, the displacement field is denoted by $\left\{u\right\}={\left\{{u}_{x}\left(x,y\right),{u}_{y}\left(x,y\right)\right\}}^{\text{T}}$, and the stress field is denoted by $\left\{\sigma \right\}={\left\{{\sigma }_{x},{\sigma }_{y},{\tau }_{xy}\right\}}^{\text{T}}$. Let $\left\{\epsilon \right\}={\left\{{\epsilon }_{x},{\epsilon }_{y},{\gamma }_{xy}\right\}}^{\text{T}}$ be the strain field. Then the strain-displacement and the stress-strain relations are given by

$\left\{\epsilon \right\}=\left[D\right]\left\{u\right\},\text{ }\left\{\sigma \right\}=\left[E\right]\left\{\epsilon \right\},$ (7)

respectively, where $\left[D\right]$ is the differential operator matrix,

$\left[D\right]=\left[\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]$

and $\left[E\right]$ is the $3×3$ symmetric positive definite matrix of material constants. Material constants are classified by the property of the material. For an isotropic elastic body,

$\left[E\right]=\frac{E}{1-{\nu }^{2}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]\text{ }\text{for}\text{\hspace{0.17em}}\text{plane}\text{\hspace{0.17em}}\text{stress,}$

$\left[E\right]=\left[\begin{array}{ccc}\zeta +2\mu & \zeta & 0\\ \zeta & \zeta +2\mu & 0\\ 0& 0& \mu \end{array}\right]\text{ }\text{for}\text{\hspace{0.17em}}\text{plane}\text{\hspace{0.17em}}\text{strain}\text{.}$

Here,

$\mu =\frac{E}{2\left(1+\nu \right)},\text{ }\zeta =\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)},$

where E is the Young’s modulus of elasticity and $\nu$ $\left(0\le \nu \le 1/2\right)$ is Poisson’s ratio.

The equilibrium equations of elasticity are

${\left[D\right]}^{\text{T}}\left\{\sigma \right\}\left(x,y\right)+\left\{f\right\}\left(x,y\right)=0,\text{ }\left(x,y\right)\in \Omega ,$ (8)

where $\left\{f\right\}={\left\{{f}_{x}\left(x,y\right),{f}_{y}\left(x,y\right)\right\}}^{\text{T}}$ is the vector of internal sources representing the body force per unit volume.

The equilibrium Equation (8) can be expressed in terms of the displacement field $\left\{u\right\}$ through the relations (7). Then we consider the following system of elliptic differential equations in terms of the displacement field,

${\left[D\right]}^{\text{T}}\left[E\right]\left[D\right]\left\{u\right\}\left(x,y\right)+\left\{f\right\}\left(x,y\right)=0,\text{ }\left(x,y\right)\in \Omega ,$ (9)

subject to the boundary conditions,

$\begin{array}{l}\left[N\right]\left\{\sigma \right\}\left(s\right)=\left\{\stackrel{˜}{T}\right\}\left(s\right)=\left\{\stackrel{¯}{T}\right\}\left(s\right)={\left\{{\stackrel{¯}{T}}_{x}\left(s\right),{\stackrel{¯}{T}}_{y}\left(s\right)\right\}}^{\text{T}},\text{ }s\in {\Gamma }_{N},\\ \left\{u\right\}\left(s\right)=\left\{\stackrel{¯}{u}\right\}\left(s\right)={\left\{{\stackrel{¯}{u}}_{x}\left(s\right),{\stackrel{¯}{u}}_{y}\left(s\right)\right\}}^{\text{T}},\text{ }s\in {\Gamma }_{D},\end{array}$ (10)

where ${\Gamma }_{N}\cup {\Gamma }_{D}=\partial \Omega$,

$\left[N\right]=\left[\begin{array}{ccc}{n}_{x}& 0& {n}_{y}\\ 0& {n}_{y}& {n}_{x}\end{array}\right],$

${\left\{{n}_{x},{n}_{y}\right\}}^{\text{T}}$ is a unit vector normal to the boundary $\partial \Omega$ of the domain $\Omega$.

For the Galerkin approximation to the equilibrium equations in terms of the displacement field (9), the variational form of (9) through (10) is:

find the vector $\left\{u\right\}$ such that ${u}_{x},{u}_{y}\in {H}^{1}\left(\Omega \right),\text{\hspace{0.17em}}\left\{u\right\}=\left\{\stackrel{¯}{u}\right\}$ on ${\Gamma }_{D}$, and

$\mathcal{B}\left(\left\{u\right\},\left\{v\right\}\right)=\mathcal{F}\left(\left\{v\right\}\right),\text{ }\text{for}\text{\hspace{0.17em}}\text{all}\text{\hspace{0.17em}}\left\{v\right\}\in {H}_{0}^{1}\left(\Omega \right),$ (11)

where

$\mathcal{B}\left(\left\{u\right\},\left\{v\right\}\right)={\int }_{\Omega }{\left(\left[D\right]\left\{v\right\}\right)}^{\text{T}}\left[E\right]\left(\left[D\right]\left\{u\right\}\right)\text{d}x\text{d}y,$

$\mathcal{F}\left(\left\{v\right\}\right)={\int }_{\Omega }{\left\{v\right\}}^{\text{T}}\left\{f\right\}\text{d}x\text{d}y+{\oint }_{{\Gamma }_{N}}{\left\{v\right\}}^{\text{T}}\left\{\stackrel{¯}{T}\right\}\text{d}s$ (12)

The finite element approximation of the solution of (11) is to construct approximations of each component of the vector $\left\{u\right\}$.

3. Mapping Method

We introduce a geometrical mapping from the parameter space $\stackrel{^}{\Omega }=\left[0,1\right]×\left[0,1\right]$ to ${ℝ}^{2}$ that generates singular basis functions  .

3.1. B-Spline Curves That Generates Singular Basis Functions

In particular, we first show how a B-spline curve $F\left(\eta \right):\left[0,1\right]×\left[0,1\right]\to ℝ$ handles effectively one-dimensional singularities. Let ${U}_{\eta }=\left\{0,\cdots ,0,1,\cdots ,1\right\}$ be an open knot vector of order ${k}^{\prime }={p}_{\eta }+1$. Then the B-spline functions ${B}_{j,{k}^{\prime }}\left(\eta \right)$ corresponding to ${U}_{\eta }$ are

${B}_{j,{k}^{\prime }}\left(\eta \right)=\left(\begin{array}{c}{p}_{\eta }\\ j-1\end{array}\right){\eta }^{j-1}{\left(1-\eta \right)}^{{p}_{\eta }-j+1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}j=1,\cdots ,{k}^{\prime }.$ (13)

Here, ${B}_{j,{k}^{\prime }}$, $j=1,\cdots ,{k}^{\prime }$, are also called the Bernstein polynomials of degree ${p}_{\eta }$. Let

${P}_{j}=\left(0,0\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}j=1,\cdots ,{k}^{\prime }-1,\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{P}_{{k}^{\prime }}=\left(0,\gamma \right)$

be control points, for a constant $\gamma$. Then the B-spline geometrical mapping

$F\left(\eta \right)=\underset{j=1}{\overset{{k}^{\prime }}{\sum }}\text{ }{B}_{j,{k}^{\prime }}\left(\eta \right){B}_{j}$ (14)

$={B}_{1,{k}^{\prime }}\left(\eta \right)\left(0,0\right)+\cdots +{B}_{{p}_{\eta },{k}^{\prime }}\left(\eta \right)\left(0,0\right)+{B}_{{k}^{\prime },{k}^{\prime }}\left(\eta \right)\left(0,\gamma \right)$ (15)

$=\left(0,\gamma {\eta }^{{p}_{\eta }}\right)$ (16)

maps the parameter space $\left[0,1\right]$ onto the physical space $\left\{0\right\}×\left[0,\gamma \right]\subset {ℝ}^{2}$ and its inverse is

$\eta ={F}^{-1}\left(0,y\right)={\left(1/\gamma \right)}^{1/{p}_{\eta }}{y}^{1/{p}_{\eta }}.$

Thus, the approximation space ${\mathcal{V}}^{h}=\text{span}\left\{{B}_{i,{k}^{″}}\circ {F}^{-1}|i=1,\cdots ,{k}^{″}\right\}$, where ${k}^{″}$ is an integer greater than or equal to ${k}^{\prime }$ and ${N}_{i,{k}^{″}}$ are the Bernstein polynomials (B-spline functions) of degree ${k}^{″}-1$ and contain the following singular as well as smooth functions:

${y}^{1/{p}_{\eta }},l=0,1,\cdots ,{k}^{″}-1.$

In other words, the geometrical mapping $F$ is able to generate the singularity of type ${r}^{\lambda }$, where $0<\lambda =1/{p}_{\eta }<1$.

For example, if ${p}_{\eta }=2$, then the Bernstein polynomials of degree 2 are

${B}_{1,3}={\left(1-\eta \right)}^{2},\text{ }{B}_{2,3}=2\eta \left(1-\eta \right),\text{ }{B}_{3,3}={\eta }^{2}.$

and

${P}_{1}=\left(0,0\right),\text{ }{P}_{2}=\left(0,0\right),\text{ }{P}_{3}=\left(0,\gamma \right)$

are control points. Then the geometrical mapping obtained by these control points and its inverse, respectively, are

$F\left(\eta \right)=\left(0,\gamma {\eta }^{2}\right)\text{ }\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }\text{ }{F}^{-1}\left(0,y\right)=\sqrt{1/\gamma }\sqrt{y}.$

Suppose ${\mathcal{S}}_{\eta }^{h}=\text{span}\left\{{B}_{j,5}|j=1,\cdots ,5\right\}$ where ${B}_{j,5}$ are the Bernstein polynomials corresponding the the open knot vector $U=\left\{0,0,0,0,0,1,1,1,1,1\right\}$ of order 5, then ${\mathcal{S}}_{\eta }^{h}$ contains $1,\eta ,\cdots ,{\eta }^{4}$. Hence the approximation space ${\mathcal{V}}_{y}^{h}=\text{span}\left\{{B}_{j,5}\circ {F}^{-1}:j=1,\cdots ,5\right\}$ for isogeometric analysis contains

$1,\sqrt{y},y,{y}^{3/2},{y}^{2}.$

3.2. NURBS Surface That Generates Singular Basis Functions

The argument which is the construction of geometrical mapping that generates singular basis functions, can be extended to NURBS surface from the parameter space $\stackrel{^}{\Omega }=\left[0,1\right]×\left[0,1\right]$ to $\Omega \subset {ℝ}^{2}$. To end this, we construct a NURBS surface to deal with monotone singularity of type ${r}^{q}\psi \left(\theta \right)$, where q is a rational number with $0, $\psi \left(\theta \right)$ is a piecewise smooth function, $\left(r,\theta \right)$ is the polar coordinates. the construction of the NURBS surface from $\stackrel{^}{\Omega }$ to the unit disk in  is generalized in this section. We refer to this reference for the details.

We now consider a NURBS surface from the parameter space $\stackrel{^}{\Omega }$ to the physical domain $\Omega$. Consider the knot vectors:

${U}_{\xi }=\left\{\underset{{p}_{\xi }+1}{\underset{︸}{0,\cdots ,0}},{\zeta }_{1},\cdots ,{\zeta }_{l},\underset{{p}_{\xi }+1}{\underset{︸}{1,\cdots ,1}}\right\},\text{ }{U}_{\eta }=\left\{\underset{{p}_{\eta }+1}{\underset{︸}{0,0,\cdots ,0}},\underset{{p}_{\eta }+1}{\underset{︸}{1,1,\cdots ,1}}\right\}.$ (17)

where ${\zeta }_{i}=\left\{{\xi }_{i,1},{\xi }_{i,2},\cdots ,{\xi }_{i,{p}_{\xi }}\right\}$, ${\xi }_{i,j}={\xi }_{i,j+1}$, $i=1,\cdots ,l$, $j=1,\cdots ,{p}_{\xi }-1$, and ${\zeta }_{1}\ne {\zeta }_{2}\ne \cdots \ne {\zeta }_{l}$.

Let m and ${m}^{\prime }$ be the number of knots in the knot vectors ${U}_{\xi }$ and ${U}_{\eta }$, respectively. Also, let k and ${k}^{\prime }$ be ${p}_{\xi }+1$ and ${p}_{\eta }+1$, respectively. Here, if the function to be approximated has a singularity of type $\mathcal{O}\left({r}^{q}\right)$ with $0, where ${n}_{q},{m}_{q}\in ℤ$, then the polynomial degree of B-spline functions corresponding to ${U}_{\eta }$ is ${p}_{\eta }={m}_{q}$.

Let ${B}_{i,k}\left(\xi \right)$, $i=1,\cdots ,m-k$ be the B-splines corresponding to the knot vector ${U}_{\xi }$ and let ${B}_{j,{k}^{\prime }}\left(\eta \right)$, $j=1,\cdots ,{p}_{\eta }+1$ be the B-splines corresponding to the knot vector ${U}_{\eta }$. Here, the B-spline functions ${B}_{j,{k}^{\prime }}$, $j=1,\cdots ,{p}_{\eta }+1$, corresponding to the open knot vector ${U}_{\eta }$ are the Bernstein polynomials of degree ${p}_{\eta }$.

Consider the control points ${P}_{i,j}$ and the weights ${w}_{i,j}$ for $1\le i\le n=m-k$, $1\le j\le {p}_{\eta }+1$, that are listed in Table 1. We now construct a NURBS surface from the parameter space $\stackrel{^}{\Omega }$ onto $\Omega$ as follows:

$F\left(\xi ,\eta \right)=\underset{i=1}{\overset{n}{\sum }}\underset{j=1}{\overset{{k}^{\prime }}{\sum }}\text{ }{R}_{i,j}\left(\xi ,\eta \right){P}_{i,j}.$ (18)

Here ${R}_{i,j}\left(\xi ,\eta \right)$, $1\le i\le n$, $1\le j\le {p}_{\eta }+1$, are NURBS basis functions defined by

${R}_{i,j}\left(\xi ,\eta \right)=\frac{{B}_{i,k}\left(\xi \right){B}_{j,{k}^{\prime }}\left(\eta \right){w}_{i,j}}{W\left(\xi ,\eta \right)},$ (19)

where

$W\left(\xi ,\eta \right)=\underset{s=1}{\overset{n}{\sum }}\underset{t=1}{\overset{{k}^{\prime }}{\sum }}\text{ }{B}_{s,k}\left(\xi \right){B}_{t,{k}^{\prime }}\left(\eta \right){w}_{s,t}.$

Since ${B}_{j,{k}^{\prime }}\left(\eta \right)$ is the Bernstein polynomial and ${P}_{i,j}=\left(0,0\right)$ unless $j={k}^{\prime }$, substituting Equations (13) into (19) the NURBS surface mapping (18) becomes

$F\left(\xi ,\eta \right)={p}_{\eta }{\eta }^{{p}_{\eta }}\left[\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){w}_{i,{k}^{\prime }}{P}_{i,{k}^{\prime }}\right]/W\left(\xi ,\eta \right).$

Table 1. Control points ${P}_{i,j}$ and weights ${w}_{i,j}$.

Now, we derive the derivative of the mapping $F\left(\xi ,\eta \right)$ by using formulas in Chapter 4.5 in  .

Let

$F\left(\xi ,\eta \right)=\frac{W\left(\xi ,\eta \right)F\left(\xi ,\eta \right)}{W\left(\xi ,\eta \right)}=\frac{A\left(\xi ,\eta \right)}{W\left(\xi ,\eta \right)},$

where $A\left(\xi ,\eta \right)$ is the numerator of $F\left(\xi ,\eta \right)$.

Denoting ${\varphi }^{\left({\alpha }_{1},{\alpha }_{2}\right)}\left(\xi ,\eta \right)=\frac{{\partial }^{{\alpha }_{1}+{\alpha }_{2}}}{\partial {\xi }^{{\alpha }_{1}}\partial {\eta }^{{\alpha }_{2}}}\varphi \left(\xi ,\eta \right)$, the derivative of $F\left(\xi ,\eta \right)$ can be expressed as

${F}^{\left({\alpha }_{1},{\alpha }_{2}\right)}={\left[\frac{A\left(\xi ,\eta \right)}{W\left(\xi ,\eta \right)}\right]}^{\left({\alpha }_{1},{\alpha }_{2}\right)}={\left[\frac{W\left(\xi ,\eta \right)F\left(\xi ,\eta \right)}{W\left(\xi ,\eta \right)}\right]}^{\left({\alpha }_{1},{\alpha }_{2}\right)}$

Then

$\begin{array}{c}A{\left(\xi ,\eta \right)}^{\left({\alpha }_{1},{\alpha }_{2}\right)}={\left[W\left(\xi ,\eta \right)F\left(\xi ,\eta \right)\right]}^{\left({\alpha }_{1},{\alpha }_{2}\right)}\\ =\frac{{\partial }^{{\alpha }_{2}}}{\partial {\eta }^{{\alpha }_{2}}}\left[\underset{i=0}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right){W}^{\left(i,0\right)}{F}^{\left({\alpha }_{1}-i,0\right)}\right]\\ =\underset{i=0}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right)\underset{j=0}{\overset{{\alpha }_{2}}{\sum }}\left(\begin{array}{c}{\alpha }_{2}\\ j\end{array}\right){W}^{\left(i,j\right)}{F}^{\left({\alpha }_{1}-i,{\alpha }_{2}-j\right)}\\ ={W}^{\left(0,0\right)}{F}^{\left({\alpha }_{1},{\alpha }_{2}\right)}+\underset{i=1}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right){W}^{\left(i,0\right)}{F}^{\left({\alpha }_{1}-i,{\alpha }_{2}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\underset{j=1}{\overset{{\alpha }_{2}}{\sum }}\left(\begin{array}{c}{\alpha }_{2}\\ j\end{array}\right){W}^{\left(0,j\right)}{F}^{\left({\alpha }_{1},{\alpha }_{2}-j\right)}+\underset{i=1}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right)\underset{j=1}{\overset{{\alpha }_{2}}{\sum }}\left(\begin{array}{c}{\alpha }_{2}\\ j\end{array}\right){W}^{\left(i,j\right)}{F}^{\left({\alpha }_{1}-i,{\alpha }_{2}-j\right)}\end{array}$ (20)

Solving the Equation (20) for $F\left(\xi ,\eta \right)$, we obtain

$\begin{array}{c}F{\left(\xi ,\eta \right)}^{\left({\alpha }_{1},{\alpha }_{2}\right)}=\frac{1}{W\left(\xi ,\eta \right)}\left[{A}^{\left({\alpha }_{1},{\alpha }_{2}\right)}-\underset{i=1}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right){W}^{\left(i,0\right)}{F}^{\left({\alpha }_{1}-i,{\alpha }_{2}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\underset{j=1}{\overset{{\alpha }_{2}}{\sum }}\left(\begin{array}{c}{\alpha }_{2}\\ j\end{array}\right){W}^{\left(0,j\right)}{F}^{\left({\alpha }_{1},{\alpha }_{2}-j\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\underset{i=1}{\overset{{\alpha }_{1}}{\sum }}\left(\begin{array}{c}{\alpha }_{1}\\ i\end{array}\right)\underset{j=1}{\overset{{\alpha }_{2}}{\sum }}\left(\begin{array}{c}{\alpha }_{2}\\ j\end{array}\right){W}^{\left(i,j\right)}{F}^{\left({\alpha }_{1}-i,{\alpha }_{2}-j\right)}\right]\end{array}$ (21)

We employ the lemma below from Chapter 3 in  in order to determine $A{\left(\xi ,\eta \right)}^{\left({\alpha }_{1},{\alpha }_{2}\right)}$ and $W{\left(\xi ,\eta \right)}^{\left({\alpha }_{1},{\alpha }_{2}\right)}$.

Lemma 1 Let ${P}_{i}^{\left(0\right)}={P}_{i}$, and $C\left(\xi \right)={C}^{\left(0\right)}\left(\xi \right)=\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}\left(\xi \right){P}_{i}^{\left(0\right)}$. Then

${C}^{\left({\alpha }_{1}\right)}\left(\xi \right)=\underset{i=1}{\overset{n-{\alpha }_{1}}{\sum }}{B}_{i,k-{\alpha }_{1}}\left(\xi \right){P}_{i}^{\left(\alpha 1\right)}$

with

${P}_{i}^{\left({\alpha }_{1}\right)}=\left(\begin{array}{ll}{P}_{i},\hfill & {\alpha }_{1}=0\hfill \\ \frac{k-{\alpha }_{1}}{{u}_{i+k}-{u}_{i+{\alpha }_{1}}}\left({P}_{i+1}^{\left({\alpha }_{1}-1\right)}-{P}_{i}^{\left({\alpha }_{1}-1\right)}\right),\hfill & {\alpha }_{1}>0\hfill \end{array}$

and the knot vector corresponding to ${C}^{\left(0\right)}\left(\xi \right)$ is

${U}^{\left({\alpha }_{1}\right)}=\left\{\underset{k-{\alpha }_{1}}{\underset{︸}{0,\cdots ,0}},{u}_{k+1},\cdots ,{u}_{m-k},\underset{k-{\alpha }_{1}}{\underset{︸}{1,\cdots ,1}}\right\}.$

Applying the lemma 1 into (21), we have

${A}^{\left({\alpha }_{1},{\alpha }_{2}\right)}=\frac{{p}_{\eta }!}{\left({p}_{\eta }-{\alpha }_{2}\right)!}{\eta }^{{p}_{\eta }-{\alpha }_{2}}\left[\underset{i=1}{\overset{n-{\alpha }_{1}}{\sum }}{B}_{i,k-{\alpha }_{1}}\left(\xi \right){P}_{i,{k}^{\prime }}^{\left({\alpha }_{1}\right)}\right],$

where

${P}_{i,{k}^{\prime }}^{\left({\alpha }_{1}\right)}=\left(\begin{array}{ll}{P}_{i,{k}^{\prime }},\hfill & {\alpha }_{1}=0\hfill \\ \frac{k-{\alpha }_{1}}{{u}_{i+k}-{u}_{i+{\alpha }_{1}}}\left({P}_{i+1,{k}^{\prime }}^{\left({\alpha }_{1}-1\right)}-{P}_{i,{k}^{\prime }}^{\left({\alpha }_{1}-1\right)}\right),\hfill & {\alpha }_{1}>0.\hfill \end{array}$

The derivative of the total weight function $W\left(\xi ,\eta \right)$, also, can be described in detail by substituting the Bernstein polynomial into ${B}_{j,{k}^{\prime }}\left(\eta \right)$.

$\begin{array}{l}W{\left(\xi ,\eta \right)}^{\left({\alpha }_{1},{\alpha }_{2}\right)}\\ =\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}{\left(\xi \right)}^{\left({\alpha }_{1}\right)}\left[\underset{j=1}{\overset{{k}^{\prime }}{\sum }}\left(\begin{array}{c}{p}_{\eta }\\ j-1\end{array}\right){\left\{{\eta }^{j-1}{\left(1-\eta \right)}^{{k}^{\prime }-j}\right\}}^{\left({\alpha }_{2}\right)}{w}_{i,j}\right]\\ =\underset{i=1}{\overset{n}{\sum }}\text{ }{B}_{i,k}{\left(\xi \right)}^{\left({\alpha }_{1}\right)}\left[\frac{{p}_{\eta }!}{\left({p}_{\eta }-{\alpha }_{2}\right)!}{\left(1-\eta \right)}^{{p}_{\eta }-{\alpha }_{2}}{w}_{i,1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }+\underset{j=2}{\overset{{p}_{\eta }}{\sum }}\left(\begin{array}{c}{p}_{\eta }\\ j-1\end{array}\right){\left\{{\eta }^{j-1}{\left(1-\eta \right)}^{{k}^{\prime }-j}\right\}}^{\left({\alpha }_{2}\right)}{w}_{i,j}+\frac{{p}_{\eta }!}{\left({p}_{\eta }-{\alpha }_{2}\right)!}{\eta }^{{p}_{\eta }-{\alpha }_{2}}{w}_{i,{k}^{\prime }}\right].\end{array}$

3.3. Numerical Example of Mapping Method for Solving an Isotropic Elasticity Containing Single Singularity

The mapping method proposed was implemented in the paper  , and the paper showed that the mapping technique using NURBS geometrical mappings constructed by an unconventional choice of control points are effective for numerical solutions of elliptic boundary value problems containing a single singularity. In this subsection, we solve an elastic problem containing a singularity to show that the proposed method is also applicable to implement elastic problems.

Throughout this paper, we measure the error $\left(u-{u}^{h}\right)$ of the computed solutions obtained by isogeomtric analysis using the proposed mapping method in the following norms:

· The relative error in the maximum norm in %:

${‖u-{u}^{h}‖}_{\infty ,\text{rel}}\left(%\right)=\frac{{‖u-{u}^{h}‖}_{\infty }}{{‖u‖}_{\infty }}×100$

· The relative error in ${L}^{2}$ norm in %:

${‖u-{u}^{h}‖}_{{L}^{2},\text{rel}}\left(%\right)=\frac{{‖u-{u}^{h}‖}_{{L}^{2}}}{{‖u‖}_{{L}^{2}}}×100$

· The relative error in energy norm in %:

${‖u-{u}^{h}‖}_{\text{eng,rel}}\left(%\right)={\left[\frac{|{‖u‖}_{\text{eng}}^{2}-{‖u‖}_{\text{eng}}^{2}|}{{‖u‖}_{\text{eng}}^{2}}\right]}^{\frac{1}{2}}×100$

Assuming that the Young’s modulus $E=1000$ and the Poisson’s ratio $\nu =0.3$ in a sector of the unit circle whose the central angle is 270˚, plane strain isotropic elastic body, we consider that the following analytic stress field,

$\begin{array}{l}{\sigma }_{x}=\lambda {r}^{\lambda -1}\left[\left(2-q\left(\lambda +1\right)\right)\mathrm{cos}\left(\left(\lambda -1\right)\theta \right)-\left(\lambda -1\right)\mathrm{cos}\left(\left(\lambda -3\right)\theta \right)\right],\\ {\sigma }_{y}=\lambda {r}^{\lambda -1}\left[\left(2-q\left(\lambda +1\right)\right)\mathrm{cos}\left(\left(\lambda -1\right)\theta \right)-\left(\lambda -1\right)\mathrm{cos}\left(\left(\lambda -3\right)\theta \right)\right],\\ {\tau }_{xy}=\lambda {r}^{\lambda -1}\left[\left(\lambda -1\right)\mathrm{sin}\left(\left(\lambda -3\right)\theta \right)+q\left(\lambda +1\right)\mathrm{sin}\left(\left(\lambda -1\right)\theta \right)\right],\end{array}$ (22)

where $\lambda =2/3$, and $q=\frac{\mathrm{cos}\left[\left(\lambda -1\right)0.75\text{π}\right]}{\mathrm{cos}\left[\left(\lambda +1\right)0.75\text{π}\right]}$. Then, the stress field (22) satisfies the equilibrium equations of elasticity on the sector shaped domain ${\Omega }_{L}$. And the displacement field has the singularity of the form ${r}^{2/3}\varphi \left(\theta \right)$ where $\varphi \left(\theta \right)$ is a smooth function.

For the design of the physical domain ${\Omega }_{L}$, we set ${p}_{\xi }=2,\text{\hspace{0.17em}}{p}_{\eta }=3$ and ${\zeta }_{1}=\left\{1/3,1/3\right\},{\zeta }_{2}=\left\{2/3,2/3\right\}$ in the knot vector (17) so that the open knot vector corresponding to $\xi$ -direction is as follows:

${U}_{\xi }=\left\{0,0,0,1/3,1/3,2/3,2/3,1,1,1\right\}$ (23)

We construct the open knot vector corresponding to $\eta$ -direction using the form of the knot vector ${U}_{\eta }$ in (17):

${U}_{\eta }=\left\{0,0,0,0,1,1,1,1\right\},$ (24)

which make Bernstein polynomials in $\left[0,1\right]$ on the parameter space. We choose $\left(0,0\right)$ for control points ${P}_{i,j},\text{\hspace{0.17em}}i=1,\cdots ,k\left(=7\right),j=1,\cdots ,{p}_{\eta }\left(=3\right)$, and set the other control points as depicted in Figure 4(a). Then the NURBS surface mapping

${F}_{L}\left(\xi ,\eta \right):{\stackrel{^}{\Omega }}_{L}↦{\Omega }_{L},\text{ }{\stackrel{^}{\Omega }}_{L}=\left[0,1\right]×\left[0,1\right],$

and the inverse of the design mapping generates the singularity of the form ${r}^{1/3}\varphi \left(\theta \right)$ along the radial direction on the ${\Omega }_{L}$ in Figure 4(a).

In order to enrich the NURBS or B-spline basis functions without failing the structure of the mapping technique, we employ refinements  ,  in the NURBS functions which are used to design the physical domain as depicted in Figure 4(a). In particular, we use p-refinement to enrich the basis functions corresponding to both the open knot vectors (24) and (23). Note that inserting new knots to increase the number of basis functions along the $\eta$ -direction may cause malfunction regarding the production of singular functions  .

Figure 4(b) and Table 2 depict the relative errors of the displacement u and v in the maximum norm(blue and red line, respectively), and in the ${L}_{2}$ norm (green and purple, respectively) versus the number of degrees of freedom. Figure 4(c) and Table 3 depict that the relative errors of the stress field $\left\{\sigma \right\}$ in the maximum norm versus the number of degrees of freedom. Both Figure 4(b) and Figure 4(c) show that the proposed mapping method captures the singularity effectively, and follows the theoretical results in  . Figure 5 exhibits the relative errors of the strain energy of the stress field.

4. Patchwise NURBS Mapping Method for Solving Elliptic Boundary Value Problems Containing Multiple Singularities

In the case of that a physical domain contains multiple cracks, we re-design the

Figure 4. (a) The NURBS surface ${F}_{L}\left(\xi ,\eta \right)$ maps from the parameter space ${\stackrel{^}{\Omega }}_{L}=\left[0,1\right]×\left[0,1\right]$ to the sector shaped domain ${\Omega }_{L}$. The coordinates of the primary control points ${P}_{i,{p}_{\eta +1}},\text{\hspace{0.17em}}i=1,\cdots ,7$ are described. (b) Relative errors of the displacement field $\left\{u\right\}$ in the maximum norm and ${L}_{2}$ norm for (22) (c) Relative errors of the stress field $\left\{\sigma \right\}$ in the maximum norm for (22). (a) The scheme of the NURBS surface ${F}_{L}$ that generates singular functions; (b) Relative errors of the displacement field $\left\{u\right\}$ ; (c) Relative errors of the stress field $\left\{\sigma \right\}$.

Table 2. The relative errors (%) of the elasticity containing singularity on the sector shaped elastic material (22): relative errors (%) of displacement field.

Table 3. The relative errors (%) of the elasticity containing singularity on the sector shaped elastic material (22): The computed strain energy and the relative errors (%) of the stress field $\left\{\sigma \right\}$. The row “ $\infty$ ” indicates the exact values.

Figure 5. Relative errors of the strain energy of the stress field (22) on the sector shaped elastic body.

geometrical mapping by using both standard NURBS mappings and the proposed mappings. Then it simplifies things to describe these sub-domains by different patches. We describe how to construct the set of global basis functions crossing interfaces between patches. Throughout the following examples, we show that the patchwise mapping method is effective in dealing with a problem containing multiple singularities.

First, We apply the mapping method for the elliptic boundary value problems with multiple singularities of type

${r}_{i}^{\lambda }{\psi }_{i}\left(\theta \right)$, where ${\Omega }_{1,i}$, and ${\psi }_{i}$ ’s are smooth functions.

Example 1. Let

${u}_{1}={r}_{1}^{1/2}\mathrm{cos}\left({\theta }_{1}/2\right)+{r}_{2}^{1/2}\mathrm{sin}\left({\theta }_{2}/2\right),\text{ }f=-\Delta {u}_{1},\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}g={{u}_{1}|}_{\Gamma }$

where

$\begin{array}{l}{\Omega }_{1}=\left[-1,1\right]×\left[0,1+\sqrt{2}\right]\\ {r}_{1}=\sqrt{{x}^{2}+{y}^{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{2}=\sqrt{{\left(x-1\right)}^{2}+{\left(y-1-\sqrt{2}\right)}^{2}}\\ {\theta }_{1}={\mathrm{cos}}^{-1}\left(\frac{x}{{r}_{1}}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\theta }_{2}=-{\mathrm{cos}}^{-1}\left(\frac{x-1}{{r}_{2}}\right)\end{array}$ (25)

Then ${u}_{1}$ is the analytic solution of the Poisson equation:

$-\Delta u=f\text{\hspace{0.17em}}\text{ }\text{in}\text{\hspace{0.17em}}{\Omega }_{1}\text{ }\text{and}\text{ }u=g\text{\hspace{0.17em}}\text{ }\text{on}\text{\hspace{0.17em}}\Gamma =\partial {\Omega }_{1}$ (26)

and has two singularities at $\left(0,0\right)$ and $\left(1,1+\sqrt{2}\right)$.

4.1. Patchwise NURBS Mappings and Interfaces

In Example 1, we divide the physical domain into three patches:

$\begin{array}{l}{\Omega }_{1,1}=\left[-1,1\right]×\left[0,1\right]\\ {\Omega }_{1,2}=\left[0,1\right]×\left[1,1+\sqrt{2}\right]\\ {\Omega }_{1,3}=\left[-1,0\right]×\left[1,1+\sqrt{2}\right]\end{array}$

Let ${\Omega }_{1,i}$ ’s are physical patches. We construct NURBS geometrical mappings ${F}_{1,1}$ and ${F}_{1,2}$ that generate singularities of the type ${r}_{1}^{1/2}\psi \left({\theta }_{1}\right)$ and ${r}_{2}^{1/2}\varphi \left({\theta }_{2}\right)$, respectively. They are also the design maps from the parameter space ${\stackrel{^}{\Omega }}_{1,i}$ to the physical patch ${\Omega }_{1,i}$, for each $i=1,2$. To build up ${F}_{1,i},i=1,2$ we use the following knot vectors: For ${F}_{1,1}$,

$\begin{array}{l}{U}_{\xi }=\left\{0,0,0,{\zeta }_{1},{\zeta }_{2},{\zeta }_{3},1,1,1\right\}\\ {U}_{\eta }=\left\{0,0,0,0.5,0.5,1,1,1\right\},\end{array}$ (27)

where ${\zeta }_{1}=\left\{0.25,0.25\right\}$, ${\zeta }_{2}=\left\{0.5,0.5\right\}$, and ${\zeta }_{3}=\left\{0.75,0.75\right\}$. For ${F}_{1,2}$,

$\begin{array}{l}{U}_{\xi }=\left\{0,0,0,{\zeta }_{1},1,1,1\right\}\\ {U}_{\eta }=\left\{0,0,0,0.5,0.5,1,1,1\right\},\end{array}$

where ${\zeta }_{1}=\left\{0.5,0.5\right\}$.

In the design mappings, we observe the following:

1) We employ control points and weights from Example 5.3 in  to build ${F}_{1,1}$, and primary control points are shown in Figure 6(a).

2) We design the NURBS geometrical mapping ${\stackrel{^}{F}}_{1,2}\left(\xi ,\eta \right)$ that generates a singularity ${\left({x}^{2}+{y}^{2}\right)}^{1/4}\psi \left({\mathrm{cos}}^{-1}\left(\frac{x}{{\left({x}^{2}+{y}^{2}\right)}^{1/2}}\right)\right)$ using the control points as depicted in Figure 6(b).

3) Using the affine transformation we define

${F}_{1,2}\left(\xi ,\eta \right)\equiv {\stackrel{^}{F}}_{1,2}\left(\xi +1,\eta +\left(1+\sqrt{2}\right)\right).$

In Figure 6, the quasi-physical patch is a physical patch translated.

4) Since a singularity does not appear in the patch ${\Omega }_{1,3}$, we employ the standard NURBS design technique to build the mapping ${F}_{1,3}\left(\xi ,\eta \right)$ from the parameter space ${\stackrel{^}{\Omega }}_{1,3}$ to ${\Omega }_{1,3}$.

4.2. Construction of Global Basis Functions over Interfaces and Approximation Space

Now, we construct an approximation space by using B-spline functions which were used in the design mapping ${F}_{1,i}$. First, we consider connectivity among B-spline functions defined on different patches and are nonzero along the interface as depicted in Figure 7(a). To obtain ${C}^{0}$ global basis functions crossing

Figure 6. Primary control points of the NURBS geometrical mapping ${F}_{1,1}\left(\xi ,\eta \right)$ and ${\stackrel{^}{F}}_{1,2}\left(\xi ,\eta \right)$ from ${\stackrel{^}{\Omega }}_{1,i}$ to (a) the physical patch ${\Omega }_{1,1}$ (b) the quasi-physical patch ${Q}_{2}$ which has the singularity at the origin. ${F}_{1,2}\left(\xi ,\eta \right)$ is constructed by composition of ${\stackrel{^}{F}}_{1,2}$ with an affine transformation, and ${\stackrel{^}{F}}_{1,2}$ maps from the parameter space, ${\stackrel{^}{\Omega }}_{1,2}$ to the physical space ${\Omega }_{1,2}$. (a) Primary control points and design of the physical patch ${\Omega }_{1,1}$ ; (b) Primary control points and design of the quasi-physical patch ${Q}_{2}$.

an interface between two different patches, we merge two B-spline local basis functions defined on different patches that have the same nonzero value on the interface between the two patches. In Figure 8(a), ${I}_{i,j}$’s represent intervals corresponding to the interface on the physical domain in Figure 7(a). In ${I}_{i,j}$, the index i means the index of the patch belonging to the interval, the other index j indicate the index of patch such that

${F}_{1,j}^{-1}\left({F}_{1,i}\left({I}_{i,j}\right)\right)={I}_{j,i}.$

Figure 7. The control points for NURBS surface with (a) three patches for the Example 1, and (b) three patches for the Example 2, respectively. (a) The control points with two singularities; (b) The control points with two singularities.

Figure 8. (a) ${I}_{j,l}$ represents the interval along ξ- or η-directions corresponding to the interface in the physical domain. We follow directions of arrows when we set the global index. (b) We merge two B-spline basis functions which are interpolants for each parameter space. (a) Parameter spaces of ${F}_{1,i}$ ’s; (b) Construction of new basis function from two interpolant B-spline basis functions defined on two distinct parameter space.

To construct global basis functions which are nonzero functions on ${I}_{i,j}\cup {I}_{j,i},\text{\hspace{0.17em}}i\ne j$, we merge the nonzero basis function in ${I}_{i,j}$ and the nonzero basis function in ${I}_{j,i}$, where $i\ne j$ such that they are reflection about the interface ${F}_{1,i}\left({I}_{i,j}\right)={F}_{1,j}\left({I}_{j,i}\right)$ in the physical domain. In Figure 8(b), for example, let ${B}_{s,{k}^{\prime }}^{\left[i\right]}\left(\eta \right)$ and ${B}_{s,{k}^{\prime }}^{\left[j\right]}\left(\eta \right),\text{\hspace{0.17em}}s=1,2,3,4$ be B-spline basis functions of $\eta$ in ${\Omega }_{1,1}$ and ${\Omega }_{1,3}$, in Figure 8(b), respectively. Note that i and j in the bracket $\left[\cdot \right]$ represent indices of patches. So $i=1$ and $j=3$. Let ${B}_{t,k}^{\left[1\right]}\left(\xi \right)$ and ${B}_{t,k}^{\left[3\right]}\left(\xi \right)$ be B-spline basis functions of $\xi$ such that

${{B}_{t,k}^{\left[1\right]}|}_{\xi \in {I}_{1,3}}\ne 0\text{ }\text{when}\text{\hspace{0.17em}}t={p}_{\xi }+1,\cdots ,2\left({p}_{\xi }+1\right)-1$

${{B}_{t,k}^{\left[3\right]}|}_{\xi \in {I}_{3,1}}\ne 0\text{ }\text{when}\text{\hspace{0.17em}}t=1,\cdots ,{p}_{\xi }+1.$

Because

· ${{B}_{4,{k}^{\prime }}^{\left[1\right]}\left(\eta \right)|}_{{I}_{1,3}}={{B}_{1,{k}^{\prime }}^{\left[3\right]}\left(\eta \right)|}_{{I}_{3,1}}=1$ and

· $\left({B}_{{t}_{1},k}^{\left[1\right]}\cdot {B}_{4,{k}^{\prime }}^{\left[1\right]}\right)\circ {{F}_{1,1}^{-1}\left(x,y\right)|}_{{F}_{1,1}\left({I}_{1,3}\right)}=\left({B}_{{t}_{2},k}^{\left[3\right]}\cdot {B}_{1,{k}^{\prime }}^{\left[3\right]}\right)\circ {{F}_{1,1}^{-1}\left(x,y\right)|}_{{F}_{1,3}\left({I}_{3,1}\right)}$,

· $\left({t}_{1},{t}_{2}\right)=\left({p}_{\xi }+1,1\right),\left({p}_{\xi }+2,2\right),\cdots ,\left(2\left({p}_{\xi }+1\right)-1,{p}_{\xi }+1\right)$,

· $\left({B}_{{t}_{1},k}^{\left[1\right]}\cdot {B}_{4,{k}^{\prime }}^{\left[1\right]}\right)=\left({B}_{{t}_{2},k}^{\left[3\right]}\cdot {B}_{1,{k}^{\prime }}^{\left[3\right]}\right)$ on the inter face ${F}_{1,1}\left({I}_{1,3}\right)\left(={F}_{1,3}\left({I}_{3,1}\right)\right)$.

We merge two B-spline basis functions $\left({B}_{{t}_{1},k}^{\left[1\right]}\cdot {B}_{4,{k}^{\prime }}^{\left[1\right]}\right)$ and $\left({B}_{{t}_{2},k}^{\left[3\right]}\cdot {B}_{1,{k}^{\prime }}^{\left[3\right]}\right)$ so that we count the new function as one global basis function. The new function has a nonzero value on two distinct patches. Here, we should carefully set $\left({t}_{1},{t}_{2}\right)$, and we apply the same degree of p-refinement into each parameter space.

For a space with the non-homogeneous boundary condition in Example 1,

${\mathcal{W}}_{1}=\left\{w\left(x,y\right)\in {H}^{1}\left({\Omega }_{1}\right):{w|}_{\partial {\Omega }_{1}}=g,{\Omega }_{1}\subset {ℝ}^{2}\right\}$ (28)

We decompose the space (28) into

${\mathcal{W}}_{1,1}=\left\{w\left(x,y\right)\in {H}^{1}\left({\Omega }_{1}\right):{w|}_{\partial {\Omega }_{1}}=0,{\Omega }_{1}\subset {ℝ}^{2}\right\}$

and

${\mathcal{W}}_{1,2}=\left\{w\left(x,y\right)\in {H}^{1}\left({\Omega }_{1}\right):{w|}_{\partial {\Omega }_{1}}=g,\text{\hspace{0.17em}}{\Omega }_{1}\subset {ℝ}^{2}\right\}.$

The finite dimensional subspace i.e. approximation space of the Poisson equation (26) is

${\mathcal{W}}_{1}^{h}={\mathcal{W}}_{1,1}^{h}\oplus {\mathcal{W}}_{1,2}^{h}=\left\{{w}_{1}+{w}_{2}:{w}_{1}\in {\mathcal{W}}_{1,1}^{h},{w}_{2}\in {\mathcal{W}}_{1,2}^{h}\right\},$

${\mathcal{W}}_{1,1}^{h}\subset {\mathcal{W}}_{1,1},\text{\hspace{0.17em}}{\mathcal{W}}_{1,2}^{h}\subset {\mathcal{W}}_{1,2},$

${\mathcal{W}}_{1,i}^{h}=\text{span}\left[{\mathcal{S}}_{i,1}\cup {\mathcal{S}}_{i,1}^{new}\cup {\mathcal{S}}_{i,2}\cup {\mathcal{S}}_{i,2}^{new}\cup {\mathcal{S}}_{i,3}\right],\text{\hspace{0.17em}}i=1,2$

${\mathcal{S}}_{1,1}=\left\{\left({B}_{i,k}^{\left[1\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[1\right]}\right)\circ {F}_{1,1}^{-1}:i=2,\cdots ,{n}_{1}-1,\text{\hspace{0.17em}}j=2,\cdots ,{{n}^{\prime }}_{1}-1\right\},$

${\mathcal{S}}_{1,2}=\left\{\left({B}_{i,k}^{\left[2\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[2\right]}\right)\circ {F}_{1,2}^{-1}:i=2,\cdots ,{n}_{2}-1,\text{\hspace{0.17em}}j=2,\cdots ,{{n}^{\prime }}_{2}-1\right\},$

${\mathcal{S}}_{1,3}=\left\{\left({B}_{i,k}^{\left[3\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[3\right]}\right)\circ {F}_{1,3}^{-1}:i=2,\cdots ,{n}_{3}-1,\text{\hspace{0.17em}}j=2,\cdots ,{{n}^{\prime }}_{3}-1\right\},$

$\begin{array}{l}{\mathcal{S}}_{2,1}=\left\{\left({B}_{i,k}^{\left[1\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[1\right]}\right)\circ {F}_{1,1}^{-1}:j=1,\cdots ,{{n}^{\prime }}_{1}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}i=1,{n}_{1}\\ \text{ }\text{ }\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}i=1,\cdots ,{p}_{\xi },{n}_{1}-\left({p}_{\xi }-1\right),\cdots ,{n}_{1}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}j={{n}^{\prime }}_{1}\right\},\end{array}$ (29)

${\mathcal{S}}_{2,2}=\left\{\left({B}_{i,k}^{\left[2\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[2\right]}\right)\circ {F}_{1,2}^{-1}:j=1,\cdots ,{{n}^{\prime }}_{2}-1\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}i=1,{n}_{2}\right\},$

$\begin{array}{l}{\mathcal{S}}_{2,3}=\left\{\left({B}_{i,k}^{\left[3\right]}\cdot {B}_{j,{k}^{\prime }}^{\left[3\right]}\right)\circ {F}_{1,3}^{-1}:j=2,\cdots ,{{n}^{\prime }}_{3}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}i=1\\ \text{ }\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\text{ }}_{\text{ }}\text{and}\text{\hspace{0.17em}}i=2,\cdots ,{n}_{3}-1\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}j={{n}^{\prime }}_{3}\right\},\end{array}$

${\mathcal{S}}_{1,1}^{new}=\left\{\left({B}_{i,k}^{\left[1\right],new}\cdot {B}_{j,{k}^{\prime }}^{\left[1\right],new}\right)\circ {F}_{1,1}^{-1}:i={p}_{\xi }+2,\cdots ,{n}_{1}-\left({p}_{\xi }+1\right),\text{\hspace{0.17em}}j={{n}^{\prime }}_{1}\right\},$

${\mathcal{S}}_{1,2}^{new}=\left\{\left({B}_{i,k}^{\left[2\right],new}\cdot {B}_{j,{k}^{\prime }}^{\left[2\right],new}\right)\circ {F}_{1,2}^{-1}:i={p}_{\xi }+2,\cdots ,{n}_{2}-1,\text{\hspace{0.17em}}i\ne {p}_{\xi }+1,\text{\hspace{0.17em}}j={{n}^{\prime }}_{2}\right\},$

${\mathcal{S}}_{2,1}^{new}=\left\{\left({B}_{i,k}^{\left[1\right],new}\cdot {B}_{j,{k}^{\prime }}^{\left[1\right],new}\right)\circ {F}_{1,1}^{-1}:i={p}_{\xi }+1,{n}_{1}-{p}_{\xi },\text{\hspace{0.17em}}j={{n}^{\prime }}_{1}\right\},$

${\mathcal{S}}_{2,2}^{new}=\left\{\left({B}_{{n}_{2},k}^{\left[2\right],new}\cdot {B}_{{{n}^{\prime }}_{2},{k}^{\prime }}^{\left[2\right],new}\right)\circ {F}_{1,2}^{-1}\right\},$

where

1) ${n}_{i}$ and ${{n}^{\prime }}_{i}$ are the number of B-spline functions in ξ- and η-direction of the patch ${\Omega }_{1,i}$, respectively.

2) ${B}_{s,k}^{\left[i\right],new},{B}_{t,{k}^{\prime }}^{\left[j\right],new}$ are new global basis functions by merging two B-spline functions in ${\Omega }_{1,i}$ and ${\Omega }_{1,j}$, respectively of $\xi$ and $\eta$, respectively.

3) ${\mathcal{S}}_{1,i}$ and ${\mathcal{S}}_{1,i}^{new}$ are the set of B-spline basis functions composition with the inverse of NURBS surface mapping ${F}_{1,i}$ on the physical domain ${\Omega }_{1}$ satisfying homogeneous boundary condition.

4) ${\mathcal{S}}_{2,i}$ and ${\mathcal{S}}_{2,i}^{new}$ are the set of B-spline basis functions composition with the inverse of NURBS surface mapping ${F}_{1,i}$ on the physical domain ${\Omega }_{1}$ satisfying non-homogeneous boundary condition.

Figure 9 shows the relative errors (%) versus DOFs. In Figure 9(a) and Table 4, we enrich the set of basis functions by p-refinement and increase the degree of polynomial ${p}_{\xi }$ and ${p}_{\eta }$ up to 14. The DOF is 3061 when ${p}_{\xi }={p}_{\eta }=15$. We can see that the proposed mapping method is effective to capture multiple singularities as well as a single crack or singularity.

Example 2. Let ${\Omega }_{2}=\underset{i=1}{\overset{3}{\cup }}{\Omega }_{2,i}$ be the unit disk, where ${\Omega }_{2,i}$ ’s are minor sectors whose central angles are 120˚ for each $i=1,2,3$, and

${u}_{2}\left(x,y\right)=\underset{i=1}{\overset{3}{\sum }}{r}_{i}^{1/2}\mathrm{cos}\left({\theta }_{i}/2\right),\text{ }f=-\Delta {u}_{2}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}g={{u}_{2}|}_{\Gamma }$

Figure 9. The relative error (%) in the maximum norm, the L2-norm, and the energy norm of the computed solutions of the Poisson equation (a) (26) (b) (31) versus number of degrees of freedom, respectively. (a) Rel error vs number of degrees of freedom; (b) Rel error vs number of degrees of freedom.

Table 4. The relative errors (%) of the Poisson equation on the rectangle (26): The computed strain energy and the relative errors (%) of the approximate solution ${u}^{h}$. The row “ $\infty$ ” indicates the exact values.

where

${r}_{i}=\sqrt{{t}_{ix}^{2}+{t}_{iy}^{2}},\text{\hspace{0.17em}}\left\{\begin{array}{c}{t}_{ix}\\ {t}_{iy}\end{array}\right\}={T}_{{\alpha }_{i}}\left(\left\{\begin{array}{c}x+{x}_{i}\\ y+{y}_{i}\end{array}\right\}\right),$

${T}_{{\alpha }_{i}}\left(\left\{\begin{array}{c}x\\ y\end{array}\right\}\right)=\left[\begin{array}{cc}\mathrm{cos}\left({\alpha }_{i}\right)& -\mathrm{sin}\left({\alpha }_{i}\right)\\ \mathrm{sin}\left({\alpha }_{i}\right)& \mathrm{cos}\left({\alpha }_{i}\right)\end{array}\right]\left\{\begin{array}{c}x\\ y\end{array}\right\},\text{ }i=1,2,3$ (30)

${x}_{1}=-0.68\mathrm{cos}{12}^{\circ },\text{\hspace{0.17em}}{y}_{1}=-0.68\mathrm{sin}{12}^{\circ },\text{\hspace{0.17em}}{\alpha }_{1}={12}^{\circ }$

${x}_{2}=-0.7\mathrm{cos}{150}^{\circ },\text{\hspace{0.17em}}{y}_{2}=-0.7\mathrm{sin}{150}^{\circ },\text{\hspace{0.17em}}{\alpha }_{2}={150}^{\circ }$

${x}_{3}=0,\text{\hspace{0.17em}}{y}_{3}=0.5\mathrm{sin}{90}^{\circ },\text{\hspace{0.17em}}{\alpha }_{3}={90}^{\circ }$

${\theta }_{i}=\left(\begin{array}{ll}{\mathrm{cos}}^{-1}\left(\frac{x+{x}_{i}}{{r}_{i}}\right)\hfill & \text{if}\text{\hspace{0.17em}}y>0,\hfill \\ -{\mathrm{cos}}^{-1}\left(\frac{x+{x}_{i}}{{r}_{i}}\right)\hfill & \text{if}\text{\hspace{0.17em}}y\le 0\hfill \end{array}$

Then ${u}_{2}$ solves the following elliptic boundary value problem:

$-\Delta u=f\text{ }\text{ }\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}{\Omega }_{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u=g\text{\hspace{0.17em}}\text{ }\text{ }\text{on}\text{\hspace{0.17em}}\Gamma =\partial {\Omega }_{2},$ (31)

and the solution ${u}_{2}$ has three crack singularities at $\left({x}_{i},{y}_{i}\right),\text{\hspace{0.17em}}i=1,2,3$.

In Example 2, we divide the unit disk into three sectors ${\Omega }_{2,i},\text{\hspace{0.17em}}i=1,2,3$ including each crack as depicted in Figure 7(b).

$\left(\begin{array}{l}{\Omega }_{2,1}=\left\{\left({r}_{1},{\theta }_{1}\right)|0\le {r}_{1}\le 1,\text{\hspace{0.17em}}{150}^{\circ }\le {\theta }_{1}\le {270}^{\circ }\right\}\\ {\Omega }_{2,2}=\left\{\left({r}_{1},{\theta }_{1}\right)|0\le {r}_{1}\le 1,\text{\hspace{0.17em}}{270}^{\circ }\le {\theta }_{1}\le {390}^{\circ }\right\}\\ {\Omega }_{2,3}=\left\{\left({r}_{1},{\theta }_{1}\right)|0\le {r}_{1}\le 1,\text{\hspace{0.17em}}{30}^{\circ }\le {\theta }_{1}\le {150}^{\circ }\right\}\end{array}$

Then, we build a design mapping ${\stackrel{^}{F}}_{2,i}\left(\xi ,\eta \right)$ from the parameter space ${\stackrel{^}{\Omega }}_{2,i}$ to a quasi-physical sector ${Q}_{i}$ using the proposed mapping method, for $i=1,2,3$. Here, we define three physical patches ${\Omega }_{2,i}$ by using quasi-physical sectors ${Q}_{i}$ as follows:

${\Omega }_{2,i}=\left\{\left(x+{x}_{i},y+{y}_{i}\right)|\left(x,y\right)\in {Q}_{i}\right\},$

which means that ${Q}_{i}$ ’s are sectors having the same radii and the central angles as these of ${\Omega }_{2,i}$ ’s through the transformation (30) but the position of the crack tip in ${Q}_{i}$ is the origin other than $\left({x}_{i},{y}_{i}\right)$ in ${\Omega }_{2,i}$ for $i=1,2,3$. A structural drawing detailed for ${Q}_{1}$ and ${Q}_{2}$ is shown in Figure 10, and ${Q}_{3}$ is designed by rotating ${Q}_{1}$ to −90˚. Finally, we define the NURBS geometrical mapping that generates singularities as follows:

${F}_{2,i}\left(\xi ,\eta \right)\equiv {\stackrel{^}{F}}_{2,i}\left(\xi ,\eta \right)+\left({x}_{i},{y}_{i}\right).$

Similar to that of Example 1, considering the continuity of the basis functions and the construction of the basis functions on interfaces, we merge two basis

Figure 10. The NURBS geometrical mapping ${\stackrel{^}{F}}_{2,i}\left(\xi ,\eta \right)$ from ${\stackrel{^}{\Omega }}_{2,i}$ to the quasi-physical patch ${Q}_{i}$ whose crack tip is the origin, generates the singularity of the type ${r}_{i}^{1/2}\psi \left({\theta }_{i}\right)$, for $i=1,2,3$. Once we design the mappings, we update them by composition with transformation so that the final NURBS geometrical mapping preserving the mapping technique in ${\stackrel{^}{F}}_{2,i}$ maps from the parameter space, ${\stackrel{^}{\Omega }}_{2,i}$ to the physical sector ${\Omega }_{i}$, for $i=1,2,3$. (a) Primary control points and design of the quasi-physical patch ${Q}_{1}$ ; (b) Primary control points and design of the quasi-physical patch ${Q}_{2}$.

Figure 11. The directions of arrows represent the order of the index t of ${B}_{t,k}^{\left[i\right]}\left(\xi \right)$ along the boundaries $\eta =1$ corresponding to interfaces on the physical domain. We need to consider these directions when we set the global index. Both right and left sides colored by red of each parameter spaces are lines corresponding to each cracks in the physical domain.

functions defined on different patches that are nonzero along the interface as depicted in Figure 7(b). Figure 11 shows the order of the index t of ${B}_{t,k}^{\left[i\right]}\left(\xi \right)$ when we set the global index. Similar to the approximation space (29), the finite dimensional subspace of the weak solution of (31) has the form

${\mathcal{W}}_{2,i}^{h}=\text{span}\left[{\mathcal{S}}_{i,1}\cup {\mathcal{S}}_{i,1}^{new}\cup {\mathcal{S}}_{i,2}\cup {\mathcal{S}}_{i,2}^{new}\cup {\mathcal{S}}_{i,3}\right],\text{\hspace{0.17em}}i=1,2$

where

1) ${\mathcal{S}}_{1,i}$ and ${\mathcal{S}}_{1,i}^{new}$ are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping ${F}_{2,i}$ on the physical domain ${\Omega }_{2}$ satisfying the homogeneous boundary condition.

2) ${\mathcal{S}}_{2,i}$ and ${\mathcal{S}}_{2,i}^{new}$ are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping ${F}_{2,i}$ on the physical domain ${\Omega }_{2}$ satisfying the non-homogeneous boundary condition.

Figure 9(b) and Table 5 show that the relative errors (%) in the maximum norm, the ${L}^{2}$ -norm, and the energy norm of the computed solution for equation (31). We can see that the relative errors in the maximum norm and ${L}^{2}$ -norm decrease exponentially, and the error in the energy norm decrease almost linearly. In Figure 9(b), we enrich the basis by p-refinement along $\xi$ and $\eta$ -directions, and the number of degrees of freedom and the degree of polynomial were increased up to 14 and 4915, respectively. Figure 12(a) and Figure 12(b) show the graph of the computed solution of Equation (26) and (31), respectively.

5 Conclusions

In this paper, the physical domains of the elliptic boundary value problems containing multiple singularities, were re-designed by the patchwise mapping methods. In the patchwise mapping method, the construction of the approximation space is different from that in the conventional mapping method  due to the use of multiple singular functions.

One of the advantages of the patchwise NURBS mapping method including the NURBS mapping technique is to not only generate singular functions but also preserve the properties of IGA. The properties are the followings  :

Table 5. The relative errors (%) of the Poisson equation on the unit disk (31): The computed strain energy and the relative errors (%) of the approximate solution ${u}^{h}$. The row “ $\infty$ ” indicates the exact values.

Figure 12. (a) The graph of approximate solution of the equation (26). The DOF was used 2675 with ${p}_{\xi }={p}_{\eta }=14$. Under these parameters, we obtained the rel. maximum norm error 2.232E−08%, rel. ${L}^{2}$ norm error 4.300E−09%, and rel. energy norm error 4.330E−04%. (b) The graph of approximate solution of the equation (31). The DOF was used 4915 with ${p}_{\xi }={p}_{\eta }=14$. Under these parameters, we obtained the rel. maximum norm error 3.022E−06%, rel. ${L}^{2}$ norm error 4.563E−07%, and rel. energy norm error 9.146E−04%. (a) Approximate solution of the equation (26); (b)Approximate solution of the equation (31).

1) The capability of more precise geometric representation of complex objects than conventional Finite Element Methods.

2) Mesh refinement without altering the geometry throughout the refinement process.

Thus, we expect that the patchwise mapping method will be effective for dealing with multiple curved  , angled, branched, or radiating cracks. Also, the proposed method can be applied to compute the stress intensity factor and energy release rate in the plate theory   . These will be introduced in the subsequent paper.

On the other hand, the drawback of the mapping method is that it is not eligible to use control points and weights imported from Computer Aided Design (CAD) whereas the conventional IGA is available. To overcome this drawback, the approximation space of the standard IGA can be enriched by the mapping method to deal with singularities  . The enrichment of IGA by the mapping method is more practical because there are several advantages in the view of engineering designers and IGA users. First, the original design mapping is not needed to be changed. The enriched NURBS approximation space in IGA can generate singular functions through the proposed mapping method on a singular zone. In the mapping method, k- and h-refinement do not produce optimal results. But k-refinement is applicable in the space of NURBS basis in the enriched IGA for improved computational solution. So the enriched IGA for multiple singularities or cracks would be the future work.

Cite this paper
Kim, H. (2019) Patchwise Mapping Method for Solving Elliptic Boundary Value Problems Containing Multiple Singularities. Journal of Applied Mathematics and Physics, 7, 1572-1598. doi: 10.4236/jamp.2019.77107.
References
   Jeong, J.W., Oh, H.-S., Kang, S. and Kim, H. (2013) Mapping Techniques in Isogeometric Analysis for Elliptic Boundary Value Problems Containing Singularities. Computer Methods in Applied Mechanics and Engineering, 254, 334-352.
https://doi.org/10.1016/j.cma.2012.09.009

   Chen, Y.Z. (1984) Multiple Crack Problems of Antiplane Elasticity in an Infinite Body. Engineering Fracture Mechanics, 20, 767-775.
https://doi.org/10.1016/0013-7944(84)90085-7

   Chen, Y.Z. and Wang, Z.X. (2012) Solution of Multiple Crack Problem in a Finite Plate Using Coupled Integral Equations. International Journal of Solids and Structures, 49, 87-94.
https://doi.org/10.1016/j.ijsolstr.2011.09.015

   Barbieri, E., Petrinic, N., Meo, M. and Tagarielli, V.L. (2012) A New Weight-Function Enrichment in Meshless Methods for Multiple Cracks in Linear Elasticity. International Journal for Numerical Methods in Engineering, 90, 177-195.
https://doi.org/10.1002/nme.3313

   Budyn, E., Zi, G., Moes, N. and Belytschko, T. (2004) A Method for Multiple Crack Growth in Brittle Materials without Remeshing. International Journal for Numerical Methods in Engineering, 61, 1741-1770.
https://doi.org/10.1002/nme.1130

   Denda, M. and Dong, Y.F. (1997) Complex Variable Approach to the Bem for Multiple Crack Problems. Computer Methods in Applied Mechanics and Engineering, 141, 247-264.
https://doi.org/10.1016/S0045-7825(96)01120-6

   Mirhosseini, S. and Fariborz, S.J. (2019) Stress Analysis in a Sheet with Multiple Cracks. Applied Mathematical Modelling, 70, 31-53.
https://doi.org/10.1016/j.apm.2019.01.005

   Mousavi, S.E., Grinspun, E. and Sukumar, N. (2011) Harmonic Enrichment Functions: A Unified Treatment of Multiple, Intersecting and Branched Cracks in the Extended Finite Element Method. International Journal for Numerical Methods in Engineering, 85, 1306-1322.
https://doi.org/10.1002/nme.3020

   Yan, X. (2006) Multiple Crack Fatigue Growth Modeling by Displacement Discontinuity Method with Crack-Tip Elements. Applied Mathematical Modelling, 30, 489-508.
https://doi.org/10.1016/j.apm.2005.05.010

   Zi, G., Song, J.-H., Budyn, E., Lee, S.-H. and Belytschko, T. (2004) A Method for Growing Multiple Cracks without Remeshing and Its Application to Fatigue Crack Growth. Modelling and Simulation in Materials Science and Engineering, 12, 901-915.
https://doi.org/10.1088/0965-0393/12/5/009

   Oh, H.-S., Kim, H. and Jeong, J.W. (2014) Enriched Isogeometric Analysis of Elliptic Boundary Value Problems in Domains with Cracks and/or Corners. International Journal for Numerical Methods in Engineering, 97, 149-180.
https://doi.org/10.1002/nme.4580

   Piegl, L. and Tiller, W. (1995) The NURBS Book. Second Edition, Springer, Berlin.
https://doi.org/10.1007/978-3-642-97385-7

   Ciarlet, P.G. (1991) Basic Error Estimates for Elliptic Problems. In: Handbook of Numerical Analysis, Vol. 2, Elsevier, Amsterdam, 17-351.
https://doi.org/10.1016/S1570-8659(05)80039-0

   Cottrell, J.A., Hughes, T.J.R. and Bazilevs, Y. (2009) Isogeometric Analysis: Toward Integration of CAD and FEM.
https://doi.org/10.1002/9780470749081

   Hughes, T.J.R., Cottrell, J.A., Hughes, T.J.R. and Bazilevs, Y. (2005) Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement. Computer Methods in Applied Mechanics and Engineering, 194, 4135-4195.
https://doi.org/10.1016/j.cma.2004.10.008

   Oh, H.-S. (2003) Accurate Mode-Separated Energy Release Rate for Delamination Cracks. Journal of Computational Physics, 193, 86-114.
https://doi.org/10.1016/j.jcp.2003.07.025

   Elfakhakhre, N.R.F., Long, N.M.A. and Eshkuvatov, Z.K. (2018) Stress Intensity Factor for an Elastic Half Plane Weakened by Multiple Curved Cracks. Applied Mathematical Modelling, 60, 540-551.
https://doi.org/10.1016/j.apm.2018.03.039

   Fartash, A.H., Ayatollahi, M. and Bagheri, R. (2019) Transient Response of Dissimilar Piezoelectric Layers with Multiple Interacting Interface Cracks. Applied Mathematical Modelling, 66, 508-526.
https://doi.org/10.1016/j.apm.2018.09.030

   Wu, K.-C. (2016) Stress Intensity Factors and Energy Release Rate for Anisotropic Plates Based on the Classical Plate Theory. Composites Part B: Engineering, 97, 300-308.
https://doi.org/10.1016/j.compositesb.2016.05.011

Top