Using Simulated Annealing to Obtain Good Nodal Approximations of Deformable Bodies

14 Pages·2007·0.54 MB·English

Using Simulated Annealing to Obtain Good
Nodal Approximations of Deformable Bodies
OliverDeussenandLeifKobbeltandPeterTu¨cke
InstituteforOperatingandDialogSystems,UniversityofKarlsruhe
76128Karlsruhe,Germany
Net:[oliver|kobbelt|tuecke]@informatik.uni-karlsruhe.de
WWW:http://i31www.informatik.uni-karlsruhe.de/[˜oliver|˜kobbelt|˜tuecke]
Abstract.Inthispaperwepresentamethodtoobtaingoodapproximationsofde-
formablebodieswithspring/masssystems.Aniterativealgorithmbasedonvoronoi
diagramsisusedtogetagoodmassdistribution.Theelasticpropertiesofthesystem
areoptimizedbysimulatedannealing.Resultsareshown,andsomeapplications
arediscussed.
Keywords:
Simulation,Spring/masslattice,Modeling,DeformableBodies,ComputerGraphics
Realtimesimulationandinteractivecontrolofdeformable bodiesareofgreatinterest
tocomputergraphicsandotherfieldslikepatternmatchingormedicalvisualization.
Spring/mass (or nodal) systems are often used to obtain efficient simulations of
deformableobjectsiflowerrequirementsforthesimulationprecisionaregiven.
Terzopoulos and Waters [15, 13] used them for modeling and simulating facial
tissue. Simulationof woven clothes was doneby Breen [4]. In an early work, Miller
[11]usednodalmodelstosimulateutricularobjects.
Inmedicalresearchsimilarsystemsareusedtosimulatehumanorgansfortraining
tasksinendoscopicsurgery[10,6],andforgeneralsurgerysimulationandplanning[5].
To reproduce a specific mechanic behaviour of a spring/mass system, parameters
likespringconstantscaneitherbepreset,oroptimizediteratively.Weshowthatwithin
somegenerallimits,suchoptimizationsyieldsystemswhichapproximateobjectswith
givenelasticpropertiesquitewell.
1 NodalApproximations
In general, a spring/mass system (SMS) is a collectionof points(or particles) with a
specificmass,thatareconnectedbysprings,dampersorplasticelements.Thediscrete
Lagrangianequationofmotionforapoint atposition andmass hasthesimple
pi ui mi
form(cf.[15]):
¨ ˙ X 1
e
miui+(cid:13)ui =F (cid:0) Fij; ( )
j2Ni
where is ageneral dampingcoefficient, and denotes an external forceactingon
e
(cid:13) F
thepoint.Theforces areexertedbyelastic,viscoticandplasticconnectorsbetween
andneighbourpoinFtsij ( :setofneighbours).
pi pj Ni
1.1 FiniteElementFormulationofNodalSystems
Theideaofapproximatingadeformablebodywithasystemofspringsissurprisingly
old.In1868,Kirsch[9]derivedthefundamentalequationsofelasticityfromasystemof
springsapproximatingasmallcube.Hisworkwasoneofthebasesofthefiniteelement
method(FEM).
Today pinjointedtrusses (assemblages of axiallyloaded elasticbar elements) are
usedinstructuralmechanicsasidealizationsformanyproblems.TheFEMapproachis
widelyusedforsuchsystems.
Thosetrussescanbeseenasspring/masssystems.Theonlydifferenceisthatwithin
trusses bars have a mass and jointsare massless whereby ina SMS massless springs
areconnectedbymass pointsthatbehavelikesphericaljoints.Inbothcases, themass
distributioncanbeexpressedbyamassmatrix .
M
Ifonlyelasticandviscoticconnectorsareused,equation(1)canbewritteninmatrix
calculus(cf.[2])by:
¨ ˙ 2
e
Mu+Du=F (cid:0)Ku: ( )
In this equation, is the vector of all point
1 1 1 T
positions1, is tuhe=dam(upinx;gumya;turixz;a:n::d;unx;udennyo;utenszt)he external forces acting on
e
D F
thenodes. isthesystemsstiffnessmatrix.Thismatrixdescribestheforcesonnodes
K
asaresultofdeformationsandcanbeseenasageneralizedspringconstant.
Incomparisontoothersystemsoffiniteelements,nodalsystemsdohavetheadvan-
tageof beingdiagonal,whichissimpleandconvenient formathematical handling
M
(cf. [8]). In addition, equation (2) can be solved by local application of equation (1)
to each node. This leads to more efficient algorithms (cf. [15]), especially if parallel
computersareused.
Table1showssomecomputationrates(calculationspersecond).Thelocalsimula-
tionmethodbasedonequation(1)isindependentofthebandwidthof .FEMmethods
K
basedonequation(2)dohavelowcomputationratesif isrecalculatedineverystep.
K
ThishastobedonetomaketheFEMmethodcomparablewiththelocalalgorithm.If
K
isupdatedeveryeighthstep,thecomputationratesforsystemswithasmallbandwidth
areclosetothatofthelocalmethod.
Forsmalldeformations(thegeometriclinearcase) canbeassumed constant.In
K
thiscasethecomputationratesarehigh,butsimulationprecisiondecreases.
The calculation of equilibriumstates is computationalexpensive because a linear
equationsystembasedon hastobesolved.
K
1Morepreciselyudescribesthedisplacementsofpointstoaninitialpositionwithzeroinnerenergyof
thesystem(indicatedby Ku 0).
j j=
2
Table1:Computationrates(calculationspersecond,SUNSparc20)
Bandwidthof 49 37 25 21
K
1)Localintegration 37 37 37 37
2)FEM( updatedeverystep) 8 9 10 10
K
3)FEM( updatedeveryeighthstep) 19 23 31 33
K
4)FEM( constant) 25 32 43 48
K
5)Equilibriumcalculation 2.5 4 7 9
Investigatedwasasystemwith84nodes,231edgesanddifferentbandwidthofK.1)computation
withlocalmethod;2)FEMsolutionwithupdatingofmatrixKineverystep;3)FEMsolution,K
isupdatedeveryeighthstep;4)FEMsolutionwithoutupdating;5)computationofequilibrium.
Asaconsequence,weuseafiniteelementmethodbasedonbandmatricestoobtain
equilibriumstateswithintheoptimizationprocess(insection5.3)andlocalmethodsfor
theefficientsolutionofequation(2)duringsimulation.
2 OptimalElastic Systems
Given adeformablebody withmass density , surface and elasticproperties
B (cid:26)(u) A
representedbyafunction (cf.[7]).Thisfunctiondeterminesthestress-strainrelation
G
3
(cid:27)=G((cid:15)); ( )
andthustheforcesneededtodeformabodyinacertainway(seesection5.1).
Inthemostgeneralcase, wehave 9 6.Forlinearproblemsandelastic
G 2 IR (cid:2)IR
homogeneousandisotropicmaterials, isasymmetric6 6matrix.Inthisformthe
G (cid:2)
relationshipisknownas Hooke’s law.Formany naturalmaterialsinthelinearelastic
case, can beexpressed bytwoindependentconstants,usually (Young’smodulus)
G E
and (Poisson’sratio).
(cid:23)
Theoptimizationprocesswillbedemonstratedintwo-dimensionalspacewithlinear
elastic loads of homogeneous materials. Extensions to 3-D and to anisotropic and
inhomogeneousmaterialsarediscussedbelow.
Togetanoptimalspring/masssystem,wehavetoapproximatethemassdistribution
bydiscretemasspoints,andelasticitybythetopologyoftheconnectionstogetherwith
theirspringconstants.Thisisdoneintwosteps.First,wefindpositionsandmassesof
thepoints,thenweinsertandoptimizetheconnections.
3 First Step:PositioningPoints
Givenanarbitrarybody withsurface ,twowaysofpointpositioningweretested.
B A
Regular grids(see below)withpredefined connectionstructures wereused as wellas
systemsbasedonrandomlygivenpoints.
3
3.1 TheIteratedVoronoiApproach
Our observation was, that systems with arbitrary given points do not lead to good
resultsinsimulatingelasticbehaviourofhomogeneousbodies.Tosolvethisproblemwe
developedamethodofmovingthepointsslightlyinordertoobtainamorehomogeneous
pointpositioning.
Foreachgivenpoint ,thecorrespondingVoronoiregion isdefined
by pi 2P V(pi;P)
3
V(pi;P)=fjq(cid:0)pij<jq(cid:0)pjj;q2IR ;pj 2P (cid:0)fpigg:
Thesetofresultingregionsforall isatesselationoftheplane,itcanbecomputed
in log time for given pointspiin two-dimensionalspace (cf. [12]). Since the
O(n n) n
boundaryofeach ispolygonal,wecaneasilyintersectitwith .
V A
(a) (b) (c) (d)
Figure1:Pointpositioning-a)randompoints;b)Voronoiregions;c)cuttedregionswithmoved
points;d)regionsafter15iterations,togetherwithpointmotion.
Next,the aremovedtothecenterofgravityofthecorrespondingVoronoiregion
(Figure1).Thpiisprocesscanbeiterated.Duringiteration,thestandarddeviationofthe
pointtopointdistancesreducesaswellasthestandarddeviationoftheVoronoi-region
sizes(Figure2).
m] 0.10 2m] 0.02
standard deviation in [ 00..00050 10 20 30 standard deviation in [ 00..00010 10 20 30
Iterations Iterations
Figure2:Standarddeviationduringiteration-a)pointtopointdistances;b)sizesofVoronoi-
regions.
Bothresultsindicateamorehomogeneouspointpositioningaftersomeiterations.
Afterthelastiterationtheweightsarecomputedaccordingtosection4(seebelow).
4
3.2 RegularGrids
Besides random positioning of points, we used regular grids together with specific
connection topologies.These configurations lead to good approximationresults with
theinvestigatedloads (see below), butmay have an anisotropicbehaviourbecause of
theirregularity(likecrystalsdo).
Figure3:a)regulargridsused;b)Suitabletopologiestoconnectthegrids.
To combine regular grids with arbitrary surfaces, tesselation is necessary. In the
contextofourquantitativestudyinthefollowingparagraphs,weomitsuchtesselations
becausetheirresultsarehardtointerpret.
4 SecondStep:CalculationoftheMass Distribution
If thepointpositions are known, we have to approximatethe mass distributionof
byassigningmasseusi : tothepoints.The approximationof bya set
Boftupels isconmsi(duerie)d=tombeigood,ifthemass moments uptoBordertwo
coincide.(Tuhie;smeim) oments determine the linear and angular accelpejrkaltions to external
forcesandtorques.
Therelevantmassmoments of arecomputedby
pjkl B
Z
2 4
j k l
pjkl = uxuyuz(cid:26)(u)du j+k+l (cid:20) : ( )
B
Substituting byafinitesetofpoints correspondstoapproximatingthisintegralby
thesum B ui
Xn 2 5
j k l
pjkl = uxiuyiuzimi j+k+l (cid:20) : ( )
1
i=
Equation(5)canbeconsideredastheapplicationofacubatureformulafortheintegrand
function overthedomain ,wherethepoints arethesample
j k l
points anfd(u) =areuxthueywuezi(cid:26)g(hut)coefficients for thBe formula. Hence, thueiapproximation
ofmimeetsallrelevantmassmomentsifandonlyifthiscubatureformulahas
a(utlie;amsti)quadBraticpolynomialprecision,i.e.,ifandonlyif
6
Am = p; ( )
where is the vector of mass coefficients,
1 1 T 000 100
m = (m ;:::;mn(cid:0) ) is the10-dimensionalvecptor=of(epxac;tpmo-;
010 001 110 101 011 200 020 002 T
p ;p ;p ;p ;p ;p ;p ;p )
5
mentsof ,and isthe 10 -matrixwith
B A ( (cid:2)n)
0 1 2 2 2 1
1 ux1 uy1 ::: uy1uz1 ux21 uy21 uz21
AT =BBB ... ux...2 uy...2 :.:..: uy2...uz2 ux...2 uy...2 uz...2 CCC:
@ A
1 2 2 2
uxn uyn ::: uynuzn uxn uyn uzn
Thesolutionoftheunderdeterminedsystem(6)withleastEuclidiannorm isof
2
theform with 10[3].Hence,themasscoefficients caknmbekfound
T T
A x x2IR m=A x
bysolving
7
T
AA x = p: ( )
Since isaVandermondematrix, isregularifthereisacollectionof10points
T
A AA
which does not lie on a polynomialsurface 0 of degree less than or
euqiualtotwo. p(x;y;z) =
Capturingthirdorhigherordermomentsof yieldsnosignificantimprovementof
B
theapproximationund thusthedegrees of freedom for thechoiceof the are used
in a way that is minimal. In this case themasses are distributedasmuiniformas
2
kmk
possible.
5 ThirdStep:OptimizingElasticity
Given pointswithpositionsand masses. The connections can be found as follows:
n
Ifthepointsaredistributedrandomly,weuseaDelaunaytriangulationtoobtainpairs
of points to be connected by springs. Two points are joined if they have a common
edgeintheVoronoidiagram.Thistriangulationmethodiswidelyusedinfiniteelement
mechanics.
For regular grids,we tested several topologiesin order tofind suitableelasticap-
proximations.GoodresultsaregivenbystructuresasthoseshowninFigure3(b).
F F
F
F
a a a a
a a a a
(1) (2) (3) (4)
Figure4:Testconfigurationswithresultingdeformationsofasurface.
Agoodapproximationoftheelasticpropertiesisfoundifthesystem’sdisplacement
due to some basic loads is similar to that of the given body . Figure 4 shows the
B
configurationsused, isasquareplate. Insituations(1)and (2),stretchingloadsare
B
applied,shearingisdoneinsituations(3)and(4).
6
5.1 GettingReferenceValues
For some simple bodies and special configurations it is possible to get an analytic
solutionoftheresultingdisplacementsduetoexternalforces(cf.[2]).2
As a reference, we calculate thedeformation of thesquare plate shown in Figure
4(1).Inthissituationauniformdistributedforceisappliedtooneside,theoppositeone
isconstrainedin -direction,butnotin -direction.
x y
σ −νσy
y y E
σx
x x E
(a) (b)
Figure5:a)Platewithuniformdistributedloadtotwooppositesides;b)resultingdeformation.
This situation is equivalent to the application of a uniform distributed load to
(cid:27)
twooppositesides, whereby thebodyisconstrainedsuch thatrigidbodymotionsare
permitted(Figure5).
We assume a linear elastic deformation and a material that can be expressed by
Poisson’sratio andYoung’smodulus .Insection2wegavethegeneralformulation
(cid:23) E
of the strain-stress relation (see equation (3)). This relation, now written in matrix
calculus,hastobesatisfied:
0 1 0 1 0 10 1
(cid:15)x 1 (cid:0)(cid:23) 1 (cid:23) 0 (cid:27)x 8
@ (cid:15)y A= E@ (cid:23)0 (cid:0)0 (cid:23) 2 1 A@ (cid:27)y A; ( )
(cid:13)xy ( +(cid:23)) (cid:28)xy
aswellastheequilibriumcondition
0 1
@0@x 0@ @@@y !@ (cid:27)(cid:27)xy A=(cid:18) Fx (cid:19) (9)
@y @x (cid:28)xy Fy
andthestrain-displacementrelation
0 @ 0 1 0 1
@ @0x @@y A(cid:18) ux (cid:19)=@ (cid:15)(cid:15)xy A: (10)
@ @ uy
@y @x (cid:13)xy
Sincerigidbodymotionsarepermitted,wehave 0 0andassume
Fx = ;Fy =
0 0
(cid:27)x =(cid:27); (cid:27)y = ; (cid:28)xy =
2Formorecomplicatedbodiesandsituationsmoreprecisefiniteelementmethodsareusedtogetreference
values.
7
whichsatisfies equation(9). Thestresses can becomputedby applicationof Hooke’s
law(equation(8)):
(cid:27) (cid:27) 0
(cid:15)x = ; (cid:15)y =(cid:0)(cid:23) ; (cid:13)xy = :
E E
Integrationaccordingtoequation(10)leadstothedesireddisplacements
(cid:27)x (cid:27)y
ux = +f(y); uy =(cid:0)(cid:23) +g(x):
E E
Weset and tozero,theproofofcorrectnessfordoingthisisomittedhere.
f(y) g(x)
5.2 InitialSMSConfiguration
TheinitialspringconstantsoftheSMSaresetaccordingtoYoung’smodulus( ).For
E
aspringofrestinglength ,thespringconstantissetheuristiciallyto 4.
0 0
s K = E(cid:3)s =
The upper part of Figure8(a) shows the calculated displacements of theinitialSMS,
approximatingthedeformationsofthegivensurface.Theyaredrawnusingdottedlines.
5.3 FastCalculationofDisplacements
Duringtheoptimizationprocess,thedisplacements oftheSMSwhichareduetothe
u
appliedforces havetobecalculatedveryoften.Withintheequilibriumstateofthe
e
F
systemtheyareobtainedbysolving
11
e b
F (cid:0)Ku=(cid:0)F ; ( )
whichisthesolutionofequation(2)withvanishing ¨ and ˙.Theadditionalvector
b
u u F
describesthebearingforcesinfixedpoints.Theequationsystemissolvedonthebase
ofbandmatrices,thenecessarynodeenumerationisdonewithanoptimizationmethod
inapreprocessingstep.Because ischangedineveryoptimizationstep,itcannotbe
K
invertedinadvance.
5.4 QualityCriteria
For each pointpositionof theSMS and each of thetest configurations,the reference
displacement iscalculatedaccordingtoparagraph5.1.TheaforementionedFEM
ref
u
approachisusedtocomputetheactualdisplacements duringoptimization.
u
One possibilityto describe the quality of an SMS approximation is the standard
deviationbetweenactualandreferencedisplacementsofallitspoints:
v
u 4
: uX(cid:0) (cid:1)2 12
Gs =t ui;ref (cid:0)ui : ( )
1
i=
Here, denotes displacements according to the test configuration of Figure 4.
i
u i
Othercriteriaareweightedstandarddeviationswhereborderandcornerpointshavea
strongerinfluence.Wedenotetheweightfactorsas ,where istheweightofinner
: :
points, thatofborderpointsand ofcornerpointGs.abc a
b c
8
5.5 OptimizationbySimulatedAnnealing
Simulated annealing(cf. [1, 14]) is used as optimizationmethod. This algorithmwas
originallydesignedfordiscreteproblemsbutwaslateradaptedtocontinuumproblems.
Themethodseemstobeadequatebecauseapuregradient-descentdeliveredbadresults,
andanalyticmethodswereveryinefficientbecauseofthehighdegreeoffreedomwithin
thesystem.Anotheradvantageisthesimplealgorithm(Figure6).
temperature:=initial temperature;
: 1034; orsomeothergoodinitialvalue :-)
Gfomriin:==0tono of temperatures-1dobeginf g
temperature:=temperature* ; over-exponential
i
r f g
good steps:=0;
k:=0;
repeat
Chooserandomlyanedge ;
e
Changerandomlyspringconstantof ;
e
Calculate accordingtoequation(12);
G
if( )thenbegin thesystemisnowabetterapproximation
G<G:=min; f g
gGomoidnstepGs:=good steps+1;
endelsebegin
Chooserandomnumber 0 1;
p2[ ; ]
if( Gmin(cid:0)G ) takeworseapproximation
p<etemperature f g
then := ;
elsepGutmbianckoGldspringconstantof ;
e
end
k:=k+1;
until(k>nlimit)or(good steps>glimit)
end
Figure6:Simulatedannealingalgorithm.
A drawback of the method is the high sensitivityto the number of temperatures,
the initial temperature, the reduction function and the number of good steps (glimit)
allowedineachtemperature.Although,wefoundgenerallyapplicablevaluesforthose
parameters(e=numberofedges):
Parameter initial temperature r nlimit glimit
Value 00001 0001 097 099 2
[ : :: : ] [ : :: : ] e(cid:3) e
Thesimulatedannealingprocessconvergesifnearlynobadstepsaredoneanymore
(duetotemperature)andtheimprovementateachtemperatureisbelowsomethreshold.
Weusedanover-exponentialtemperaturereductiontogetconvergence.
Anotherpossibilityistostopthealgorithmafteraspecificnumberoftemperatures(e.g.
no of temperatures=2*e)asdoneinFigure6.Inthiscase, thetimecomplexityofthe
algorithmis 2 .
O(e )
9
6 Extensionto3-D
Because equations (2) and (11) are formulated independently to the dimension, the
methodcanbeextendedto3-Dwithoutgeneralchanges.Insteadofoptimizingthefour
situationsofFigure4,nowwehavetoconsiderninebasicloads(Figure7).
F
F
F
(1) (2) (3)
F
F
F
(4) / (5) (6) / (7) (8) / (9)
Figure7:Basicloadsinthreedimensionalspace.
Referencevaluescanbeobtainedsimilartothe2-Dmethoddemonstratedinsection
5.1.TheannealingprocessworksthesamewayasshowninFigure6.
Some more effort is needed to calculate the Voronoi regions (now polyhedrons)
withinthemasspointiteration.ThenumberofedgesobtainedbytheDelaunaytriangu-
lationis ifworkingin2-D,but 32 inthreedimensionalspace(cf.[12]).This
O(n) O(n )
isofinterestbecausetheoptimizationhasquadratictimecomplexity.
7 Results
In the upper part of Figure 8, initial displacements of a system with 25 nodes and
connection structure 1 (see Figure 3(b)) are shown, in the lower part the system is
optimized. The constants for horizontal, vertical and diagonal springs rsp. are set to
thesamevalue,whichisconvenientforhomogeneousbodiesandreducesoptimization
effort.ThesamesystemwithallconstantsoptimizedseparatelyisshowninFigure9.
Rectangularcellsareusefultosimulateanisotropicbehaviour.Thiscan beseen in
Figure10,wheredifferentelasticitiesin and directionaremodeled.
x y
AninhomogeneousmaterialissimulatedinFigure11,theresultsareyetinneedof
improvement,butdemonstratethecapabilityoftheoptimization.
10