Elliptic Grid Generation with B-Spline Collocation
PhilippLamby
Karl-HeinzBrakhage
InstituteforGeometryand AppliedMathematis
RWTH AahenUniversity
Templergraben55,52078Aahen, Germany
flamby,brakhagegigpm.rwth-a ahe n.de
Abstract
Werevisitthelassialtehniqueofelliptigridgenerationwithharmoni
mappings. Forthedeterminationoftheontrolfuntionsweusetheframe-
work developedbySpekreijse [1℄. However,insteadwithnite dierenes
we disretize the underlying partial dierential equation with a B-spline
olloationmethodinordertoworkdiretlywiththenativedatarepresen-
tations of ourCAGD system. This way weanmake useof the sparsity
and aurayof theB-spline boundaryrepresentationsand guaranteethe
geometrionsistenyofourCADmodels. Inthispaperwewillsummarize
theunderlyingalgorithmsandpresentsomerstappliationexamples.
Introduction
Inthe ontextof thedevelopment ofanewadaptiveNavier-Stokessolver
quadflow whih aims at the simulation of uid-struture interation at
airplanewings,f. [3℄,anewgridgenerationmodulehasbeenimplemented
whihisbasedontherepresentationofthegeometrywithparametrimap-
pings. AsinmanyommonCADsystemsfreeformurvesandsurfaesare
represented by B-splines. From this parametri representations adaptive
gridsareomputedbyfuntion evaluation. Conretely,withinthissystem
urvesarerepresentedbyB-splinesoftheform
x(u)= N
X
i=0 p
i N
i;p;U
(u); (1)
planargridsandsurfaesaremodeledbybivariateB-splinetensorproduts
x(u;v)= N
X M
X
p
i;j N
i;p;U (u)N
j;q;V
(v) (2)
x(u;v;w)= N;M;L
X
i;j;k =0 p
i;j;k N
i;p;U (u)N
j;q;V (v)N
k ;r;W
(w): (3)
HereU;V;W arenon-dereasingandnon-stationaryknotsequenes,i.e.,
U =(u
i )
N+p 2
i=0
: u
0 u
1
:::u
N+p 2
; u
i+p
<u
i
; (4)
andN
i;p;U
denotesthei-thnormalizedB-spline-funtionoforderpdened
bythereursion
N
i;1;U (u)=
(
1 ifu
i
u<u
i+1
0 otherwise
; (5)
N
i;p;U (u)=
u u
i
u
i+p 1 u
i N
i;p 1 (u)+
u
i+p u
u
i+p u
i+1 N
i+1;p 1
(u): (6)
B-splineurvesarepieewisepolynomialsofdegreep 1. Usuallywehoose
p = 4, i.e. ubi splines, and knot sequenes with p-fold knots at the
interval ends. This has the advantage, that the knot sequene interval
oinides with theparameterintervalof the urve, that therst and last
ontrol pointoinide withthestartandend pointof theurve,and that
therstandlastspanoftheontrolpolygonaretangentialtotheurveat
thestartandendpointoftheurve. Inmanypratialasesnon-uniform
knotsequenes are onstruted, for instane in the ourse of anadaptive
approximationor interpolationproedure. Multiple interiorknots an be
usedtomodelnon-smoothfeaturesofanobjet.
Grid Generation Equations
Inordertogeneratesmoothgrids,alltheommontehniquesofstrutured
gridgeneration,namelytransniteinterpolationandmethodsbasedonthe
solutionofpartialdierentialequationsareapplied. Inthispaperwewant
to integrate anellipti grid generation tehniqueinto our CAGD system.
OurhoiewasforSpekreijse'sapproahwhihanbeverybrieysumma-
rizedas follows: letx(s) be aharmoni mappingfrom the d-dimensional
parameterspaeP ontothephysialdomainCands()beaso-alledon-
trolmappingfromtheomputationaldomainContotheparameterdomain
P. Thentheompositemapping
x(s()): C !D (7)
fulllsadierentialequationoftheform
L(x)= d
X
g
ij
2
x
i
j +
d
X
P
k x
k
=0
(8)
P
k
= d
X
i;j=1 J
2
g ij
P k
ij
; (9)
J =detx 0
() is the Jaobianof the omposite mapping, the g
ij and g
ij
aretheovariantandontravariantmetritensorsdenedby
g
ij
= x
i
x
j
; d
X
k =1 g
ik g
k j
=Æ
ij
; (10)
theP k
ij
aretheomponentsofthevetor
P
ij
= T
1
2
s
i
j
; (11)
and T =s 0
() istheJaobianmatrixoftheontrol mapping. The aimof
thispaperistosolvethisPDEwithaB-splineolloationmethod.
B-Spline Collocation
Thegeneralideaofolloationistodetermineafuntion,sothatitexatly
satisesthedierentialequationatertainpoints,theolloationpoints. In
awayolloationissimilartointerpolation,butinontrasttointerpolation
we donotmath funtion valuesbut ertain ombinationof funtion and
derivative values. In order to simplify the notation we onentrate onto
the bivariate ase from now on and denote the Cartesian oordinates of
the omputational domain, the unit square, with = (u;v) and of the
parameterdomainwiths=(s;t). Hene,wesearhafuntionoftheform
(2)whihfullls
Lx(^u
i
;v^
j
)=0; i=1;:::;N 1; j=1;:::;M 1 (12)
atertainolloationpointsu^
i ,^v
j
andDirihletboundaryonditionsforthe
ontrol pointsp
0;j , p
N;j
,j =0;:::;M andp
i;0 ,p
i;M
, i=0;:::;N. The
task is nowto hoose appropriateolloation pointsfor theonguration
athand.
The most popular B-spline olloation sheme is Gauÿ olloation with
ubiHermite-splines. Itsappliationtoelliptigridgenerationhasalready
beeninvestigatedbyManke[2℄. Herein eahknotintervaltheolloation
points arethe absissaeof Gaussian quadraturerules. An obviousaveat
isthat in apreproessing steponehastomakeallinteriorknots two-fold
by knot-insertion, and that therefore the resulting grid will be only C 1
.
u
i
= i+p
X
k =i+1 u
k
; (13)
asolloationpoints. ThishoieismotivatedbytheShoenberg-Whitney
theorem,seereferene[4℄,whihsaysthattheinterpolationproblemx(^u
i )=
f
i
is well posed if, and only if, everyu^
i
lies in the support of the i th
B-splinefuntion,i.e., ifN
i (^u
i
)>0. Asoneaneasilyverify,theGreville
absissae always give a set of as many distint points, as the spline has
ontrol points and they fulll the onditions of the Shoenberg-Whitney
theorem. Thisolloationshemefreesusfrom theneessitytoinsertad-
ditionalknotsinourtensorprodutrepresentations. Adisadvantage,how-
ever,isthatolloationattheGrevilleabsissaedoesnothavetheoptimal
onsistenyorder.
The Shoenberg-Whitney theorem is also the reason why we do not use
astandardnite diereneodefollowedbyaninterpolationalgorithm in
order to generateellipti B-spline grids: a typialnite dierene ode is
basedon the assumption that thegrid points x
i;j
are numerialapproxi-
mationsof regular spaed valuesx(ih
u
;jh
v
). However,depending onthe
strutureof theunderlying splineit ouldbeomeneessaryto work with
unevenlyspaedgridsinordertofulllthestipulationsoftheShoenberg-
Whitneytheorem duringtheinterpolationproess.
Application Example
Theafore-mentionedolloationshemeshavebeenimplementedandtested
forplanargrids,surfaesandvolumegrids. Inordertosolvethesystemswe
just followthestandardapproah anduseaxedpoint iteration,freezing
themetrioeientsinEquation(8)inordertogetalinearsysteminevery
single iteration. Then we apply the olloation sheme to the linearized
equations. Thearisingsparselinearsystemsaresolvedwithadiretsolver.
This kind of implementation is not well suited to solvebig systems with
maximumeieny, buttheaimofthe urrentstudywasnottoompare
theeieny of theimplementation (thishasalreadybedonein [2℄),but
tostudytheprinipalvalidityofthemethod.
Asarstappliationwepresentthegridinablokthatistakenfromagrid
foradual-bellonguration,seeFigure1. Theboundariesareapproxima-
tivelyparameterizedbyarlength,sothat weanusetheidentimapping
asontrolmapping. Hene,allontrolfuntionsP k
ij
arezeroandtheresult-
inggridmappingisharmoni. However,thespaingof theontrolpoints,
resolvethedierent featuresof thenozzle ontourandfrom the neessity
tomutuallyinserttheknotswhiharenotpresentintherepresentationof
the opposite boundary in order to build a tensor produt. However, the
resultingnumerialgrid, whih is omputedby evaluation of theB-spline
funtion hasthedesiredsmoothnessproperties.
Figure1. Control pointsand harmoni meshfor thedual bell.
Boundary Orthogonality
InSpekreijse'sapproahthere remains theproblem to determinesuitable
ontrol funtions in order to inorporate desired features into the grids.
In order to nd a ontrol mapping that ensures boundary orthogonality
Spekreijseproposestoproeedasfollows. Letusassumethatafolding-free
gridx()isalreadyavailable. Thisgridmaybe,forinstane,atransnite
interpolantorthesolutionofthepurelyharmonigridgenerationsystem.
Thenin therststepwesolvethetransformedLaplae equation
div(Agrads)=0 (14)
where
A=J
g 11
g 12
g 12
g 22
= 1
J
g
22 g
12
g g
: (15)
ditions, in partiularwerequires=n=0at theboundariesx(u;0) and
x(u;1) and t=n=0at the boundariesx(0;v) andx(1;v) ofthe physi-
al domain. Thesolution ofthis problem givesus aone-to-one boundary
mappingC !P. Intheseondstepweomplete thisboundarymap-
pingtoasuitableontrolmappingthatfulllstheorthogonalityonditions
t=u= 0along the boundaries u= 0and u= 1 and s=v = 0along
theboundariesv=0andv=1usinganalgebraigridgenerationmethod.
Forthedetailsofthismethod wehavetoreferto[1℄.
Whereasthedisretization ofEquation8byolloationisstraightforward
itismoreonvenienttodisretizeequation14withanitevolumemethod.
For this we observe that for any ontrol volume in the omputation
domaintheequation
Z
(grads;An)d=0 (16)
holds. Of ourse, asontrol volumes we will hoose intervalsof the form
[u
i
;u
i+1
℄[v
j
;v
j+1
℄. Againwe wanttorepresent theontrol mappingas
tensor produtB-spline. Therefore, theintegralovertheboundaryofthe
ontrolvolumeisomposedofsegmentsoftheform
Z
ui+1
u
i
(grads;An)d= Z
ui+1
u
i
s
u
s
v
;A
0
1
d
= X
i;j p
ij
N
j (v)
Z
ui+1
ui N
0
i (u)
g
12
J
du+N 0
j (v)
Z
u2
u1 N
i (u)
g
11
J du
and
Z
vj+1
vj
(grads;An)d= Z
vj+1
vj
s
u
s
v
;A
1
0
d
= X
i;j p
ij
N
i (u)
Z
vj+1
v
j N
0
j (v)
g
12
J
dv+N 0
i (u)
Z
vj+1
v
j N
j (v)
g
22
J du
:
The integralsin the brakets enter the matrixof the disretized problem
and anheaplybe evaluated by quadratureformulas. Inorder to get as
many equations as ontrol points weenter the ontrol volumina around
theGrevilleabsissaebyhoosing
u
i
= 1
2 (u
i 1 +u
i
); i=1;:::;N 1; u
0
=0; u
N
=1
v
i
= 1
2 (v
j 1 +v
j
); i=1;:::;N 1; v
0
=0; v
M
=1
(17)
Theboundaryonditions=ntransformsto(grads;An)=0attheorre-
evaluationoftheresultingorthogonalgrid,Figure3showstheorrespond-
ingontrolmappingandadetailviewatthenozzlethroat.
Figure2. Control points andorthogonal meshfor thedualbell.
Figure3. Control Mapand DetailView
Convergence and Stability Matters
Theaboveexampleshowsthatinprinipletheolloationmethodpresents
aviablemethodto integrateelliptigrid generationstrategyinto aspline
basedCAGDsystem. However,itturnsoutthattherearealsosomeom-
pliations. Firstsignsfor these problemsalreadyreveal themselvesin the
followingonvergenestudy.
Figure 4. Cosine Testase
Weonsider asimpleretanglewherethelowersidehasbeenreplaedby
a osine-like ar, see Figure 4, and ompare the L
2
-residual of the fully
onvergedsolutions,i.e.,thesolutionwegetwhenthexedpointiteration
doesnotimprovethesolutionanymore.
Gauÿ Greville
N r
N
:=kL(x)k
2 r
N =2
=r
N r
N
:=kL(x)k
2 r
N =2
=r
N
10 4.832e+2 4.757e+2
20 3.138e+2 1.54 2.916e+2 1.63
40 1.233e+1 2.55 1.480e+1 1.97
80 6.146 2.01 5.939 2.49
160 2.570 2.39 1.988 2.99
Figure 5. Convergenebehavior ofthe osinetestase.
Firstofallweobservethat theGauÿolloationshemedoesnotprodue
better onvergene faster than the Greville olloation sheme, although
theorypreditsfourthorderonvergenefortheformerandonlyseondor-
derforthelatter. (WeanindeedobservetheseratesforthelinearLaplae
the Greville sheme. One reason for this might be the very bad ondi-
tion of the olloation matries whih typially goes with N o
where o is
the order of the sheme. However,from the geometri pointof viewthe
resultis notentirelysatisfatory either. Asis wellknown, harmonigrid
generation systems have the tendeny to push away the grid lines from
onaveboundaries. Espeiallyifonetriestoapplytheolloationsheme
withveryoarseontrolnetsthistendeny seemstobeevenaentuated.
Figure 4,for example,showsthe resulting1010 ontrol pointgrid and
theorrespondinggridevaluation. Thisdefetdisappearsin theourseof
extensivegridrenement, but only veryslowly. Forinstane, in theeven
moreextremeexampleofFigure6oneneedsmorethan160ontrolpoints
ineahoordinatediretionbeforethegridlinesstarttoonvergetowards
theboundary. TheFigureitselfshowstheresultofthedisretizationbased
on4040ontrolpoints.
Figure6. Disretization with4040ontrol points
Atthisplaeitis interestingto note,that thesolutionof thetransformed
Laplaeequation(14)analsobeusedtoomputeforanygivengridmap-
pingx(u;v)aorrespondingontrolmappings(u;v). Thisanbedoneby
replaing theNeumannboundaryonditionsby standardDirihletondi-
tions. Figure7showstheontrolmappingthatorrespondstothegridsin
Figure6. Theaberrationofthisgridfromtheidentimappingisobviously
ausedbythedisretizationerror.
Conclusion
The here proposed olloation sheme presents a useful method to real-
ize ellipti grid generation methods diretly on B-spline representations.
FollowingSpekreijsethiswouldaordtheomputationofappropriateon-
trolmappings viasolutionofthebiharmoniequation. This, however,has
notyetbeenimplementedforB-spline grids. Sinein [2℄Mankedoesnot
reportsimilarproblems,itmightwellbethatmethodswhihomputethe
ontrol funtions iteratively are a reasonable alternative in this ontext,
too.
Acknowledgments
ThisworkhasbeenandperformedwithfundingbytheDeutsheForshungs-
gemeinshaft in theCollaborativeResearh Center SFB401 "Flow Modu-
lation and Fluid-Struture Interation at Airplane Wings"of theRWTH
Aahen,UniversityofTehnology,Aahen,Germany.
References
[1℄ StefanSpekreijse.ElliptiGridGenerationBasedonLaplae-Equations
andAlgebraiTransformations.J. Comp. Phys.,118:3861,1995.
[2℄ J.W. Manke. A Tensor ProdutB-Spline Method forNumerialGrid
Generation. J. Comp. Phys.,108:1526,1993.
[3℄ FrankBramkamp, Philipp Lamby, and Siegfried Müller. An adaptive
multisalenitevolumesolverforunsteadyandsteadystateowom-
putations. J.Comp. Phys.,197:460490,2004.
