New Perspective to Isogeometric Analysis: Solving Isogeometric Analysis Problem by Fitting Load Function

1  Introduction

Isogeometric analysis (IGA) [1] is a computational mechanics technology that advances the seamless integration between computer-aided design (CAD) and computer-aided engineering (CAE). Although finite element analysis (FEA) [2] gains widespread applications in physical simulation, IGA uses the same nonuniform rational B-spline (NURBS) bases [3] as the physical domain space to construct the numerical solution domain space. Based on these NURBS basis functions, IGA can directly deal with the NURBS-based CAD model. This condition avoids tedious and time-consuming mesh generation in classical FEA when CAD model simulation is performed. Therefore, IGA saves a considerable amount of computation and greatly improves the numerical accuracy of the solution.

IGA has been initially implemented by Galerkin formulations (isogeometric Galerkin (IGA-G)) [1], thereby raising the efficiency issue due to the numerical integration. Isogeometric collocation (IGA-C) [4] focuses on solving the strong form of partial differential equations (PDEs) [5] rather than the weak form in IGA-G to promote the accuracy-to-computational-cost ratio. Substituting the derivatives of the NURBS numerical solution at some collocation points into the PDE leads to a linear system. The selection of these collocation points is the key to accuracy and convergence rate.

When PDEs are solved over a physical domain defined by a NURBS geometry, typical IGA commonly takes the NURBS as the bases of the approximation solution space, which is the same as the NURBS bases that construct the physical domain space. The NURBS bases are determined totally by the knot vector [3]. This fact indicates that the knot vectors for the NURBS bases of the solution field are directly inherited from the knot vectors of the geometry representation. However, few previous works investigated the influence of the knot vector on the approximating accuracy and convergence rate of the numerical solution. On the other hand, in the area of the data fitting in geometric design, the feature points that depict the shape and characteristics of the NURBS function to be fitted are vital to the knot vector that constructs the NURBS bases [3,6,7]. A theorem stating that the spline numerical solution and its k-order derivative have the same breakpoint sequence has been proposed and proved in [8]. Thus, we can treat the left-hand side and right-hand side of the PDE as a fitting problem, that is, approximating the load function defined by the NURBS geometry (right-hand side) with the derivatives of the NURBS numerical solution (left-hand side).

A new method of IGA is proposed in this study to improve the accuracy when PDEs are solved over the physical domain defined by NURBS geometry. We exploit the idea that the knot vector of the NURBS numerical solution is decided together by the load function and geometry representation. The feature points of the load function over the physical domain are detected in 1D, 2D, and 3D cases. Then, we merge these feature points into the initial knot vector of the spline representation of the physical domain. Given the great potential of the computational efficiency of IGA-C, we focus on the performance of the IGA-C framework with our generated knot vector. Thus, we name it as IGA-C by fitting load function (IGA-CL) method. In numerical experiments, we implement our IGA-CL framework to solve 1D, 2D, and 3D PDEs at Greville abscissae [4] and superconvergent (SC) collocation [9]. We demonstrate the efficiency of our method in terms of numerical accuracy in analysis by comparing it with the typical IGA-C method, where the knot vector of the numerical solution is the same as the knot vector of the physical domain. We also show the superiority of the knot selection scheme applied in the IGA-G framework over the typical IGA-G, which we call IGA-G by fitting load function (IGA-GL).

The rest of this paper is structured as follows: Section 1.1 provides the related work. The preliminary of the IGA is provided in Section 2. We present an overview of our IGA-CL framework in Section 3. Section 4 presents the strategy to detect feature points of load functions of the PDEs in 1D, 2D, and 3D cases. Section 4.4 introduces the method to generate the integrated knot vectors and proposes our IGA-CL framework. Section 5 conducts experiments for solving PDEs in 1D, 2D, and 3D cases and applying the IGA-G method. The comparisons exhibit the influences of numerical accuracy in terms of our knot selection. Finally, the conclusions and perspectives on future works are detailed in Section 6.

1.1 Related Work

IGA-G and applications. IGA was first introduced by [1] in 2005 to bridge the gap between CAD and CAE. The core idea of IGA is to utilize the same smooth and high-order NURBS bases for both the geometry parameterization and solution approximation. The finite element method (FEM) [2] in analysis has grown strong sufficiently. However, FEM always suffers from the time-consuming mesh generation and mesh approximation error. IGA possesses some additional merits over the traditional FEM [1,10], such as high continuity, exact representation of geometry, and low computational cost, because IGA directly takes the NURBS-expressed physical domain as the computational domain and computes on the control mesh.

The standard IGA method is implemented by Galerkin formulation (IGA-G) based on the weak form of the problem. The stiffness matrix assembly involves many numerical integrations, such as the commonly used Gauss quadrature. The integration efficiency in IGA has been first considered by [11]. Then, Gauss-Galerkin quadrature rules [12] and optimal and reduced quadrature rules [13] were developed. A promising approach for reusing the basis function evaluations in the numerical integration has been recently presented to improve the accuracy and time cost in [14]. In addition to tensor product B-splines, the isogeometric Galerkin matrices with hierarchical B-splines were efficiently computed through the interpolation, look-up and sum factorization techniques [15,16]. In the past 20 years, IGA has been widely applied in physical simulations such as structural mechanics [17–19], fluid mechanics [20–22], elastic mechanics [23–25], and shape optimization problems [26–28]. Its overview and computer implementation aspects were presented by [29].

