Isogeometric Least-squares Collocation Method with Consistency and Convergence Analysis

1.2 MB·

by  Hongwei Lin

Isogeometric Least-squares Collocation Method with Consistency
and Convergence Analysis
HongweiLina,∗,YunyangXionga,QianqianHub
aDepartmentofMathematics,StateKeyLab.ofCAD&CG,ZhejiangUniversity,Hangzhou,310027,China
bDepartmentofMathematics,ZhejiangGongshangUniversity,Hangzhou,310018,China
6
1
0
Abstract
2

n Isogeometric analysis (IGA) approximates the unknown solution of a boundary value problem
a (orinitialvalueproblem)usinganon-uniformrationalbasisspline(NURBS)functionwithnon-
J
linearNURBSbasisfunctions.IftheorderoftheNURBSfunctionishighenough,theboundary
7
valueproblemcanbesolvedusinganisogeometriccollocation(IGA-C)methodthatinterpolates
2
the real differential operator at collocation points. In this paper, we present the isogeometric

] least-squarescollocation(IGA-L)method. IGA-Ldeterminesthenumericalsolutionbymaking
A the approximate differential operator fit the real differential operator in a least-squares sense.
N Therefore, while IGA-C requires the number of collocation points to be equal to that of the
. unknown coefficients of the numerical solution, the number of collocation points employed in
h
IGA-L can be larger than that of the unknown coefficients. IGA-L has some advantages over
t
a IGA-C. First, if the number of coefficients of the numerical solution is fixed, a small increase
m
in the number of collocation points in IGA-L leads to a large improvement in the accuracy of
[ itsnumericalsolutioncomparedtoIGA-C.Second,becausethenumberofcollocationpointsin

1 IGA-Lisvariable,itismoreflexibleandmorestablethanIGA-C.Wealsoshowtheconsistency
v andconvergencepropertiesofIGA-Linthispaper.
4
4 Keywords: Isogeometricanalysis,collocationmethod,least-squaresfitting,NURBS,
2 consistencyandconvergence
7
0
.
1
0
1. Introduction
6
1
: WhileclassicalFiniteElementAnalysis(FEA)methods,widelyemployedinphysicalsimula-
v
tions, are based on element-wise piecewise polynomials, computer-aided design (CAD) mod-
i
X els are usually represented by non-uniform rational basis splines (NURBSs) with non-linear
r NURBS basis functions. Therefore, the first task in a CAD model simulation is to transform
a
the non-linear NURBS-represented CAD model into a linear mesh representation. This mesh
transformationisaverytediousoperation, andithasbecomethemosttime-consumingtaskof
thewholeFEAprocedure. Toavoidthemeshtransformationandadvancetheseamlessintegra-
tionofCADandcomputer-aidedengineering(CAE),isogeometricanalysis(IGA)wasinvented
∗Corresponding author: phone number: 86-571-87951860-8304, fax number: 86-571-88206681, email:
[email protected]
PreprintsubmittedtoElsevier January28,2016
byHughesetal.[1]. IGAisbasedonthenon-linearNURBSbasisfunctions,henceitcandeal
with NURBS-represented CAD models directly. In this way, IGA not only saves a significant
amountofcomputation,italsogreatlyimprovesthenumericalaccuracyofthesolution.
InIGA,theanalyticalsolutionT toaboundaryvalueproblemisapproximatedbyaNURBS
functionT withunknowncoefficients. (Forbrevity,weonlymentiontheboundaryvalueprob-
r
leminthispaper. However, thismethodisalsosuitablefortheinitialvalueproblem.) Solving
theboundaryvalueproblemisthenequivalenttodeterminingtheunknowncoefficientsofT . If
r
theorderofT ishigherthanthatofthedifferentialoperatorDoftheboundaryvalueproblem,
r
DT canberepresentedexplicitlybyaNURBSderivativeformula. Therefore,collocationmeth-
r
odscanbeappliedtothestrongformoftheboundaryvalueproblemtodeterminetheunknown
coefficientsofT .
r
In Ref. [2], an isogeometric collocation (IGA-C) method was presented. Suppose n is the
number of unknown coefficients of T . IGA-C first samples n values DT(η ), j = 1,2,··· ,n,
r j
and then generates a linear system of equations by making DT interpolate the n values, i.e.,
r
DT (η ) = DT(η ), j = 1,2,··· ,n. In this way, the unknown coefficients of T can be deter-
r j j r
minedbysolvingthelinearsystem.
Essentially, IGA-Cacquirestheunknowncoefficientsbyinterpolation, wherethenumberof
thecollocationpointsmustbeequaltothatoftheunknowncoefficients.Inthispaper,wepropose
theisogeometricleast-squarescollocation(IGA-L)method, whichallowsthenumberofcollo-
cationpointstobelargerthanthatoftheunknowncoefficients. Ityieldssomeadvantagesover
IGA-C.Insteadofinterpolation,IGA-Lmakesuseofapproximationtodeterminetheunknown
coefficientsofT . Specifically,IGA-LfirstsamplesmvaluesDT(η ), j = 1,2,··· ,m,wherem
r j
isgreaterthanthenumberofunknowncoefficients,i.e.,m>n. Thecoefficientsoftheunknown
solutionT arethendeterminedbysolvingtheleast-squaresfittingproblem
r
(cid:88)m (cid:13) (cid:13)
min (cid:13)(cid:13)DT(ηj)−DTr(ηj)(cid:13)(cid:13)2.
j=1
TherearetwoadvantagesofIGA-LoverIGA-C.First,numericaltestspresentedinthispaper
show that a small increase in the computation of IGA-L leads to a large improvement in the
numericalaccuracyofthesolution,eventhoughthecomputationalcostofIGA-Lisonlyslightly
morethanthatofIGA-C.Second,IGA-LismoreflexibleandmorestablethanIGA-C,because
IGA-Lallowsavariablenumberofcollocationpoints,whilethenumberofcollocationpointsin
IGA-Cisfixedtobeequaltothenumberofcontrolpoints.
Thestructureofthispaperisasfollows. InSection1.1,somerelatedworkisbrieflyreviewed.
In Section 2, the generic formulation of IGA-L is presented. In Section 3, we show the con-
sistencyandconvergencepropertiesofIGA-L.AfterthoroughlycomparingIGA-LandIGA-C,
bothintheoryandwithnumericalexamplesinSection4,thispaperisconcludedinSection5.
1.1. RelatedWork
Least-squares collocation methods: Although collocation-based meshless methods are ef-
ficient, equilibrium conditions are satisfied only at collocation points. Thus, collocation-based
meshlessmethodscanresultinsignificanterror. Toimprovecomputationalaccuracy, Zhanget
al. developedaleast-squarescollocationmethod[3],whereequilibriumconditionsholdatboth
collocation points and auxiliary points in a least-squares sense. In order to generate a better
conditioned linear algebraic equations, a least-squares meshfree collocation method was pro-
posed, based on first-order differential equations [4]. In Ref. [5], a point collocation method
2
was invented that employs the approximating derivatives based on the moving least-square re-
producing kernel approximations. Moreover, several schemes using least-squares collocation
methodsweredevelopedfortwo-andthree-dimensionalheatconductionproblems[6],transient
andsteady-statehyperbolicproblems[7],andadaptiveanalysisproblemsinlinearelasticity[8].
However,theconsistencyandconvergencepropertiesoftheaforementionedleast-squarescol-
locationmethodswereonlyvalidatedusingnumericalexamples,andtheoreticalnumericalanal-
yseswerenotreported. Inthispaper,wenotonlydevelopanisogeometricleast-squarescollo-
cationmethod,butalsoshowitsconsistencyandconvergencepropertiesintheory.
Isogeometric analysis: NURBS is a mathematical model for representing curves and sur-
facesbyblendingweightedcontrolpointsandNURBSbasisfunctions. Hence,theshapeofthe
NURBS curves and surfaces can be easily modified by moving their control points. Because
ofthedesirabletraitsofNURBSbasisfunctions,NURBSscurvesandsurfaceshavemanygood
properties,suchasconvexhull,affineinvariance,andvariationdiminishing.Moreover,NURBSs
can represent conic sections and quadric surfaces accurately. Therefore, NURBSs have been
widelyusedinCAD,computer-aidedmanufacturing(CAM),andCAE,andhavebecomeapart
of numerous industry standards, such as IGES, STEP, ACIS, and PHIGS. For more details on
NURBSs,thereareexcellentbooksonthesubject[9,10].
WhileaNURBSemploysnon-linearbasisfunctions,classicalFEAisbasedonelement-wise
piecewisepolynomials. Hence,whenanalyzingNURBS-representedCADmodelsusingclassi-
calFEA,theCADmodelshouldbediscretizedintoameshmodel. Themeshgenerationproce-
durenotonlyyieldsanapproximation,itisalsoverytedious,andhencehasbecomeabottleneck
inFEA.Toovercomethesedifficulties,Hughesetal. inventedtheIGAtechnique[1]. Because
IGA is based on NURBS basis functions, it can handle NURBS-represented CAD models di-
rectly without generating a mesh. Moreover, because NURBSs can represent complex shapes
(andphysicalfields)withsignificantlyfewercontrolpointsthanameshrepresentation,thecom-
putationalefficiencyandnumericalaccuracyofIGAarehigherthantheclassicalFEAmethod.
In addition, because of the knot insertion property of NURBSs, the original shape of the CAD
modelcanbemaintainedexactlyintherefinementprocedure[1].
Ingeometricdesign,theNURBSrepresentationisusuallyemployedtomodelcurvesandsur-
faces. There are a few effective methods in geometric design for modeling spline solids. To
modelaNURBSsolidforIGAapplications,Zhangetal. proposedasolidconstructionmethod
fromaboundaryrepresentation[11]. InthestudyofIGA,makingthegeometricrepresentation
moresuitableforanalysispurposesisakeyresearchgoal. Cohenetal. presentedtheanalysis-
awaremodelingtechnique[12]. Moreover,T-spline[13,14],trimmedsurfaces[15],andsubdi-
visionsolids[16]werealsoemployedintheIGAmethodtomodelthephysicaldomain.
Currently, the IGA method has been successfully applied in various simulation problems,
including elasticity [17, 18], structure [19, 20, 21], and fluids [22, 23, 24]. On the other hand,
someworkconcernsthecomputationalaspectoftheIGAmethod, forinstance, toimprovethe
accuracyandefficiencybyreparameterizationandrefinement[25,26,27,28,29].
Isogeometriccollocationmethods:Thecollocationmethodisasimplenumericalmethodfor
solvingdifferentialequationsthatcangenerateasolutionthatsatisfiesthedifferentialequation
at a set of discrete points, called collocation points [30]. If the order of the unknown NURBS
functionthatisemployedtoapproximatethesolutionofthedifferentialequationishighenough,
the collocation method can be applied to the strong form of the differential equation. In this
way,theIGA-CmethodwaspresentedinRef.[2]. However,IGA-Cisgreatlyinfluencedbythe
collocation points. Recently, a comprehensive study on IGA-C discovered its superior behav-
iorovertheGalerkinmethodintermsofitsaccuracy-to-computational-timeratio[31], andthe
3
consistencyandconvergencepropertiesoftheIGA-CmethodweredevelopedinRef.[32].
Moreover,basedonthelocalhierarchicalrefinementofNURBSs,adaptiveIGA-Cswerede-
velopedandanalyzed[31]. Meanwhile, IGA-Chasalsobeenextendedtomulti-patchNURBS
configurations, various boundary and patch interface conditions, and explicit dynamic analy-
sis [33]. Recently, IGA-C was successfully employed to solve the Timoshenko beam prob-
lem[34]andspatialTimoshenkorodproblem[35],showingthatmixedcollocationschemesare
locking-free,independentlyofthechoiceofthepolynomialdegreesfortheunknownfields.
2. GenericFormulationofIGA-L
Considerthefollowingboundaryvalueproblem,

DT = f, inΩ∈Rd,
(1)
GT =g, on∂Ω,
whereΩisthephysicaldomaininRd,Disadifferentialoperatorinthephysicaldomain,GT =g
is the boundary condition, and f : Ω → R and g : ∂Ω → R are given functions. Suppose d
1
isthemaximumorderofderivativesappearinginD : V → W,whereV andW aretwoHilbert
spaces,andtheanalyticalsolutionT ∈Cd2(Ω),d2 >d1 ≥1.
InIGA,thephysicaldomainΩisrepresentedbyaNURBSmapping:
F:Ω →Ω, (2)
0
where Ω is the parametric domain. By replacing the control points of F(Ω ) with unknown
0 0
controlcoefficients,therepresentationoftheunknownnumericalsolutionT isgenerated.
r
SupposetherearenunknowncontrolcoefficientsintheunknownnumericalsolutionT . We
r
first sample m points θ inside the parametric domain Ω that correspond to m values η =
1 i 0 1 i
F(θ), i = 1,2,··· ,m ,insidethephysicaldomainΩ. Furthermore,wesamplem pointsθ on
i 1 2 j
theparametricdomainboundary∂Ω thatcorrespondtom valuesη = F(θ ), j=m +1,m +
0 2 j j 1 1
2,··· ,m +m , on the physical domain boundary ∂Ω. The total number of these points, i.e.,
1 2
m = m +m ,isgreaterthanthenumberofunknowncoefficientsofthenumericalsolutionT ,
1 2 r
namely,m = m +m > n. JustasinIGA-C,thesesamplingpointsarealsocalledcollocation
1 2
points.
Substitutingthesecollocationpointsintotheboundaryvalueproblem(1), asystemofequa-
tionswithmequationsandnunknownsisobtained(wherem=m +m >n),
1 2

DTr(ηi)= f(ηi), i=1,2,···m1,
(3)
GT (η )=g(η ), j=m +1,m +2,··· ,m +m .
r j j 1 1 1 2
Arranging the unknowns of the numerical solution T into an n × 1 matrix, i.e., X =
r
[x x ··· x ]t,thesystemofequations(3)canberepresentedinmatrixformby
1 2 n
AX =b.
Becausethenumberofequationsisgreaterthanthenumberofunknowns,thesolutionissought
intheleast-squaressense,i.e.,
min(cid:107)AX−b(cid:107)2. (4)
X
4
Computation of the least-squares problem (4): The least-squares problem (4) is very im-
portantinpractice,andtherehavebeendevelopedlotsofrobustandefficientmethodsforsolving
it[36]. Onefrequentlyemployedmethodistosolvethenormalequation(5),
AtAX = Atb. (5)
Althoughtheconditionnumberofthematrix ATAisthesquareofthatofthematrix A, ATAis
a symmetric positive definite matrix, and the normal equation (5) can be solved efficiently by
Choleskydecomposition[36]. Moreover,HouseholderorthogonalizationandGivenorthogonal-
ization are also two efficient methods [36] which often employed in solving the least-squares
problem(4). Formoremethodsonsolving(4),pleasereferto[36].
In the following, we will show the consistency and convergence properties of IGA-L, and
compareitwithIGA-C,bothintheoryandwithnumericalexamples.
3. NumericalAnalysis
IntheIGA-Lmethoddevelopedabove,aNURBSfunctionT isemployedtoapproximatethe
r
analyticalsolutionT oftheboundaryvalueproblem(1),andhencetherealdifferentialoperator
DT isapproximatedbyDT . InRef.[32],itwasshownthatbothDT andT aredefinedonthe
r r r
sameknotintervals,i.e.,
Lemma1. IfT isaNURBSfunctionandDisadifferentialoperator,thenbothDT andT are
r r r
definedonthesameknotintervals[32].
Moreover,givenasetΦ,thediameterofΦisdefinedas
diam(Φ)=max{d(x,y),x,y∈Φ},
whered(x,y)istheEuclideandistancebetweenxandy.Inthefollowing,wesupposetheNURBS
functionT tobedefinedonaknotgridTh,wherehistheknotgridsize:Intheone-dimensional
r
case,Th isaknotsequence,andh = maxi{diam([ui,ui+1])};inthetwo-dimensionalcase,Th is
arectangulargrid,andh=maxij{diam([ui,ui+1]×[vj,vj+1])};andinthethree-dimensionalcase,
Thisahexahedralgrid,andh=maxijk{diam([ui,ui+1]×[vj,vj+1]×[wk,wk+1])}.
Inthissection,westudytheconsistencyandconvergencepropertiesofIGA-L,i.e.,whenh→
0, notonlywillapproximatedifferentialoperatorDT tendtorealoperatorDT, butnumerical
r
solutionT willalsotendtoanalyticalsolutionT.
r
3.1. Consistency
In this section, we explore the consistency property of the IGA-L method. The theorem for
consistencyispresentedasfollows.
Theorem1. IntheIGA-Lmethod,if
(1.1) eachknotintervaloftheNURBSfunctionT definedonknotgridTh containsatleastone
r
collocationpoint
(1.2) thedegreeofeachvariableinT islargerthanthemaximumorderofthepartialderivatives
r
tothevariablesappearinginD(seeEq.(1))
5
(1.3) theleast-squaresfittingerrore isbounded,i.e.,e < M¯,whereM¯ isapositiveconstant
h h
thenapproximatedifferentialoperatorDT tendstorealdifferentialoperatorDT intheL2norm
r
whenh→0.
Proof: We only show the theorem in the two-dimensional case. The proof for the one- and
three-dimensionalcasesissimilar.
Inthetwo-dimensionalcase, supposethetensorproductNURBSfunctionT (u,v)ofdegree
r
l ×l isdefinedontheknotsequences
u v
{u ,u ,··· ,u ,u ,··· ,u ,u ,u ,··· ,u },
(cid:124)0(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)0(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)0 1 nu−1 (cid:124)n(cid:32)(cid:32)u(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)n(cid:32)(cid:32)u(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)n(cid:125)u
lu+1 lu+1 (6)
{v ,v ,··· ,v ,v ,··· ,v ,u ,u ,··· ,u }.
(cid:124)0(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)0(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)0 1 nv−1 (cid:124)n(cid:32)(cid:32)v(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)n(cid:32)(cid:32)v(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)n(cid:125)v
lv+1 lv+1
Thenthecorrespondingknotgridis
Th ={[ui,ui+1]×[vj,vj+1], i=0,1,··· ,nu−1, j=0,1,··· ,nv−1}. (7)
Denote
R(u,v)=(DT(u,v)−DT (u,v))2,(u,v)∈[u ,u ]×[v ,v ].
r 0 nu 0 nv
Note that T (u,v) is generated by least-squares fitting the values of DT(u,v) at the collocation
r
pointsϑ =(η ,ξ ),i.e.,DT(ϑ ),d =1,2,··· ,D,(D≥n n ). Andsupposethefittingerroris
d d d d u v
(cid:88)D (cid:88)D
e = R(ϑ )= (DT(ϑ )−DT (ϑ ))2. (8)
h d d r d
k=1 k=1
First, based on Lemma 1, the numerical solution T (u,v) and the approximate differential
r
operatorDT (u,v)havethesameknotintervals,
r
[ui,ui+1]×[vj,vj+1], i=0,1,··· ,nu−1, j=0,1,··· ,nv−1.
Now,considertheerrorbetweenDT(u,v)andDT (u,v)intheL2norm,
r
(cid:90) vnv (cid:90) unu n(cid:88)v−1n(cid:88)u−1(cid:90) vj+1(cid:90) ui+1
(cid:107)DT(u,v)−DT (u,v)(cid:107)2 = R(u,v)dudv= R(u,v)dudv.
r L2
v0 u0 j=0 i=0 vj ui
Sinceeachknotinterval[ui,ui+1]×[vj,vj+1]containsatleastonecollocationpoint,wesuppose
ϑd =(ηd,ξd)∈[ui,ui+1]×[vj,vj+1].Usingtheleftandrightrectangleintegralformularepeatedly,
6
weget
(cid:90) vj+1(cid:90) ui+1 (cid:90) vj+1 (cid:32)(cid:90) ηd (cid:90) ui+1 (cid:33)
R(u,v)dudv= dv R(u,v)du+ R(u,v)du
vj ui vj ui ηd
=(cid:90) vj+1(ηd−ui)R(ηd,v)+(ui+1−ηd)R(ηd,v)+(ηd−ui)2R(cid:48)u(µ(i12)(v),v) +(ui+1−ηd)2R(cid:48)u(µ(i22)(v),v)dv
vj
=(ui+1−ui)(cid:90) vj+1R(ηd,v)dv+(ηd−ui)2(cid:90) vj+1 R(cid:48)u(µ(i1)(v),v)dv+(ui+1−ηd)2(cid:90) vj+1 R(cid:48)u(µ(i2)(v),v)dv
2 2
vj vj vj
 
=(ui+1−ui)(ξd−vj)R(ηd,ξd)+(vj+1−ξd)R(ηd,ξd)+(ξd−vj)2R(cid:48)v(ηd2,ω(i1j)) +(vj+1−ξd)2R(cid:48)v(ηd2,ω(i2j))
+(ηd−ui)2(cid:90) vj+1 R(cid:48)u(µ(i1)(v),v)dv+(ui+1−ηd)2(cid:90) vj+1 R(cid:48)u(µ(i2)(v),v)dv
2 2
vj vj
 
=(ui+1−ui)(vj+1−vj)R(ηd,ξd)+(ui+1−ui)(ξd−vj)2R(cid:48)v(ηd2,ω(i1j)) +(vj+1−ξd)2R(cid:48)v(ηd2,ω(i2j))
+(ηd−ui)2(cid:90) vj+1 R(cid:48)u(µ(i1)(v),v)dv+(ui+1−ηd)2(cid:90) vj+1 R(cid:48)u(µ(i2)(v),v)dv,
2 2
vj vj
whereµ(i1)(v),µ(i2)(v)∈(ui,ui+1)andω(i1j),ω(i2j) ∈(vj,vj+1).Bythemeanvaluetheoremofintegral,
thereexist(µ¯(i1j),ω¯(i1j)),(µ¯(i2j),ω¯(i2j))∈[ui,ui+1]×[vj,vj+1]suchthat
(cid:90) vj+1 R(cid:48)u(µ(1)(v),v)dv=(vj+1−vj)R(cid:48)u(µ¯(i1j),ω¯(i1j)), and
2 2
vj
(cid:90) vj+1 R(cid:48)u(µ(2)(v),v)dv=(vj+1−vj)R(cid:48)u(µ¯(i2j),ω¯(i2j)),
2 2
vj
Therefore,
(cid:90) vj+1(cid:90) ui+1
R(u,v)dudv=(ui+1−ui)(vj+1−vj)R(ηd,ξd)
vj ui
 
+(ui+1−ui)(ξd−vj)2R(cid:48)v(ηd2,ω(i1j)) +(vj+1−ξd)2R(cid:48)v(ηd2,ω(i2j))
 
+(vj+1−vj)(ηd−ui)2R(cid:48)u(µ¯(i1j2),ω¯(i1j)) +(ui+1−ηd)2R(cid:48)u(µ¯(i2j2),ω¯(i2j)).
Moreover,wedenoteΞ=[u ,u ]×[v ,v ]. Itiseasytoshowthat
0 nu 0 nv
(cid:12) (cid:12) (cid:12) (cid:12)
min (cid:12)(cid:12)(cid:12)R(cid:48)(u,v)(cid:12)(cid:12)(cid:12)≤n(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i2j))(cid:12)(cid:12)(cid:12) ≤ max (cid:12)(cid:12)(cid:12)R(cid:48)(u,v)(cid:12)(cid:12)(cid:12),
(u,v)∈Ξ v j=0 i=0 (unu −u0)(vnv −v0) 2 (u,v)∈Ξ v
7
(cid:12) (cid:12) (cid:12) (cid:12)
min (cid:12)(cid:12)(cid:12)R(cid:48)(u,v)(cid:12)(cid:12)(cid:12)≤n(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)u(µ¯(i1j),ω¯(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)u(µ¯(i2j),ω¯(i2j))(cid:12)(cid:12)(cid:12) ≤ max (cid:12)(cid:12)(cid:12)R(cid:48)(u,v)(cid:12)(cid:12)(cid:12).
(u,v)∈Ξ u j=0 i=0 (unu −u0)(vnv −v0) 2 (u,v)∈Ξ u
Andthen,basedontheintermediatevaluetheorem,thereexist(η(1),ξ(1))∈Ξand(η(2),ξ(2))∈Ξ,
suchthat
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12)(cid:12)(cid:12)R(cid:48)(η(1),ξ(1))(cid:12)(cid:12)(cid:12)=n(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i2j))(cid:12)(cid:12)(cid:12),
v (u −u )(v −v ) 2
j=0 i=0 nu 0 nv 0
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12)(cid:12)(cid:12)R(cid:48)(η(2),ξ(2))(cid:12)(cid:12)(cid:12)=n(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)u(µ¯(i1j),ω¯(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)u(µ¯(i2j),ω¯(i2j))(cid:12)(cid:12)(cid:12).
u (u −u )(v −v ) 2
j=0 i=0 nu 0 nv 0
Asaresult,
(cid:90) vnv (cid:90) unu n(cid:88)v−1n(cid:88)u−1(cid:90) vj+1(cid:90) ui+1
(cid:107)DT(u,v)−DT (u,v)(cid:107)2 = R(u,v)dudv= R(u,v)dudv
r L2
v0 u0 j=0 i=0 vj ui
n(cid:88)v−1n(cid:88)u−1
= (ui+1−ui)(vj+1−vj)R(ηd,ξd)
j=0 i=0
 
+n(cid:88)v−1n(cid:88)u−1(ui+1−ui)(ξd−vj)2R(cid:48)v(ηd2,ω(i1j)) +(vj+1−ξd)2R(cid:48)v(ηd2,ω(i2j))
j=0 i=0
 
+n(cid:88)v−1n(cid:88)u−1(vj+1−vj)(ηd−ui)2R(cid:48)u(µ¯(i1j2),ω¯(i1j)) +(ui+1−ηd)2R(cid:48)u(µ¯(i2j2),ω¯(i2j))
j=0 i=0
(cid:12) (cid:12) (cid:12) (cid:12)
≤h2(cid:88)D R(ηd,ξd)+hn(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)v(ηd,ω(i2j))(cid:12)(cid:12)(cid:12)
2
d=1 j=0 i=0
(cid:12) (cid:12) (cid:12) (cid:12)
+hn(cid:88)v−1n(cid:88)u−1(ui+1−ui)(vj+1−vj)(cid:12)(cid:12)(cid:12)R(cid:48)v(µ¯(i1j),ω¯(i1j))(cid:12)(cid:12)(cid:12)+(cid:12)(cid:12)(cid:12)R(cid:48)v(µ¯(i2j),ω¯(i2j))(cid:12)(cid:12)(cid:12)
2
j=0 i=0
(cid:16)(cid:12) (cid:12) (cid:12) (cid:12)(cid:17)
=h2eh+h(unu −u0)(vnv −v0) (cid:12)(cid:12)R(cid:48)v(η(1),ξ(1))(cid:12)(cid:12)+(cid:12)(cid:12)R(cid:48)u(η(2),ξ(2))(cid:12)(cid:12)
(cid:16)(cid:12) (cid:12) (cid:12) (cid:12)(cid:17)
<h2M¯ +h(unu −u0)(vnv −v0) (cid:12)(cid:12)R(cid:48)v(η(1),ξ(1))(cid:12)(cid:12)+(cid:12)(cid:12)R(cid:48)u(η(2),ξ(2))(cid:12)(cid:12) .
Ontheotherhand,sincethedegreesofuandvinT (u,v)arebothlargerthanthemaximum
r
ordersofthepartialderivativestouandvappearinginD(referto(1)),respectively,similaras
theone-dimensionalcase,R(cid:48)(u,v)andR(cid:48)(u,v)arebothcontinuous,andthenboundedonΩ∪∂Ω,
u v
i.e.,
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12)(cid:12)R(cid:48)(η(1),ξ(1))(cid:12)(cid:12)≤ Mˆ, and, (cid:12)(cid:12)R(cid:48)(η(2),ξ(2))(cid:12)(cid:12)≤ Mˆ,
v u
whereMˆ isapositiveconstant.
Inconclusion,
(cid:107)DT(u,v)−DT (u,v)(cid:107)2 ≤h2M¯ +2h(u −u )(v −v )Mˆ,
r L2 nu 0 nv 0
8
andTheorem1isproved. (cid:50)
3.2. Convergence
Based on the consistency property of IGA-L, we can show that IGA-L is convergent if the
differentialoperatorisstableorstronglymonotonic,similarlyasinRef.[32].
LetV andW betwoHilbertspacesand(cid:107)·(cid:107) and(cid:107)·(cid:107) betwonormsdefinedonV andW,re-
V W
spectively.Suppose(cid:107)·(cid:107) and(cid:107)·(cid:107) areequivalenttotheL2norm,i.e.,thereexistpositiveconstants
V W
α , β , α ,andβ suchthat,
V V W W
α (cid:107)·(cid:107) ≤(cid:107)·(cid:107) ≤β (cid:107)·(cid:107) , (9)
V L2 V V L2
α (cid:107)·(cid:107) ≤(cid:107)·(cid:107) ≤β (cid:107)·(cid:107) . (10)
W L2 W W L2
Wefirstgivethedefinitionsforstableandstronglymonotonicoperators.
Definition1(Stabilityestimateandstableoperator[37]). LetV,WbeHilbertspacesandD:
V →W beadifferentialoperator. IfthereexistsaconstantC >0suchthat
S
(cid:107)Dv(cid:107) ≥C (cid:107)v(cid:107) , forallv∈ D(D), (11)
W S V
where D(D) represents the domain of D, differential operator D is called the stable operator,
andtheinequality(11)iscalledthestabilityestimate.
Definition2(Stronglymonotonicoperator[37]). LetV beaHilbertspaceandD∈L(V,V(cid:48)).
OperatorDissaidtobeastronglymonotonicoperator,ifthereexistsaconstantC > 0,such
D
that
(cid:104)Dv,v(cid:105)≥C (cid:107)v(cid:107)2 , forallv∈V. (12)
D V
For every v ∈ V, element Dv ∈ V(cid:48) is of a linear form. The symbol (cid:104)Dv,v(cid:105), which denotes the
applicationofDvtov∈V,iscalledadualitypairing.
Clearly,ifadifferentialoperatorDisstronglymonotonic,itisstable.
Lemma2. Let V be a Hilbert space and D ∈ L(V,V(cid:48)) be a continuous strongly monotonic
linear operator. Then, there exists a constant C > 0 such that D satisfies the stability esti-
D
mate(11)[37].
TheproofcanbefoundinRef.[32].
Therefore,wehavetheconvergencepropertyofIGA-Lasfollows.
Theorem2. SupposeNURBSfunctionT ,definedonknotgridTh,isthenumericalsolutionto
r
theboundaryvalueproblem(1),generatedbyIGA-L,andtheconditionspresentedinTheorem1
aresatisfiedinone,two,andthreedimensions,respectively. IfdifferentialoperatorD:V →W
in(1)isastableoperator,T willconvergetoanalyticsolutionT whentheknotgridsizeh→0.
r
Proof: Differential operator D in (1) is a stable operator, so there exists a constantC > 0,
S
suchthat
(cid:107)D(T −T )(cid:107) ≥C (cid:107)T −T (cid:107) .
r W S r V
9
Anditisequivalentto
1
(cid:107)T −T (cid:107) ≤ (cid:107)DT −DT (cid:107) .
r V C r W
S
Duetotheequivalenceof(cid:107)·(cid:107) and(cid:107)·(cid:107) (9),wehave,
L2 W
1 β
(cid:107)T −T (cid:107) ≤ (cid:107)DT −DT (cid:107) ≤ W (cid:107)DT −DT (cid:107) .
r V C r W C r L2
S S
BecauseoftheconsistencyofIGA-L(Theorem1), (cid:107)T −T (cid:107) willconvergeto0whenh → 0.
r V
Andthistheoremisproved.
(cid:50)
Moreover,Theorem2andLemma2leadtothedirectcorollary.
Corollary1. SupposeNURBSfunctionT definedonknotgridTh isthenumericalsolutionto
r
theboundaryvalueproblem(1),generatedbyIGA-L,andtheconditionspresentedinTheorem1
are satisfied in one, two, and three dimensions, respectively. Additionally, suppose norm (cid:107)·(cid:107)
L2
boundsnorm(cid:107)·(cid:107) .IfdifferentialoperatorDin(1)isastronglymonotonicoperator,thenT will
V(cid:48) r
convergetoanalyticsolutionT whenh→0.
Itiswellknownthatawideclassofellipticdifferentialoperatorsarestableorstronglymono-
tonic. Hence, IGA-L is convergent for equations that have these elliptic differential operators.
The examples of a PDE whose differential operators are strongly monotonic can be found in
Refs.[37,32].
4. ComparisonswithIGA-C
4.1. Theoreticalcomparisons
Inthefollowing,wecompareIGA-LwithIGA-Cintermsoftheircomputationalefficiencyat
solvingascalarproblem(Laplaceequation[31])andvectorproblem(elasticityequation[31]).
Weconsidermodeldiscretizationsinone,two,andthreedimensionsthatarecharacterizedbythe
degreeofthebasisfunctionsandthenumbersofcontrolandcollocationpointsineachparametric
direction. Forthesakeofsimplicity,weassumethatthemodeldiscretizationsintwoandthree
dimensionshavethesamenumberofcollocationpoints(andcontrolpoints)ineachparametric
direction,andthedegreesofbasisfunctionsinone,two,andthreedimensionsare p, p×p,and
p×p×p,respectively.
First, the costs in floating point operations (flops) for the formation at one collocation point
arethesameforbothIGA-LandIGA-C.ThesecostsarelistedinTable1. Itshouldbepointed
out that there are some slight flaws in the formulae employed to calculate the costs in flops in
Ref.[31]. Hence,theformulaepresentedinTable1areabitdifferentfromthosein[31].
Second,wecomparedthecosts(inflops)ofsolvingthelinearsystemofequationsinIGA-C
andIGA-L,andpresentthecomparisoninTable2.Inthistable,thefirstcolumnisthedimension
oftheproblemsolved,thesecondcolumnisthenumberofcontrolpoints(equaltothenumberof
collocationpointsinIGA-C),andthethirdcolumnisthecostinflopstosolvethelinearsystem
of equations using Gaussian elimination in IGA-C [36]. In addition, the fourth column is the
number of control points, and the fifth column is the number of collocation points in IGA-L.
Finally,thelastcolumnisthecostinflopstosolvenormalequation(5)usingCholeskydecom-
positioninIGA-L[36].WecanseethatthecosttosolveEq.(5)inIGA-Llinearlyincreaseswith
thenumberofcollocationpoints(m, m×m, m×m×m).
10