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.
In this section, we introduce B-spline, NURBS, and design geometrical mappings referring to  .
A knot vector is a nondecreasing sequence of real numbers in the parameter space , and the components are called knots. An open knot vector of order is a knot vector that satisfies
in which the first and the last knots are repeated.
The B-spline functions of order corresponding to the knot vector are piecewise polynomials of degree p which are constructed recursively by the Cox-de Boor recursion formula:
for For example, the piecewise quadratic polynomial B-spline functions corresponding to the open knot vector
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:
where and are control points that make B-spline functions draw a desired curve as shown in Figure 2(a).
Let be an open knot vector and let and , respectively, be the polynomial degree and order of B-spline functions . Then a B-spline surface is defined by
Figure 1. B-spline functions , of order corresponding to the knot vector .
Figure 2. (a) B-spline curve and control points on the open knot vector . (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 and 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
where the are the control points, the are the weights, , and the are the pth-degree B-spline basis functions defined on the nonperiodic (and non-uniform) knot vector
We assume that , and for all i. Setting
allows us to rewrite Equation (1) in the form
The are the rational basis functions; they are piecewise rational functions on .
2.2.2. NURBS Surface
A NURBS surface of degree in the direction and degree in the direction is a bivariate vector-valued piecewise rational function of the form
The form a bidirectional control net, the are the weights, and the and are the nonrational B-spline basis functions defined on the knot vectors
Introducing the piecewise rational basis functions
the surface Equation (4) can be written as
An example of the NURBS surface is shown in Figure 3.
2.3. Weak Solution in Sobolev Space
Let be a connected open subset of . We define the vector space
Figure 3. An example of B-spline surface with control net in three dimensional space.
to consist of all those functions which, together with all their partial derivatives of orders , are continuous on . A function is said to be a -function. If is a function defined on , we define the support of as
For an integer , we also use the usual Sobolev space denoted by . For , the norm and the semi-norm, respectively, are
Suppose we are concerned with an elliptic boundary value problem on a domain with Dirichlet boundary condition along the boundary . Let
The variational formulation of the Dirichlet boundary value problem can be written as: Find such that
where is a continuous bilinear form that is -elliptic (  ) and 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
Let , 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: such that
We can write the Galerkin form (a discrete variational equation) of (5) as follows: Given , find , where , such that
which can be rewritten as: Find the trial function such that
2.4. Variational Formulation of Equilibrium Equations of Elasticity
In elasticity, the displacement field is denoted by , and the stress field is denoted by . Let be the strain field. Then the strain-displacement and the stress-strain relations are given by
respectively, where is the differential operator matrix,
and is the symmetric positive definite matrix of material constants. Material constants are classified by the property of the material. For an isotropic elastic body,
where E is the Young’s modulus of elasticity and is Poisson’s ratio.
The equilibrium equations of elasticity are
where 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 through the relations (7). Then we consider the following system of elliptic differential equations in terms of the displacement field,
subject to the boundary conditions,
is a unit vector normal to the boundary of the domain .
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 such that on , and
The finite element approximation of the solution of (11) is to construct approximations of each component of the vector .
3. Mapping Method
We introduce a geometrical mapping from the parameter space to 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 handles effectively one-dimensional singularities. Let be an open knot vector of order . Then the B-spline functions corresponding to are
Here, , , are also called the Bernstein polynomials of degree . Let
be control points, for a constant . Then the B-spline geometrical mapping
maps the parameter space onto the physical space and its inverse is
Thus, the approximation space , where is an integer greater than or equal to and are the Bernstein polynomials (B-spline functions) of degree and contain the following singular as well as smooth functions:
In other words, the geometrical mapping is able to generate the singularity of type , where .
For example, if , then the Bernstein polynomials of degree 2 are
are control points. Then the geometrical mapping obtained by these control points and its inverse, respectively, are
Suppose where are the Bernstein polynomials corresponding the the open knot vector of order 5, then contains . Hence the approximation space for isogeometric analysis contains
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 to . To end this, we construct a NURBS surface to deal with monotone singularity of type , where q is a rational number with , is a piecewise smooth function, is the polar coordinates. the construction of the NURBS surface from 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 to the physical domain . Consider the knot vectors:
where , , , , and .
Let m and be the number of knots in the knot vectors and , respectively. Also, let k and be and , respectively. Here, if the function to be approximated has a singularity of type with , where , then the polynomial degree of B-spline functions corresponding to is .
Let , be the B-splines corresponding to the knot vector and let , be the B-splines corresponding to the knot vector . Here, the B-spline functions , , corresponding to the open knot vector are the Bernstein polynomials of degree .
Consider the control points and the weights for , , that are listed in Table 1. We now construct a NURBS surface from the parameter space onto as follows:
Here , , , are NURBS basis functions defined by
Since is the Bernstein polynomial and unless , substituting Equations (13) into (19) the NURBS surface mapping (18) becomes
Table 1. Control points and weights .
Now, we derive the derivative of the mapping by using formulas in Chapter 4.5 in  .
where is the numerator of .
Denoting , the derivative of can be expressed as
Solving the Equation (20) for , we obtain
We employ the lemma below from Chapter 3 in  in order to determine and .
Lemma 1 Let , and . Then
and the knot vector corresponding to is
Applying the lemma 1 into (21), we have
The derivative of the total weight function , also, can be described in detail by substituting the Bernstein polynomial into .
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 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 %:
· The relative error in norm in %:
· The relative error in energy norm in %:
Assuming that the Young’s modulus and the Poisson’s ratio 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,
where , and . Then, the stress field (22) satisfies the equilibrium equations of elasticity on the sector shaped domain . And the displacement field has the singularity of the form where is a smooth function.
For the design of the physical domain , we set and in the knot vector (17) so that the open knot vector corresponding to -direction is as follows:
We construct the open knot vector corresponding to -direction using the form of the knot vector in (17):
which make Bernstein polynomials in on the parameter space. We choose for control points , and set the other control points as depicted in Figure 4(a). Then the NURBS surface mapping
and the inverse of the design mapping generates the singularity of the form along the radial direction on the 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 -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 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 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 maps from the parameter space to the sector shaped domain . The coordinates of the primary control points are described. (b) Relative errors of the displacement field in the maximum norm and norm for (22) (c) Relative errors of the stress field in the maximum norm for (22). (a) The scheme of the NURBS surface that generates singular functions; (b) Relative errors of the displacement field ; (c) Relative errors of the stress field .
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 . The row “ ” 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
, where , and ’s are smooth functions.
Example 1. Let
Then is the analytic solution of the Poisson equation:
and has two singularities at and .
4.1. Patchwise NURBS Mappings and Interfaces
In Example 1, we divide the physical domain into three patches:
Let ’s are physical patches. We construct NURBS geometrical mappings and that generate singularities of the type and , respectively. They are also the design maps from the parameter space to the physical patch , for each . To build up we use the following knot vectors: For ,
where , , and . For ,
In the design mappings, we observe the following:
1) We employ control points and weights from Example 5.3 in  to build , and primary control points are shown in Figure 6(a).
2) We design the NURBS geometrical mapping that generates a singularity using the control points as depicted in Figure 6(b).
3) Using the affine transformation we define
In Figure 6, the quasi-physical patch is a physical patch translated.
4) Since a singularity does not appear in the patch , we employ the standard NURBS design technique to build the mapping from the parameter space to .
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 . 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 global basis functions crossing
Figure 6. Primary control points of the NURBS geometrical mapping and from to (a) the physical patch (b) the quasi-physical patch which has the singularity at the origin. is constructed by composition of with an affine transformation, and maps from the parameter space, to the physical space . (a) Primary control points and design of the physical patch ; (b) Primary control points and design of the quasi-physical patch .
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), ’s represent intervals corresponding to the interface on the physical domain in Figure 7(a). In , the index i means the index of the patch belonging to the interval, the other index j indicate the index of patch such that
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) 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 ’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 , we merge the nonzero basis function in and the nonzero basis function in , where such that they are reflection about the interface in the physical domain. In Figure 8(b), for example, let and be B-spline basis functions of in and , in Figure 8(b), respectively. Note that i and j in the bracket represent indices of patches. So and . Let and be B-spline basis functions of such that
· on the inter face .
We merge two B-spline basis functions and 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 , 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,
We decompose the space (28) into
The finite dimensional subspace i.e. approximation space of the Poisson equation (26) is
1) and are the number of B-spline functions in ξ- and η-direction of the patch , respectively.
2) are new global basis functions by merging two B-spline functions in and , respectively of and , respectively.
3) and are the set of B-spline basis functions composition with the inverse of NURBS surface mapping on the physical domain satisfying homogeneous boundary condition.
4) and are the set of B-spline basis functions composition with the inverse of NURBS surface mapping on the physical domain 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 and up to 14. The DOF is 3061 when . 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 be the unit disk, where ’s are minor sectors whose central angles are 120˚ for each , and
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 . The row “ ” indicates the exact values.
Then solves the following elliptic boundary value problem:
and the solution has three crack singularities at .
In Example 2, we divide the unit disk into three sectors including each crack as depicted in Figure 7(b).
Then, we build a design mapping from the parameter space to a quasi-physical sector using the proposed mapping method, for . Here, we define three physical patches by using quasi-physical sectors as follows:
which means that ’s are sectors having the same radii and the central angles as these of ’s through the transformation (30) but the position of the crack tip in is the origin other than in for . A structural drawing detailed for and is shown in Figure 10, and is designed by rotating to −90˚. Finally, we define the NURBS geometrical mapping that generates singularities as follows:
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 from to the quasi-physical patch whose crack tip is the origin, generates the singularity of the type , for . Once we design the mappings, we update them by composition with transformation so that the final NURBS geometrical mapping preserving the mapping technique in maps from the parameter space, to the physical sector , for . (a) Primary control points and design of the quasi-physical patch ; (b) Primary control points and design of the quasi-physical patch .
Figure 11. The directions of arrows represent the order of the index t of along the boundaries 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 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
1) and are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping on the physical domain satisfying the homogeneous boundary condition.
2) and are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping on the physical domain satisfying the non-homogeneous boundary condition.
Figure 9(b) and Table 5 show that the relative errors (%) in the maximum norm, the -norm, and the energy norm of the computed solution for equation (31). We can see that the relative errors in the maximum norm and -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 and -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.
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 . The row “ ” indicates the exact values.
Figure 12. (a) The graph of approximate solution of the equation (26). The DOF was used 2675 with . Under these parameters, we obtained the rel. maximum norm error 2.232E−08%, rel. 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 . Under these parameters, we obtained the rel. maximum norm error 3.022E−06%, rel. 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.
 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.
 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.
 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.
 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.
 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.
 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.
 Yan, X. (2006) Multiple Crack Fatigue Growth Modeling by Displacement Discontinuity Method with Crack-Tip Elements. Applied Mathematical Modelling, 30, 489-508.
 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.
 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.
 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.
 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.
 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.
 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.