IGA-C and applications. Although the mainstream IGA-G method based on the discretization of the weak form has the optimal convergence of 𝒪(hp+1) in the L2 norm [9], the computational efficiency relying on the integration is relatively low. Inspired by this, Auricchio et al. [4] came up with the IGA-C method. This method simply chooses some collocation points in the computational domain, substitutes the NURBS numerical solution at these collocation points into the strong form of the PDE, and finally generates a linear system. Unlike IGA-G, IGA-C considers both the superior accuracy and smoothness of NURBS basis functions and the low computational cost, which is a compromise between the accuracy and the computational time. The assembly of the stiffness matrix becomes simple and easy, and the bandwidth of the matrix becomes narrow [30]. Therefore, IGA-C is more competitive than IGA-G to some extent with respect to the accuracy-to-computational-cost ratio [30,31].

The research of IGA-C is mostly limited by the options of the ideal collocation points. In IGA-C, the commonly used collocation points are Greville abscissae [4] and Demko abscissae [30,32]. Demko abscissae are obtained from the extreme points of Chebyshev splines, and Greville abscissae are acquired from the averaging operation of the knot vector. However, their convergence orders are 𝒪(hp−1) for the odd degree p of the NURBS bases and 𝒪(hp) for the even degree p. On the basis of superconvergent theory, Anitescu et al. developed SC points [9] to advance the convergence rate. They improved the convergence order of the IGA-C solution to the optimal 𝒪(hp+1) for odd degree p. It has been proven that Cauchy-Galerkin points (CG points) exist such that the IGA-C solution can be the same as the IGA-G solution, but the exact locations of CG points are hard to determine [33]. The approximated collocation points are verified to have the convergence order of 𝒪(hp) for both odd and even p. Montardini et al. found a subset of SC points and obtained the optimal convergence order 𝒪(hp+1) for odd p [34]. Wang et al. established a superconvergent IGA-C method collocating at Greville points that utilized the recursive gradients to replace the conventional gradients of the bases [35]. Moreover, Atroshchenko et al. [36] attempted to decouple the geometry and simulation and yielded satisfactory results. Timoshenko beam and rod simulations are solved by Beirão et al. with a locking-free collocation approach [37,38]. IGA-C was also applied in fluid-structure interaction problems [39] and poromechanics problems [40]. For nearly-incompressible elasticity problems, Morganti et al. [41] employed IGA-C to achieve the same accuracy results as compressible elasticity via the mixed displacement-pressure method. An overview of IGA-C methods with convergence rates can be found in [42].

Data point sampling and fitting. In geometric design, the quality of curve/surface fitting is related to the quality of data point sampling. Efficient sampling methods enable the reconstruction of a curve/surface with a limited amount of points [43]. Filip et al. [44] estimated the number of isometric parameter points using the bounds on derivatives. Park [45] utilized an offset envelope of a line segment and presented the approximation of error-bounded polygons. This approximation plays an important role in Hausdorff distance between the control segments and the corresponding curves. Razdan [46] discussed the uniform sampling of arc length or total curvature. Then an iteration algorithm was proposed to sample arc length or weighted arc length uniformly and the bending energy [47]. Pagani et al. [43] have acquired certain sampled points in terms of the mean of the arc length and total curvature. The reconstruction of the curve and surface outperforms the uniform sampling method. Instead of the uniform sampling method, Lu et al. [48] focused on the sampling feature points on planar parameter curves to retain the shapes of the original curves and weigh the arc length and the total curvature to sample points. This method is different from uniform sampling but closely related to the feature points.

2  Preliminary

2.1 Isogeometric Galerkin Method

In CAD, the physical domain is expressed by a NURBS mapping from a parametric domain Ω^, that is, F:Ω^→Ω. The physical domain Ω can be written as

x(u)=∑i=0NNi,p(u)Pi,(1)

where {Pi}i=0,…,N are control points, and {Ni,p}i=0,…,N are the p-degree NURBS basis functions defined on the open knot vector U={u0,u1,…,uN+p+1}. The ui are called knots satisfying

0=u0=⋯=up≤up+1≤⋯≤uN+1=⋯=uN+p+1=1.

We name {ui0,ui1,…,uiK},uik<uik+1,k=0,…,K−1 the breakpoint sequence. Then, the ith NURBS basis function Ni,p(u) is defined as ([3])

Bi0(u)={1ui≤u<ui+1,0else,Bip(u)=u−uiui+p−uiBip−1(u)+ui+p+1−uui+p+1−ui+1Bi+1p−1(u),Ni,p(u)=Bip(u)ωi∑i′=0nBi′p(u)ωi′,(2)

where {ωi}i=0,…,N are weights. For high dimensions, the basis functions are calculated by their tensor product. In particular, when all the weights are 1, {Ni,p}i=0,…,N are B-spline basis functions {Bip}i=0,…,N.

In particular, let the following Laplace equation over the physical domain defined by a NURBS geometry Ω be

ΔT+f=0inΩ,T=gonΓD,∇T⋅n=honΓN,βT+∇T⋅n=ronΓR,(3)

where the Dirichlet boundary ΓD, Neumann boundary ΓN, and Robin boundary ΓR form the boundary ∂Ω. f, g, h, r are all given functions, and T is the unknown function. We call Eq. (3) the strong form of the PDE.

The IGA-G method ([1]) focuses on the weak form of the PDE, which is given by

∫Ω∇ω⋅∇TdΩ+β∫ΓRωTdΓ=∫ΩωfΩ+β∫ΓNωhdΓ+β∫ΓRωrdΓ,(4)

where the weighting function ω∈𝒱,𝒱={ω|ω∈H1(Ω),ω|γD=0}, the trial solution T∈𝒮,𝒮={T|T∈H1(Ω),T|γD=g}, and H1 is a Sobolev space. The weak form (4) is then rewritten by a bilinear operator a(⋅,⋅) and a linear form L(⋅) as

a(ω,T)=L(ω).(5)

Galerkin’s method aims to construct the finite-dimensional approximations of the spaces 𝒮,𝒱, i.e., 𝒮h∈𝒮,𝒱h∈𝒱. In the IGA-G method, these spaces can be spanned by the same NURBS bases as the physical domain Ω. Let the unknown numerical solution Tr be a linear combination of these basis functions with N+1 unknown coefficients, that is,

Tr(u)=∑i=0NNi,p(u)ti.(6)

Then the weak form can be transformed into

a(ωr,Tr)=L(ωr).(7)

By assembling the system of the weak form (7), we finally solve the following linear equations and obtain the solution Tr:

Kt=b.(8)

2.2 Isogeometric Collocation Method

The IGA-C method ([4]) is also proposed to solve the PDE. In general, the boundary value problem (BVP) can be summarized as follows:

{LT=finΩ,GT=gon∂Ω,(9)

where L represents a scalar differential operator, T is the unknown solution, and the given right-hand-side f is called the load function. The boundary vector operator GT represents the boundary conditions.

IGA-C solves the equations by introducing a set of collocation points {u¯k}k=0,…,Nint in the parametric domain Ω^ and {u¯k}k=Nint+1,…,Nint+Nbd (Nint+Nbd≥N) on the boundary of parametric domain ∂Ω^. The corresponding points of physical domain are x¯k=F(u¯k),k=0,…,Nint+Nbd. Thus, replacing the unknown T of the BVP with the expression of Tr at {x¯k}i=0,…,Nint+Nbd leads to a linear system:

{LTr(x¯k)=f(x¯k)k=0,…,Nint,GTr(x¯k)=g(x¯k)k=Nint+1,…,Nint+Nbd.(10)

When Nint+Nbd>N, the linear system can be solved in least-squares fashion. In this study, we only discuss the common case, where Nint+Nbd=N, which has a unique solution of {ti}i=0,…,N. The choice of the collocation points is the crucial part of the procedure, affecting directly the accuracy and convergence of the numerical solution. Some possible options exist, such as Greville abscissae ([4]), SC points ([9]), CG points ([33]), and clustered SC points ([34]). The collocation points of these schemes are generated based on the knot vector U. For instance, the typical selection, Greville abscissae in 1D, determines the collocation points as u¯i=1p(ui+ui+1+⋯+ui+p−1),i=1,…,N+1. The SC collocation points are generated proportionally from SC points of the spline basis of the degree p on the interval [−1,1] to the knot interval [ui,ui+1],i=p,…,N−p−1.

3  Algorithm Overview

The IGA-G and IGA-C methods involve the derivatives of the NURBS basis functions when the numerical solution Tr is brought into the PDE. Without loss of generality, a 1D second-order differential equation −d2T/dx2=f(x) is taken as an example. On the one hand, NURBS basis functions and their derivatives have the same breakpoint sequence; thus, the breakpoints of Tr should be the same as the breakpoints of −d2Tr/dx2. For instance, Fig. 1 demonstrates that a cubic NURBS basis and its derivatives share the same breakpoints. On the other hand, we can treat the left-hand side and right-hand side of the PDE as a fitting problem. In particular, we approximate the load function f(x) on the NURBS-expressed physical domain Ω with the derivatives of the NURBS numerical solution −d2Tr/dx2. In NURBS curve/surface fitting referring to [3], the knot vector for the basis should apply the technique of averaging, which is generated by the averages of the parameters assigned to the data points.

Figure 1: The cubic NURBS basis Ni,3 and its scaled first-order and scaled second-order derivatives Ni,3′,Ni,3″. The basis functions and their derivatives have the same breakpoints ui,ui+1,ui+2,ui+3,ui+4

Therefore, the feature points on the load function f(x) are essential to the knot vector for the left-hand-side fitted function and for the numerical solution Tr. The NURBS basis functions constructing the numerical solution are decided by the knot vector; thus, utilizing the information of feature points on the load function to generate the knot vector is helpful in improving the computational accuracy.

Given the great potential of computational efficiency of IGA-C, this study focuses on the performance of IGA-C framework with our generated knot vector. We propose the IGA-CL method, which combines the representation of the physical domain and the load function characteristics to construct the knot vector for the numerical solution in the IGA-C framework. The main idea is illustrated in Fig. 2. First, uniform parameterization of the physical domain defined by the NURBS geometry is done in the parametric domain to obtain a set of parameters, and the corresponding coordinates in the physical domain are calculated by NURBS mapping. Second, via feature point detection operation of the load function presented in Section 4, we choose the parameters assigned to these feature points, which serve as the knots to insert into the initial knot vector for the physical domain. Finally, on the basis of the integrated knot vector, we construct the NURBS bases for the numerical solution that are different from the initial bases for the physical domain. We employ the IGA-C framework to obtain the numerical solution to the BVP, where both the basis functions and collocation points depend on our integrated knot vector.

Figure 2: The main idea of IGA-CL: The integrated knot vector combines the initial knot vector of the physical domain Ω and the feature points of the load function f(Ω). Then, it serves as the knot vector for the NURBS numerical solution Tr. The dotted box on the right presents the detailed flow of the dotted box on the left showing how to obtain the knots

4  Generation of Knot Vectors by Detecting Feature Points of Load Functions

In this section, we introduce the detecting method of feature points from the load functions of the PDEs defined on a physical domain represented by NURBS geometry. We develop a method of knot generation to construct the knot vector once the feature points are detected.

4.1 One-Dimensional Detection

We begin with the detecting method of the load function of the PDE in a 1D case. Comparing with uniform sampling methods using arc length and curvature information ([43,47]), Lu et al. [48] devised a nonuniform point sampling method for fitting the parametric curves with B-spline curves. This method captures the salient features on the curve and produces noticeable fitting results. It begins with the uniform parameterization of the parametric curve; thus, we employ the uniform parameterization to our NURBS-expressed physical domain. Therefore, we obtain the coordinates in the physical domain by NURBS mapping corresponding to the generated parameters and the load function values at these coordinates. Suppose the load function f(x),x∈Ω=[a,b], where

x=∑i′=0NPi′Ni′,p(u),u∈[0,1],

and Ni′,p is NURBS basis defined on the initial knot vector U, then the parametric curve can be written as R(u)=(x(u),f(x(u))),u∈Ω^=[0,1]. Seed points are acquired, including endpoints, inflection points, and local extreme curvature points. Auxiliary points should be added to generate the detected points until the required number is satisfied, which is realized by a characteristic function, thereby better capturing the feature of the curve.

4.1.1 Seed Points

First, the physical domain Ω=[a,b] is recognized as an expression of one-to-one NURBS mapping from the parametric domain Ω^=[0,1]. We parameterize the domain [0, 1] uniformly to obtain M small intervals and M+1 parameters u^i,i=0,…,M. Thus, the corresponding M+1 coordinates x^i in domain [a, b] are calculated by x^i=∑i′=0NPi′Ni′,p(u^i), where a=x^0≤x^1≤…≤x^M=b. Then, the point set of the curve at these coordinates is denoted as {R(u^i)=(x^i,f(x^i))}i=0,1,…,M. This M is expected to be sufficiently large so that the interior of each small interval has at most one seed point. The curve curvature is denoted as k^i:=k(u^i), defined as

k(u)= det(R′(u),R′′(u))‖R′(u)‖3= f′′(x)(1+f′2(x))3/2.(11)

The seed points include endpoints, inflection points, and extreme curvature points, which are determined by referring to [48]. In the present study, we briefly summarize the rules of determination as follows:

Inflection points If k^i=0 and k^i−1k^i+1≤0, then point R(u^i) is an inflection point. If k^i≠0 and k^ik^i+1≤0, then an inflection point exists in the interior of the interval (u^i,u^i+1). We denote the parameter of this inflection point as u¯i. Let u^i+0.5=12(u^i+u^i+1) and k^i+0.5:=k(u^i+0.5) be the curvature on u^i+0.5. When k^i+0.5 is extremely close to zero, we let u¯i=u^i+0.5=12(u^i+u^i+1). Otherwise, the parabolic function k^(u)=a2u2+a1u+a0 is considered to interpolate points (u^i,k^i), (u^i+0.5,k^i+0.5), and (u^i+1,k^i+1). The coefficients of the parabolic function are determined uniquely on the basis of the interpolation conditions. As long as we let u¯i∈(u^i,u^i+1) be the root of function k^(u), R(u¯i) is an approximating inflection point. Then, u¯i is adopted as a detected point. Note that in our implementation, when |k^i|<10−5, k^i is taken as 0.

Extreme curvature points If the conditions of k^i−1<k^i and k^i+1<k^i inside the interval (u^i−1,u^i+1) are satisfied, then a local maximum curvature point must lie on it. On the contrary, if k^i−1>k^i and k^i+1>k^i, then a local minimum curvature point exists in (u^i−1,u^i+1). R(u¯i) is denoted as the extreme curvature point. We seek u¯i through parabolic interpolation method, where k^(u)=b2u2+b1u+b0 is used to interpolate points (u^i−1,k^i−1), (u^i,k^i), and (u^i+1,k^i+1). Let k^′(u)=0, and we obtain the following:

u¯i=−b12b2=(u^i2−u^i−12)(k^i−k^i+1)+(u^i+12−u^i2)(k^i−k^i−1)2(u^i−u^i−1)(k^i−k^i+1)+2(u^i+1−u^i)(k^i−k^i−1).(12)

Then u¯i is adopted as a detected point.

The seed point detection is concluded in Algorithm 1. As illustrated, the domain [0, 1] of the parametric curve of the load function R(u)=(x(u),f(x(u))) is first divided into small intervals, so a set of points in the physical domain [a, b] are generated. Given the calculation of the curvatures at these points, the next step is to push the two endpoints into the initial seed point set. Then, we loop all the curvatures to determine the inflection and extreme curvature points according to the schemes above and push them to the current seed point set. Finally, we obtain the corresponding detected point set {u¯i}i=0,…,m.

:

4.1.2 Auxiliary Points

The seed point set {R(u¯i)}i=0,…,m of curve R(u) is assumed to be already generated, including two endpoints, inflection points, and extreme curvature points. To achieve better detection performance, we need more points to depict the characteristics of the load function.

Define λ(u) as the characteristic function of curve R(u):

λ(u)=(1−γ)L(u)L(1)+γK(u)K(1),u∈[0,1],(13)

where γ is the weight, and L(u) and K(u) are the arc length and total curvature from R(0) to R(u), respectively.

L(u)=∫0u‖R′(s)‖ds,K(u)=∫0u|k(s)|dL(s)=∫0u|k(s)|‖R′(s)‖ds.(14)

Apparently, λ(u) is a strictly increasing function, where λ(0)=0 and λ(1)=1 hold naturally. The integrals Eq. (14) can be numerically computed by discretization. For instance, L(u)=∑i=0NL−1∫s1si+1||R′(s)||ds, where si∈[0,uNL,2uNL,…,u]. For this matter, Gauss numerical integration is commonly considered an advanced choice to solve the inner ∫sisi+1‖R′(s)‖ds because of the difficulty in computing the integral. The property of the characteristic function depends only on weight γ.

The feature point detection is presented in Algorithm 2. Algorithm 1 returns m+1 seed points as the initial feature point set. Then the parametric curve is initially split into m curve segments divided by the acquired seed points. We determine the salient curve segment by calculating the difference between the characteristic function values of each segment at the two endpoints, i.e., λ(u¯i+1)−λ(u¯i),i=0,…,m−1, to obtain auxiliary detected points. Here, salient refers to the curve segment with the largest difference among all curve segments (denoted as j-th curve segment). We seek the u¯ that reaches the middle of this salient curve segment, i.e., λmid(u¯)=12(λ(u¯j)+λ(u¯j+1)). By employing the bisection method through setting a threshold ε, we find the middle sampling point u¯mid such that the characteristic function value at this point is nearest the median characteristic value of the salient curve segment, i.e., (λ(u¯mid)−λmid(u¯))<ε. Then the u¯mid is the approximated u¯ that we adopt as a detected point. The search for auxiliary points is an iteration process until the required number of detected points is satisfied.

4.2 Two-Dimensional Detection

2D case is regarded as the extension of the 1D case along the u− and v−directions. [43] presented the mixed parameterization weighting arc length and curvature to the reconstruction of curves and surfaces similar to [47] rather than the square of the curvature. We exploit the surface sampling of [43] to detect the feature points of the load functions of PDEs here.

A load function f(x, y) with parametric surface is written as R(u,v)=(x(u,v),y(u,v),f(x(u,v),y(u,v))), (x,y)∈Ω, where Ω is the physical domain represented by NURBS surface

S(u,v)=∑i=0N1∑j=0N2Pi,jNi,j;p,q(u,v),u,v∈[0,1],

where Ni,j;p,q(u,v)=Bip(u)Bjq(v)ωij∑i′=0N1∑j′=0N2Bj′p(u)Bj′q(v)ωi′j′. Detection is conducted in u-direction and v-direction because NURBS surfaces are tensor product surfaces. In this study, the seed point set of surfaces here only consists of endpoints. We detect feature points on the parametric surface mainly by constructing the characteristic function. The characteristic function contains the terms of area and curvature. We define the area of the parametric surface as below:

A=∫01∫01‖Ru(u,v)×Rv(u,v)‖dudv.(15)

Ru(u,v) and Rv(u,v) are partial derivatives in Eq. (15). Thus, the marginal cumulative areas along u- and v-directions are computed as

Au(u)=∫0uds∫01‖Rs(s,v)×Rv(s,v)‖dv,Av(v)=∫0vds∫01‖Ru(u,s)×Rs(u,s)‖du.(16)

The mean curvature of the two principle curvatures is good for reflecting the surface shape. Nevertheless, we maintain the curvature computed in each direction without loss of generality. Fix v=v0 and u=u0, respectively, and we obtain the curvatures in u- and v-directions:

ku(u,v0)=|f″(u,v0)|(1+f′2(u,v0))32,kv(u0,v)=|f″(u0,v)|(1+f′2(u0,v))32.(17)

The curvature of the load function f at point (u0,v0) takes the average of the two directional curvatures:

k(u0,v0)=ku(u0,v0)+kv(u0,v0)2.(18)

Subsequently, we express the marginal cumulative integrals of curvature along the u- and v-directions:

ku(u)=∫0uds∫01k(s,v)‖Rs(s,v)×Rv(s,v)‖dv,kv(v)=∫0vds∫01k(u,s)‖Ru(u,s)×Rs(u,s)‖du.(19)

Thus, the marginal cumulative characteristic functions along the u- and v-directions are constructed as follows:

λu(u)=(1−γu)Au(u)Au(1)+γuku(u)ku(1),u∈[0,1],λv(v)=(1−γv)Av(v)Av(1)+γvkv(v)kv(1),v∈[0,1].(20)

Based on λu(u), we generate the detected points {u¯i}i=0,…,n1 in u−direction in accordance with Algorithm 2. Based on λv(v), we also attain the detected points {v¯j}j=0,…,n2 in v−direction.

:

4.3 Three-Dimensional Detection

For the 3D case, we construct the characteristic functions in a way similar to that in the 2D case. The load function is f(x, y, z) and the parametric solid is R(u,v,w)=(x(u,v,w),y(u,v,w),z(u,v,w),f(x,y,z)), (x,y,z)∈Ω, where Ω is the physical domain defined by

S(u,v,w)=∑i=0N1∑j=0N2∑k=0N3Pi,j,kNi,j,k;p,q,r(u,v,w),u,v,w∈[0,1],

where Ni,j,k;p,q,r(u,v,w)=Bip(u)Bjq(v)Bkr(w)ωijk∑i′=0N1∑j′=0N2∑k′=0N3Bi′p(u)Bj′q(v)Bk′r(w)ωi′j′k′. Detected points are computed from u-, v-, and w-directions.

The volume of the parametric solid is presented as below:

V=∫01∫01∫01|(Ru×Rv)⋅Rw|dudvdw.(21)

Thus, the marginal cumulative volumes in three directions can be described as

Vu(u)=∫0uds∫01∫01|(Rs×Rv)⋅Rw|dvdw,Vv(v)=∫0vds∫01∫01|(Ru×Rs)⋅Rw|dudw,Vw(w)=∫0wds∫01∫01|(Ru×Rv)⋅Rs|dudv.(22)

On the other hand, the curvatures in u-, v-, and w-directions are

ku(u,v0,w0)=|f″(u,v0,w0)|(1+f′(u,v0,w0)2)32,kv(u0,v,w0)=|f″(u0,v,w0)|(1+f′(u0,v,w0)2)32,kw(u0,v0,w)=|f″(u0,v0,w)|(1+f′(u0,v0,w)2)32.(23)

Then, the curvature of the load function f at point (u0,v0,w0) is

k(u0,v0,w0)=ku(u0,v0,w0)+kv(u0,v0,w0)+kw(u0,v0,w0)3.(24)

The marginal cumulative integrals of the curvature in u-, v-, and w-directions are

ku(u)=∫0uds∫01∫01k(s,v,w)|(Rs×Rv)⋅Rw|dvdw,kv(v)=∫0vds∫01∫01k(u,s,w)|(Ru×Rs)⋅Rw|dudw,kw(w)=∫0wds∫01∫01k(u,v,s)|(Ru×Rv)⋅Rs|dudv.(25)

Therefore, the marginal cumulative characteristic functions of the load function f in three directions are constructed as follows:

λu(u)=(1−γu)Vu(u)Vu(1)+γuku(u)ku(1),u∈[0,1],λv(v)=(1−γv)Vv(v)Vv(1)+γvkv(v)kv(1),v∈[0,1],λw(w)=(1−γw)Vw(w)Vw(1)+γwkw(w)kw(1),w∈[0,1].(26)

Given that the process is similar to surface detection, we omit the details here.

4.4 Generation of Knot Vector

Many numerical approaches were proposed to solve PDEs, among which IGA-C is an efficient method, which is described in Section 2.2. However, the knot vector of the IGA-C method can influence not only the NURBS bases for the numerical solution but also the collocation points.

Given the 1D scenario, the other two cases are treated similarly. Since the point detection is done, we obtain the feature point set {R(u¯i)}i=0,…,n and the corresponding parameter set {u¯i}i=0,…,n. Determining the collocation points is vital to solve a PDE in IGA-C. To make good use of feature points that depict the shape and characteristics of the load function, the parameter set {u¯i}i=0,…,n of the feature points will be served as knots to be inserted into the knot vector U for the NURBS bases of the physical domain. Then, the knot vector U∗ for the NURBS bases of the solution field that we constructed through insertion can be deemed as the derivation of the initial knot vector U of Ω. Hence, assuming the initial open knot vector of the physical domain represented by NURBS is

U={u0,u1,…,uN+p+1},(27)

we generate additional knots by averaging operation ([3]) based on Schoenberg-Whitney theorem ([49]) from the feature points for each direction in curve/surface fitting:

uj+p∗=1p∑i=jj+p−1u¯i,j=1,2,…,n−p.(28)

Moreover, Hosseini et al. summarized three knot placement techniques for geometric construction from the data points to decide the positions of the internal knots of the knot vector ([50]), namely, uniform knot placement, De Boor’s algorithm, and Piegl and Tiller’s algorithm. Here, the generated knots representing the characteristics of the load function are inserted into the knot vector U. The knots {uj+p∗}j=1,…,n−p that are extremely close to the internal knots {ui}i=p+1,…,N in U should be abandoned, because the multiplicity of the knot leads to unnecessary decrease of the continuity here, and the integrated knot vector should be sorted again. These processes are needed to ensure the nondecreasing knot vector and maintain the continuity of the splines. We denote the nondecreasing integrated knot vector as

U∗={u0,u1,…,uN∗+p+1}.(29)

The spline numerical solution is then expressed as Tr(u)=∑0NtiNi,p∗(u) to connect it with IGA-C and solve the PDE.

Employing collocation schemes such as Greville abscissae, superconvergent points, and Cauchy-Galerkin points yields a set of collocation points. Replacing the unknown function T(u) in the PDE with the numerical solution Tr(u) at collocation points leads to a linear system. Such a knot vector constructed for numerical solution combines the information of the physical domain and the load function. We name it as IGA-CL. It reflects how feature points are distributed, and distinguished from the knot vector grabbed directly from the knot vector of geometry expression in typical IGA-C. We only need to calculate the knot vectors U∗ and V∗ from u- and v-directions for the 2D cases and the knot vector W∗ from the additional w-direction for the 3D cases.

5  Experimental Results

In this section, we validate the influence of computing accuracy with our isogeometric analysis by fitting load functions in solving PDEs through numerical experiments. The typical IGA-C and our IGA-CL are compared in Sections 5.1–5.3. We verify that the IGA-CL method outperforms IGA-C by contrasting relative errors and absolute errors, particularly for the case when the solution has obvious characteristics. In addition, we extend to exploit the IGA-GL method, which displays the same superiority of the accuracy over the typical IGA-G method through the comparisons in Section 5.4. For the following examples, relative error er is used, which is defined as

er=∫Ω(T−Tr)t(T−Tr)dΩ∫ΩTtTdΩ,(30)

where t stands for the transpose, the analytical solution and numerical solution are denoted as T and Tr, respectively. When the solution is a scalar-valued function, er=∫Ω(T−Tr)2dΩ∫ΩT2dΩ.

5.1 One-Dimensional PDE

Example I: A 1D Poisson problem with Dirichlet boundary condition:

{−d2Tdx2=f1(x)inΩ=[0,2],T|∂Ω(x)=T1(x)on∂Ω,(31)

where

f1(x)=13500(2((x−0.3)2+0.02)2−2(2x−0.6)2((x−0.3)2+0.02)3)

and the analytical solution is

T1(x)=13500((x−0.3)2+0.02).

Given the advantage of our IGA-CL method, we utilize the feature points of the load function and construct an integrated knot vector. To verify the accuracy improvement, we apply our IGA-CL framework to solve Eq. (31) and compare it with the typical IGA-C at Greville abscissae ([4]) and other collocation schemes such as SC points ([9]) of the same degrees of freedom.

Fig. 3a displays the detecting points of the load function f1(x) according to the point detection. The weight in the characteristic function is set to 0.1, and the number of auxiliary points is set to 7. The knot vector is then constructed with our knot generation method. The NUBRS bases are calculated iteratively based on the knot vector. Fig. 3b shows an instance of cubic NUBRS bases from our integrated knot vector. This instance of cubic NURBS bases is distinguished from the NURBS bases generated only using the knot vector of the physical domain. The distribution of these bases appears dense in the area where the load function is steep; otherwise, the distribution of these bases becomes sparse.

Figure 3: Our detected points on the load function and the corresponding cubic NURBS bases with knot generation method

Fig. 4 shows the comparisons of convergence of the relative error by exploiting our IGA-CL framework and the typical IGA-C framework at Greville abscissae with the same degrees of freedom under different degrees of the bases. The convergence of solving this 1D problem by IGA-CL with the integrated knot vector is expressed by solid lines, whereas the convergence by typical IGA-C is displayed by dashed lines. We denote the abscissa as h∗=log10⁡(1/dofs), where dofs is the abbreviation of degrees of freedom, which are obtained by h-refinements that enrich the basis without changing the physical domain geometrically or parametrically. In our implementation, if the original knot vector is U∗={u0,u1,⋯,uN∗+p+1} where u0=u1=⋯=up=0 and uN∗+1=uN∗+2=⋯=uN∗+p+1=1, after once h-refinement, then we obtain the refined knot vector U1∗={u0,u1,…,up,12(up+up+1),up+1,12(up+1+up+2),up+2,…,uN∗,12(uN∗+uN∗+1),uN∗+1,uN∗+2,…,uN∗+p+1}. In other words, we refer to each “h-refinement” as inserting the midpoint of each nonzero knot interval into the current knot vector each time.

Figure 4: Convergence comparison of IGA-CL/IGA-C at Greville collocation points under different degrees of Example I

Let the ordinate be the log10⁡(er). We describe the convergence rates with the help of triangles, on which the numbers stand for the convergence rates. Besides Greville abscissae, SC collocation points are also applied to IGA-CL in Fig. 5. Notice that the number of SC collocation points is almost twice the degrees of freedom, we solve it in the way of least-squares.

Figure 5: Convergence comparison of IGA-CL/IGA-C at SC points under different degrees of Example I

Our method demonstrates its superiority over typical IGA-C in terms of accuracy under all the degrees p of the NURBS basis functions. Because of the knot selection from the load function, the relative errors of small dofs at the beginning of IGA-CL are much lower than those of the typical IGA-C that uses the knot vector only from the physical domain. The general convergence rates of IGA-CL remain unchanged or slightly lower than those of IGA-C. For Greville collocation points, the convergence orders of IGA-C and IGA-CL are both 𝒪(hp) for even p and 𝒪(hp−1) for odd p. For SC collocation points, the convergence orders of IGA-CL and IGA-C are both 𝒪(hp) for even p; however, the convergence order of IGA-CL is slightly lower than 𝒪(hp+1) for a high odd p. Given the desirable accuracy, IGA-CL still outperforms the IGA-C.

For PDEs under other boundary conditions, we present the following example:

Example II: A 1D ODE with Dirichlet and Neumann boundary conditions:

{−d2Tdx2+T=f2(x)inΩ=[0,1],T(0)=T2(0),T′(1)=T2′(1).(32)

where

f2(x)=(101−10000(x−0.3)2)e−(x−0.3)20.02

and the analytical solution is

T2(x)=e−(x−0.3)20.02.

In this case, we set the parameter γ as 0.1 and the number of auxiliary points as 4. The detection as well as the convergence of relative error at Greville are displayed in Fig. 6.

Figure 6: Detected points on the load function and the convergence results of relative error

5.2 Two-Dimensional PDE

We compare IGA-CL with IGA-C at Greville abscissae and SC points in solving the following two 2D Examples III and IV. The main idea of the IGA-CL is detecting the feature points on the load function and then linking them to IGA-C for solving PDEs. Whether the accuracy of the IGA-CL solution outperforms IGA-C principally depends on the quality of the feature points detected. The weights in the characteristic functions can be controlled and adjusted, thereby influencing the detection effect. We show the convergence of examples with a relatively satisfying choice of weight.

Example III: A 2D Poisson problem on a quarter of an annulus Ω (Fig. 7a) with Dirichlet boundary condition is given as follows:

−ΔT=f3(x,y)inΩ,T|∂Ω(x,y)=T3(x,y)on∂Ω.(33)

Figure 7: Example III. The analytical solution and the load function have nearly the same characteristics

The load function is

f3(x,y)=π2sin⁡(πy)A+4sin⁡(πy)+2πcos⁡(πy)(2y−5)A2−2sin⁡(πy)((2x−2)2+(2y−5)2)A3+(π2−2−(2x−2)2)sin⁡(πy)e(x−1)2−5

where A=(x−1)2+(y−2.5)2+0.06 such that the analytical solution is

T3(x,y)=(e(x−1)2−5+1A)sin⁡(πy).

Example IV: An elliptic PDE on a quarter of a plate with a circular hole Ω (Fig. 8a) with Dirichlet condition is given as follows:

{−ΔT+T=f4(x,y)inΩ,T|∂Ω(x,y)=T4(x,y)on∂Ω.(34)

Figure 8: Example IV. The analytical solution and the load function have nearly the same characteristics

The load function is

f4(x,y)=B[sin⁡(πx)(π2+41−(20x+20)2−(20y+20)2)+2πcos⁡(πx)(20x+20)]+C[sin⁡(πx)(−41−π2+(20x+50)2+(20y+70)2)−2πcos⁡(πx)(20x+50)]

where B=e−10(x+1)2−10(y+1)2, C=e−10(x+2.5)2−10(y+3.5)2 and

T4(x,y)=(B−C)sin⁡(πx).

To begin with, we present Figs. 7 and 8 displaying the analytical solutions and the load functions to demonstrate the effectiveness of our method. For the convenience of comparing errors, the analytical solutions for all experiments in this study are assumed to be known. In fact, in “real settings”, the “exact” solutions could be computed with a dense mesh. It is observed that the load functions roughly retain the characteristics of the analytical solutions of the PDEs. Consequently, detecting the feature points of the load function is a helpful guidance when solving PDEs.

In order to solve the elliptic equations Examples III and IV, the IGA-CL and IGA-C methods at Greville and SC collocation points are employed and compared in terms of absolute error. The initial number of detected points is 15×15 and 17×15 for Examples III and IV, respectively. Fig. 9 displays the comparison of IGA-C solution and IGA-CL solution of Example III, while Fig. 10 displays the comparison of Example IV. Both figures are shown under the degrees of p=q=6 at Greville. The selection of the weight γ in the characteristic function is vital for detecting. Based on various experiments, if the load function has a curvature distribution that is mostly the same or varies slowly, a large weight γ is suggested; otherwise, a small weight γ is suggested.

Figure 9: Comparison of IGA-CL/IGA-C solution of Example III after three times of h-refinement

Figure 10: Comparison of IGA-CL/IGA-C solution of Example IV after three times of h-refinement

However, in these two cases, considering the extreme nonuniformity after h-refinements when a small weight is adopted, we take a large weight in practical experiments to ensure that the flat part of the shape of load function has fairly and evenly detected points. For simplicity, the weights we adopted here are γu=γv=0.8 in the two directions for Example III and γu=γv=0.95 for both directions for Example IV.

Figs. 9d and 10d show that the distribution of the collocation points in the physical domain is very related to the solution. This finding indicates that the initial number of the detected points adopted in the area where the load function reaches the local maximum or the minimum is more than the number of detected points in the area where the load function appears steady. Through the comparisons of absolute error between IGA-C and IGA-CL methods, our IGA-CL method has a much lower error than IGA-C method.

With the appropriate weights, we exhibit the convergence comparisons of the relative error of IGA-CL and IGA-C at Greville abscissae and SC points under different degrees of NURBS bases of the numerical solution of Examples III and IV in Figs. 11 and 12, respectively. Our IGA-CL outperforms IGA-C distinctly under the same dofs. The convergence orders of the IGA-CL solution appear a little lower than those of IGA-C because the difference between the size of knot intervals of IGA-C and IGA-CL diminishes as the dofs increase. However, IGA-CL solutions still show an obvious increase in accuracy, which is conspicuously considerable.

Figure 11: Convergence comparison of IGA-CL/IGA-C at Greville/SC collocation under different degrees of Example III

Figure 12: Convergence comparison of IGA-CL/IGA-C at Greville/SC collocation under different degrees of Example IV

5.3 Three-Dimensional PDE

Example V: A 3D elliptic PDE, on a unit cube with Dirichlet boundary condition is given as follows:

{−ΔT+T=f5(x,y,z)inΩ,T|∂Ω(x,y,z)=T5(x,y,z)on∂Ω,(35)

where

f5(x,y,z)=(12π2−4(x+2)2−1)e(x+2)2sin⁡(2πx)sin⁡(2πy)sin⁡(2πz)−4π(2x+4)e(x+2)2cos⁡(2πx)sin⁡(2πy)sin⁡(2πz),

and the analytical solution is

T5(x,y,z)=e(x+2)2sin⁡(2πx)sin⁡(2πy)sin⁡(2πz).

In this case, the initial number of feature points that we sampled in our IGA-CL method is 6 for each direction, and the weights in characteristic functions are γu=γv=γw=0.1. Fig. 13c demonstrates the numerical solution of our IGA-CL method for solving Example V after two times of h-refinements when the degrees are p=q=r=4. These figures are visualized on ParaView 5.9.0-RC3.

Figure 13: IGA-CL solution of Example V after two times of h-refinements

Then, we compare the distribution of absolute errors between IGA-CL and IGA-C at Greville abscissae with the same dofs at an arbitrary degree of the numerical solution after two times of h-refinements, as shown in Fig. 14. The absolute errors of IGA-CL and IGA-C methods are shown when p=q=r=4, with the same dofs of 20×20×20. The maximum absolute error of IGA-C solving the PDE is 0.4390, and the relative error is 1.3465×10−4. However, the maximum absolute error of our IGA-CL method is 0.0775, and the relative error is 2.1933×10−5.