Designing a class library for interactive simulation of rigid body
dynamics
Citation for published version (APA):
Barenbrug, B. G. B. (2000). Designing a class library for interactive simulation of rigid body dynamics.
Technische Universiteit Eindhoven. https://doi.org/10.6100/IR531480
DOI:
10.6100/IR531480
Document status and date:
Published: 01/01/2000
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be
important differences between the submitted version and the official published version of record. People
interested in the research are advised to contact the author for the final version of the publication, or visit the
DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page
numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.
• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please
follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at:
openaccess@tue.nl
Designing a Class Library
for Interactive Simulation
Barenbrug, Bart G.B.
Designing aclass libraryforinteractive simulationof rigidbody dynamics / byBart G.B. Barenbrug. - Eindhoven: Eindhoven Universityof Technology,2000.
Proefschrift. -ISBN 90-386-0751-2 NUGI 855
Subjectheadings: animationtechniques/ computer simulation;physics / software libraries/object orientedprogramming/ computationalgeometry/ numerical algorithms
CR SubjectClassi cation(1998): I.3.5,D.2.2, D.1.5, I.3.7,I.6.8, F.2.2, G.1.0
PrintedbyUniversityPressFacilities,Eindhoven Coverdesign byBen Mobach
Coverillustrations byBart Barenbrug
Stan Ackermans Instituut,CentrumvoorTechnologischOntwerpen
c
April 2000,Eindhoven UniversityofTechnology
Niets vandeze uitgavemag worden vermenigvuldigden/ofopenbaar gemaakt doormiddelvandruk,fotokopie,micro lmof op welke andere wijze danookzonder voorafgaandeschriftelijketoestemming van de auteur.
No part of thispublication maybe reproducedortransmittedinany form orbyanymeans, electronicormechanical,includingphotocopy, recording, oranyinformationstorageand retrievalsystem, without permissionfrom thecopyright owner.
Simulation of Rigid Body Dynamics
PROEFONTWERP
ter
verkrijging
van
de
graad
van
doctor
aan
de
Technische
Universiteit
Eindhoven,
op
gezag
van
de Rector Magnificus,
prof.dr. M. Rem,
voor een
commissie
aangewezen
door
het
College
voor
Promoties
in
het
openbaar
te
verdedigen
op
donderdag 27 april 2000 om 16.00 uur
door
Bart Gerard Bernard Barenbrug
prof.dr.dipl.ing. D.K.Hammer en
prof.dr. R.M.M. Mattheij
Copromotor:
Almost six years ago, I was looking for a subject in the eld of computer graphics for a graduation project for my master's degree. Kees van Overveld suggested an assignment in the eld of dynamics in computer animation, a subject that needed some explanation by Kees, because Ihad neverheard ofit. Little didIknow thatsome sixyears laterI wouldbe writing a Ph.D. thesis on thetopic. I've learneda lot inthose six years, and for theirhelp, support and insight during the time of the project that led to this thesis, I owe thanks to manypeople.
The work described inthisthesisgot startedas a nal project fortheOOTI course. During that time, Ihad support frommyfellow(nowex-)OOTIs and Corryand Marloes, aswellas themembersoftheTechnischeToepassingengroup. The supportofthelatteronlyincreased in the last two years when I was working in that group as an AIO. I had a very nice time bothasan OOTI andasa TTer.
Ialso hada verypleasanttimewiththeMadymocrashvalidationteamat TNO,wherePaul van Bemmelen, Ronald Schouten and Erwin Poeze and the rest of the team made me feel verywelcome.
At the university, I owe thanks to Richard and Isabelle for TAO, along with the AIOs and programmers who were part of it. Gratitude also to Marja for keeping TT as vibrant a working environment as it was, to Tineke for her help during the nal period of nishing thiswork,to Maarten foralwaysbeinggamely,and toAlex, Pauland Elsaforthewonderful holidays.
Butnexttoholidays,therewasalsothepleasureofwork,whichbringsmetoallthepeoplewho colaborated on the project itself: Arnold Reusken and Mischa Courtinwhogave their input on someof thenumericalproblems. Huuband Elisabethfortheirhelp withtheintricacies of theGDP.Wim forhelpingreinvent thewheel,andKeesHuizing forhis insightson thetopic of component basedsoftwareengineering. Ginoforthediscussionsaboutinterfacing collision detectionandcollisionhandlingideas,andAlexTeleafortheinsightinthegeneralworkingsof simulationenvironments. PhapNguyenand CristianPaufortheirwork onthevisualization of forces andlookinginto Java'sEAI respectively. Themembersofthereadingcommiteefor providingthecommentsthatmadethisthesisthatwhichitis. Mostimportantly,abigthank you to Kees van Overveld for all the slight nudges and big shoves (as needed) in the right direction, forthecommentsand suggestions for improvingthisthesis and the work which it describes, formaking theproject possibleinthe rstplace, and nallyforintroducingmeto the wonderfultopic of dynamicsincomputer animationsixyears ago.
Finally,moreappreciationthancanbeputintowordstomyloveAndreaforhersupportever sincewe gottogether, and to myfamily,myparentsRiet and Ericin particular,foralltheir support over a period of a lot longer than six years in the past, and hopefully many many more yearsin thefuture.
1 Introduction 1
1.1 Previouswork . . . 2
1.1.1 Dynamics at theTUE . . . 5
1.2 Organization . . . 5
2 The environment and requirements 6 2.1 Simulation/Animationenvironments . . . 6
2.2 Theassumptionsaboutthehost system . . . 7
2.3 TherequirementsofDynamo . . . 8
3 The design 11 3.1 Methodology . . . 11
3.1.1 Views and notation. . . 11
3.1.2 Software designphases . . . 11
3.1.3 Objectorientation . . . 12
3.2 Forwarddynamics . . . 13
3.2.1 Physics . . . 13
3.2.2 The forwarddynamicsclasses . . . 16
3.2.3 Integration ofthemotion equations . . . 19
3.2.4 Forwarddynamics classesrevisited . . . 23
3.2.5 Summary . . . 27
3.3 Controllers . . . 28
3.4 Inverse dynamics . . . 32
3.4.1 Introduction . . . 32
3.4.2 The interfaces oftheconstraints and theconstraint manager . . . 34
3.4.3 The inverse dynamicsalgorithm: overview . . . 37
3.4.5 Determining @C @R . . . 42 3.4.6 SolvingR from C= @C @R R . . . 50
3.4.7 Positionalconstraints issues . . . 54
3.5 Additionalfeatures . . . 56
3.5.1 Visualizingforces . . . 56
3.5.2 Friction . . . 58
4 Specializations of the constraints and the controllers 61 4.1 Speci c constraints . . . 61
4.1.1 The empty constraint . . . 61
4.1.2 The point-to-point constraint . . . 62
4.1.3 The point-to-curve constraint . . . 63
4.1.4 The point-to-surface constraint . . . 66
4.1.5 The line-hingeconstraint . . . 67
4.1.6 The cylinderconstraint . . . 69
4.1.7 The connectorconstraint . . . 70
4.1.8 The orientationconstraint . . . 71
4.1.9 The prismconstraint . . . 72
4.1.10 The planeconstraint . . . 73
4.1.11 The barconstraint . . . 74
4.1.12 The multibarconstraint . . . 75
4.1.13 The wheelconstraint . . . 75
4.1.14 The collisionconstraint . . . 79
4.2 Speci c controllers . . . 84
4.2.1 The springand dampercontroller . . . 84
4.2.2 The PIDcontroller . . . 84
5 Examples 87 5.1 Asimpleexample. . . 87
5.2 Thebicycle . . . 89
5.4 Thedoorspring. . . 93
5.5 Therollercoaster . . . 94
6 The design process 101 6.1 Introduction. . . 101 6.2 Designphases . . . 101 6.2.1 Interfacedesign . . . 102 6.2.2 Algorithm design . . . 103 6.2.3 Code design . . . 104 6.3 Planning. . . 104
7 Conclusions and future work 106 7.1 Conclusions . . . 106
7.2 Future work . . . 107
7.2.1 Extensions . . . 107
7.2.2 Internalchanges . . . 108
A Dyna's partial derivatives 109 A.1 Introduction. . . 109
A.2 Partial derivativesfor centralforces. . . 110
A.2.1 Determining @p t+h @Ft analytically . . . 110 A.2.2 Determining @p_ t+h @F t analytically . . . 110
A.3 Partial derivativesfor torques . . . 111
A.3.1 Determining @p t+h @M t analytically . . . 111 A.3.2 Determining @p_ t+h @M t analytically . . . 112 A.3.3 Determining @d t+h @M t analytically . . . 113 A.3.4 Determining @ _ d t+h @Mt analytically . . . 114
A.4 Partial derivativesfor (non-central)forces . . . 114
A.4.1 Determining @p t+h @f p 0 t analytically . . . 114 A.4.2 Determining @p_ t+h @f p 0 t analytically . . . 114
A.4.3 Determining @d t+h @f p 0 t analytically . . . 114 A.4.4 Determining @ _ d t+h @f p 0 t analytically . . . 115
A.5 Partial derivativesfor impulses . . . 115
A.5.1 Determining @p_ t + @ t analytically . . . 115 A.5.2 Determining @ _ d t + @ t analytically. . . 116 A.5.3 Determining @p t+h @t analytically . . . 116 A.5.4 Determining @d t+h @ t analytically . . . 118 A.5.5 Determining @p_ t+h @ t analytically . . . 118 A.5.6 Determining @ _ d t+h @t analytically . . . 119
B Discrete time motion integration 120 C Velocity terms 124 D The tensor of inertia and one dimensional dynas 126 E List of symbols 127 F OMT notation 128 G Dynamo classes overview 129 H Constraint satisfaction 130 H.1 Theconstraint satisfactionproblem . . . 130
H.2 Domainspeci cissues . . . 130
H.3 Iterative constraint correction . . . 130
H.4 Constraint correction considerations . . . 131
I Source of an example 132
References 134
Summary 140
Samenvatting 141
The success of motion pictures such as Jurassic Park, The Phantom Menace and A Bug's Life show that computer animation is becoming more popular than ever. While the big movie studios such as ILM, Disney and Pixar can a ord a lot of time, money and e ort to get every motion right in these movies, in general it is very hard to make convincingly looking animation, especially animation supporting interactivity. Interactive animation is also becomingmoreandmorecommonplaceinareas suchascomputergames,simulationand virtual realityapplications.
Oneofthereasonsthatitishardtogenerateconvincinglylookinganimationsisthesuitability ofthetechniquesthatareoftenusedtospecifyandcalculatethemotionsofthemovingobjects in computer animations. Commonlyused techniquesare keyframing and kinematics. Using keyframing, an animator has to specify the exact positions and orientations of some of the objects he/she wants to animate for several moments in time, and the computer is used to calculate the positions and orientations of the objects in the in-between moments. Using kinematics, an animatordirectly speci esrotations and translations ofthe objects (possibly relativeto eachother)to specifytheirmotions. Boththese techniquesrequireaskilledartist to make themotionslooknatural,butrequirerelativelylittlecomputationale ort.
Dynamics
Withthecontinuouslygrowingpowerofcomputers, thetechniqueof usingdynamics calcula-tionshasbecomefeasibleforapplicationininteractiveanimation. So-calledforwarddynamics entails programmingtheanimationprograminawaythatitcan simulatethephysics(asfar as motions are concerned) that we see around us every day, thereby making the motions lookingnatural. In other words, dynamics calculations result in motionsby calculating the e ects of inertiaandof forces.
Furthermore,itispossibletonotonlycalculatethee ects offorces, butalso to derivemany of theforcesthat occur ina given system, based on some notionof what theire ect should be (inverse dynamics). Usingso-called constraints, one can easily specifyforexample joints betweenmovingobjects(pin-joints,ball-jointsandsuch),causingthemtoformmorecomplex articulated bodies. Such joints partially limit the relative motions of the objects that are connected, and as a result the motions of the whole are far more complex, yet can still be calculated automaticallyvia theoccuringforces.
Software structure
To use themore complex dynamics algorithmswithina software structure, attention hasto be given to theway thealgorithms are incorporated into the software. Ofmain importance is themannerinwhichthedi erentpartsofthesoftwareinterfacewitheachother. It should bepossibleto useandsteerthedynamicssubsystemwithoutbeing botheredbymanyofthe
detailsoftheunderlyingphysicsandmathematicsinamannerwhichcomesasnaturalasthe resultingmotionswilllook.
The main goal of the work presented in this thesis is to provide a software structure for performingand steeringthedynamicscalculations. Emphasisisnotonlyonthecalculations themselves, but also on how they are distributed over the components of a well structured pieceofsoftware. ThissoftwareismadeavailableintheformofaclasslibrarycalledDynamo (the DYNAmicsMOtion library,see[Baren 99 ]).
Software design
Determining a good software structure and interface is the main design challenge in the creationofDynamo. Makingthecomplexdynamicsalgorithmsaccessible,andharnessingthe approximately 20.000 lines of source code of Dynamo in a structure that is understandable and extendible takes a lot of extra e ort next to the algorithm design which would be the focus of a regular Ph.D. thesis. This is why this thesis has taken the form of a designer's Ph.D.thesis.
Bene ciaries
Thebene ciariesofthisworkare rstofalltheusersofDynamo. Dynamocanbeusedinthe eld of computer animation,but because of thesimulation approach taken, the application area of the software described in this thesis is widenedfrom just application in the eld of computer animation: since themotions are (close to) physically correct, it can also be used to analyse mechanical systems. For a more detailed description of the pro le of a user of Dynamo,seeSection2.1.
Secondofall,peopleinterestedinsoftwaredesigncanbene tfromthisworkasanexampleof how arelativelylargesystem suchas Dynamocan bedeveloped ina way thatuses thetools provided by object orientation to keep its structure clear, and its interface understandable. PartofthisworkhasbeensponsoredbytheTNORoad-VehiclesResearchInstitute,wherethe Madymosystem is developed and exploited. Madymo provides (non-interactive) simulation ofmotions, andis aimedtowards reliablesimulation,sosafety issuesregardingroadvehicles can be accurately studied. Madymoalready incorporatesdynamics calculations, soDynamo itselfisnotsomuchofinterestofTNO,butthesoftwaredesignaspectsoftheworkdescribed hereare.
1.1 Previous work
Forseveral decades,dynamicscalculationshave beenappliedinthe eldofsimulation. Only recently,thecomputingpowerhasbeenavailabletoapplydynamicsininteractivesimulations. Several distinctionscan be madewith respect to software inthe eld ofmotion dynamics.
(i) Modular versus monolithic architectures
The existingsystems di erconsiderablyinarchitecture: inmostsystems,thedynamics soft-ware is embedded within a speci c simulation or animation program (the dynamics within Madymo is such an example), and together with that program forms a rather monolithic structure, whichcan give problemswithmaintenance andextendibility. Inrecentyears there hasbeenatrend totry todevelopsoftwareinamore modularfashionintheformof compo-nents. Fromthat pointofview,apieceof softwarewhichprovidesdynamicscanbeseenasa component ([Huiz 97 ]). This requiresa more modulararchitecture notonlyof thedynamics software, butalso of thesoftwaresysteminto whichit isembedded.
(ii) Rigid moving objects versus exible moving objects
Most methodsassumerigidnessof theobjectsthat have dynamicmotion. But work hasalso beendoneon exiblemodels,for examplein[Platt88 ].
(iii) Constraints restricting the degrees of freedom of the moving objects
Thewaythedegreesoffreedomofthemovingobjectsarerepresentedisalsoamajordi erence betweenexisting solutions. A collection of n moving rigidbodieshas6n degrees of freedom (a positionand anorientationforeach rigidbody). Supposemdegrees offreedom areto be restrictedbyconstraints.
The dynamics and inverse dynamicscalculationscan be combined byexpressingthedegrees of freedom of an object in well-chosen generalized coordinates (as was done for Madymo, see [Wijck96 ]). Inthat way there areonly6n m variables inthesystem to be calculated. Thisapproachworkswellinsituationswheretheconstraintsconnecttherigidbodiesin tree-like structures (loops aregenerally notsupportedin thisapproach), and insituations where the connectionsbetweentherigid bodiesremain the same overthe courseof the animation. This approachis oftenused inrobotics.
A more general approach is to represent the degrees of freedom of each object in the 6n cartesian coordinates. There are then two ways in which the constraint problem can be solved.
The rst is to impose the constraints by directly taking the algebraic constraint equations into account when solving the forwarddynamics problem([Haug 89 ]). This way, a problem of 6n variableshas to be solved under theconditionthat m constraintshold.
Thesecondistointroducemextravariablesfortheconstraintforceswhichrestrictthedegrees of freedom of the movingobjects. First, the m constraint variables are calculated, and the resulting forces arethen inputforthe calculationsthat evaluatethe 6nrigid body variables over time. This decouples the forward dynamics calculations from the inverse dynamics calculations, astheconstraint forcesarejustthesameto theforwarddynamicssubsystemas anyother externally suppliedforces.
Note thaton thelatter casethem variablescorrespond to therestricteddegreesof freedom, whilethevariablesinthe very rst approach correspond to thefree degreesof freedom.
Thedecouplingoftheforwardandinversedynamicscalculationsinthelatterapproachallows fortheuseofthesameinversedynamicsmethodologyincasetheforwarddynamicssubsystem isreplaced byforexamplea subsystemwhich allows for exible(non-rigid)objects.
(iv)Accuracy versus speed
Forapplicationinsimulation,themainfocusisontheaccuracy ofthedynamicscalculations, sincesuchsimulationsareforexampleusedtogather knowledgeaboutthesafetyof mechan-ical systems (Madymo). Forinteractive animation, accuracy can sometimes be sacri cedin exchangeforspeed. Afastresponsecan insuch applicationsbe farmore important than ex-actness, aslongasthe motionsstill \looknatural". A major factor inthisregardistheway constraintsaresatis ed. Thiscan bedonethroughapplicationofconstraintforces(which ac-curatelymodeltheunderlyingphysics),orby`fakingit'bydirectlycalculatingdisplacements and rotations([Gasc 94], anapproachwhichispartlybasedon [Overv 91 ]). Theconnections (otherwise modeled byconstraints) can be also be approximated by spring damper systems ([Kell 96]). For very approximative methods, [Chenn98 ] describeshow only themotions of objectswhicharevisiblecan becalculated, usingstatisticalmodelsto estimate theirmotion when they reappear.
(v) Impulse-based versus force-based
Some approaches take motion discontinuities such as collisions as a starting point ([Faure 96 ], [Faure 99 ]) instead of (continuous) constraints. This can give rise ([Mirti96 ]) to an impulse-based approach in which objects follow a ballistic trajectory in between im-pulseexchanges. Thetimingmodelindicedbythisapproachisnotverysuitableforinteractive animations,becauseofthelargevariationsinthesizesofthetimestepstaken,whichdoesnot correspond to theconstant timebetween framesof an animation. Impulse-basedapproaches require extra constructs to supportarticulated rigid bodysystems, because the rigidbodies then donot move ballisticallyinbetweencollisions.
(vi) Di erentiating the constraints or not
Constraints give rise to algebraic di erentialequations (ADEs). These can be solved using numerical methods, or by di erentiating them ([Baraf 92 ], [Baraf 96]) to arrive at easier-to-solve di erential equations. The latter approach su ers the risk of numerical drift, which can cause problems for longer-running animations/simulations. Such drift then has to be compensated(see forexample[Baum 72 ]).
1.1.1 Dynamics at the TUE
At theTUE,research towards addingdynamicsto thearsenalofknownmotionspeci cation and calculation tools has started quite a while ago. Loosely based on the model described in[Barz88 ](whichwasproposedforassemblingarticulatedbodies,andnotanimatingthem) the ideas described in [Overv 93a ], [Overv 93b] and [Baren 94 ]) have been implemented in the animation program Walt. The successor of Walt is the GDP: the Generalized Display Processor, a system which uses the direct modi cation and direct manipulation paradigms to allow applicationsand/or usersto interact ina 3Denvironment (see [Peet 95 ]). Some of the results pertaining to dynamics obtained from Walt have been introduced in the GDP (see [Krijg94]), but by far not all, and more recent studies (see [Overv 95 ]) have shown that the algorithms used (relaxation, Euler integration, see [Baren94 ]) should be replaced bybetter ones. The Walt implementationwas also a monolithicapproach where hardly any attention waspaid to theinterfaces,resulting incode thatwas notvery extendable.
Thesewereamong the reasons forthework described inthisthesis. It hasbeencarried out intwo parts: the rstpartcomprisedthecreationof thebasicdesign,and implementingthis aspartoftheGDP system. Thiswasdoneintheformofa nalprojectfortheOOTIcourse (the postgraduate programme for Software Technology at the TUE/SAI), and is described in[Baren 96 ]. Inthesecondpart,thiswork wasexpandeduponbyseparatingthe implemen-tationfromtheGDPandprovidingitasastandalonelibrary(makingthesoftwareusablefor other programs),andaugmentingitwithnewfeaturesandupgradingtheinternalalgorithms for moreversatilityand robustness. The attention towards thedesignaspects isalso evident in thecooperation withTNO.
1.2 Organization
Theremainderofthisthesisisstructuredinthefollowingmanner. InChapter2ofthisthesis, animation/simulationsystemsingeneralaredescribed,andtherequirementsforthesoftware are presented. The design of the software is then described in Chapter 3. In that chapter, the forward dynamics subsystemis discussed rst (Section 3.2), followed by an overview of the controllers which can be used to steer thedynamicallymovinggeometries (Section 3.3). Then theinverse dynamicssubsystemisintroduced(Section 3.4), which allowsthemodeling of more complex articulated body systems, and which can be used to perform for example collision handling. In Chapter 4, theactual constraintsand controllers aredescribed, which are speci c applicationsofthemore general theorydescribedin Chapter3. Examples which showsomeusesofthelibraryarethenshowninChapter5. Thedesignprocessissubsequently presented inChapter 6, and conclusions and recommendations forfuture work are given in Chapter 7.
requirements
2.1 Simulation/Animation environments
In [Telea 99m ] and [Telea 99j ], three categories of users of animation/simulation tools are distinguished:
The end user An enduser runsanimations and simulationsandinteractswith them.
The application designer An applicationdesignercreates applicationswhichan enduser can run. He/she does not build every application from scratch, but uses prede ned functionality: the application is assmbled from functional modules, so-called compo-nents. An applicationdesignerdecideswhichcomponentstouseand howthey interact with each otherand withthe user.
The component designer A componentdesignercreates moduleswhichprovide function-alityand services to the application designer. A component usually concentrates on a speci c type of functionality, such as rendering, or user-interfacing, or nite element analysis.
Sincecomponentshave to becombined,implicitly,there is afourth category:
The framework designer A frameworkdesignercreatesasystemwhichcanbeusedbyan applicationdesignertoeasilylinkthedi erentcomponentsdevelopedbythecomponent designers. Examples of frameworks are the GDP, the Java virtual machine, AVS and Vission([Telea99m ]).
AVS for example is a simulation framework (see [AVS]). In AVS, an application designer assembles applications by instantiating components and connecting their ports, creating a data- ow diagram which represents the application. In for example the GDP (provided by frameworkdesignerEricPeeters,see[Peet 95 ]),thecomponentdesignerprovidesthedi erent classlibraries(for rendering,userinterfacing,kinematicmotion calculation,andsoon). The applicationdesignerwrites Looksapplications, usingthe classesprovided bythe component designer and the framework designer, and deciding how an end user (who runs the Looks application)can interactwiththesimulation.
The focus of this thesis is on the component design phase. A library is built which o ers the functionality of dynamic motion calculations. As any component, the Dynamo library hasapplication designersas immediateusers, and end usersasindirect users. It also means that Dynamo library assumes the presence of a framework (also called host system) which, throughtheuseofother components,takescareofall aspects oftheanimationsthatarenot
directly related to the dynamic motion calculations (such as renderingthe moving objects, and high level control). This enables the focus of this thesis to be solely on the dynamic motion calculations. The next section focuses on the host system and what the Dynamo component expectsfrom it.
2.2 The assumptions about the host system
Dynamo aimsto providethemotionsforsomeof themovingobjectsinan animation. These objects willbecalled dynas from nowon (see Section3.1.1 forgeneralmarks about nomen-clatureand notation). Themotionofotherobjectscan becalculatedusingdi erent methods (kinematics, motion captureetc.). All moving objects are visualized by the host system, or one oftheother componentsitemploys. Thehost systemtherefore alreadyhasa representa-tionforthem,which,amongstothers,capturestheirshapes. Foreachframeoftheanimation, the output of the dynamics calculations is a new set of values for position, linear velocity, orientation, and angular velocity for each dyna. These are properties that Dynamo keeps track of. A set of these properties (position, orientation, and their rst time derivatives) is called a motion state from now on. This assumes rigid bodies: for non-rigid bodiesthe de-formations oftheobjects mustalso betaken into account andshouldbereturned asa result of thedynamics calculations.
The host system shouldprovidethe informationto initialize themotion state of each of the dynas(i.e. a functionsuch asreadinggeometryinformationfrom a le,and therebyreading the initialpositionandorientation,is theresponsibilityof thehostsystem). Fortheobjects that arenotundercontrolofDynamo,butparticipatinginone ofDynamo'sconstraints,the host system should provide the motion state for each frame to enable Dynamo to evaluate the constraints.
Of all the other geometrical and topological dataassociated witha movingobject, onlythe associated mass and mass distribution are of importanceto the dynamics calculations (see Section 3.2.1). Providing this information (based on the geometrical details of the moving objects, for which Dynamo does not know the representations) is also the responsibility of the hostsystem.
Setting the motion state and mass (distribution) for each dyna provides all the required initialization. Afterwards,thecalculationsperformed byDynamo shouldbe interleavedwith all theother tasksof thehost system (such asrendering a frame, handlinguser inputetc.). Tosupportthis,Dynamoprovidesafunctionwhichcanbecalledbythehostsystemtocause Dynamo to calculatethepositionsand orientations forthenext frame. Callingthisfunction at appropriate timesis theresponsibilityofthe hostsystem.
Since Dynamo maintainsconstraints, the host systemshould take care that no actions that can violatetheseconstraintsareperformedbetweenthecallto Dynamoand therenderingof the correspondingnew frame. Forexample: moving akinematically controlledobject which is relatedto adynamicallycontrolledobjectbyaconstraintcouldinvalidate thatconstraint. Anyobjects notundercontrolofDynamo,butparticipatinginone ofDynamo's constraints, shouldalready have themotionstate forthenewframebeforeDynamoiscalled,soDynamo can match thenew motionstates of thedynasto ful lltheconstraints.
Simulatedtime isa propertyglobal to thewholeanimationsystem. The hostsystem should decidehow much the simulation time shouldbe advanced between one frame and the next, andprovidethisinformationtoDynamo,beforemakingthecallthattellsDynamotoprovide themotionstates forthenext frame. Thiswillensure thatthetimesteps taken byDynamo match the time steps taken by the objects that are not under Dynamo's control. Dynamo providesa methodto specifythistime stepinits interface.
Thefollowingpseudocodesummarizestheresponsibilitiesandthebehaviorofthehostsystem towards Dynamo:
for eachdyna do
set initialmotionstate and massdistribution od;
whileanimating do
performhostspeci ctasks (suchashandling userinput); set thetimestep forthisframe;
updateall motionstates notunderDynamo's control, but participating inDynamo's constraints;
callDynamo;
renderthenew frame od
DuringDynamo's development,theGDP (see[Peet 95 ])wasusedasahostsystem, enabling the use of its framework and already available libraries for testing, and creating the demo examplesdescribedin Chapter5.
2.3 The requirements of Dynamo
Dynamohasbeendeveloped intwostages. The rst partwasperformedasa nalprojectof theOOTIcourse, andinthispartthebasicdesignwasmade. Inthesecondpart,thisdesign was expandedupon through theaddition of several new features and re nement of existing features.
Withrespectto thecriteriapresentedinSection 1.1(previouswork),thebasicdesign devel-oped inthe rst part wassubjectto thefollowingrequirements:
(i)a It shouldhave a modulararchitecture: i.e. itshouldnotbepart ofa monolithic struc-ture.
(i)b The designshouldallowfor exibilityand easeofextendibility,andits interfaceshould provideeasy accessto its functionality.
(ii) Rigidnessofthemovingobjects canbeassumed,althougheasyextendibilityto exible geometries shouldbe taken into account.
(iii) To support the the exibility of being able to add and delete constraints on the y and notlimitingthe constraint con guration to a tree structure, therepresentation of
the degrees of freedom should be in cartesian coordinates (as opposed to well-chosen generalizedcoordinates;seethetextunderthethirdheaderinSection1.1). Thismeans that therigid objects are all positioned inworld coordinates, and not relative to one-another.
(iv)a The design should be suitable for interactive animations: i.e. speed is very impor-tant, andassuch,accuracymaybesacri cedinexchangeformore speed(if absolutely necessary).
(iv)b However,theconstraintsshouldbesatis edusingexplicitlycalculatedconstraintforces. `Faking it' such asdescribed inSection1.1 underthe header ofaccuracy versus speed, wouldlikelynotallowforsuÆcientaccuracyforapplicationinmoresimulationoriented areas. Having explicitlycalculated constraint forces allows the inverse dynamics sub-component to make decisions based on the occurringconstraint forces (which can be interpreted in a very natural way), and for the application designer to inspect them. The availability of all forces arising in a simulated system (including those resulting from constraints and controllers)allowsfor agreater insightinthe dynamicsof such a system.
(v) The dynamics systemshouldnotbepurely impulse-driven,butfullysupportforceson the dynasto supportarticulated rigidbodysystems.
(vi) The algebraic di erential equations arising from the constraints should be solved nu-merically to avoiddriftinlong animations(suchas invirtualrealityapplications).
Furthermore,thebasicdesignshouldalso:
incorporate techniques, where ever possible, that improve upon the versatility of the earlier TUE designs, to alleviate the problems identi ed with the work described in[Baren 94 ]:
{ accuracy islimitedbythesoleuseof Eulerintegration (thisrelates tothe require-ment of exibility: requirement (i)b)
{ theearlierdesign(orlackthereof) doesnotallowadditionof newconstrainttypes without the need for modifying all existing constraint types (this relates to the requirementof extendibility: requirement(i)b)
{ the useof relaxation (a method for solving constraints locally)does not support very long chains of objects (this relates to requirement (iii), since relaxation se-riously limitssupported constraint con gurations because it can nothandle long chains).
{ theuseofconstraintstoenforce rigidityaddsalargenumberofconstraints, which are expensive to solve (this relates to requirement (iv)a: speed is an important issue forDynamo).
providetheinterfacethatthehostsystemneedsto ful llits requirements(seeprevious section).
Tofacilitaterapidprototyping, thedynamicsfunctionalitywasincorporateddirectlyintothe GDP inthe rstphase (keepinga future separationinmind). Dynamo wasuncoupled from the GDP, before the work on the actual extensions and improvements in the second phase couldcommence. Thenthefollowingextensionswererequired:
1. A user manual forthestand-alonelibrary.
2. More examples of possible applications of the library, adding constraints (such as the wheel constraint forthebicycleexample) asneeded.
3. Improvement of the collision handling provided by the basic implementation (for the collision constraint that will be presented in Section 4.1.14, this entails the addition of the positional terms that prevent the positional aliasing when collisions are usedin prolonged contact).
4. Additionof awayto modelfriction.
5. Addition of more curves and surface types for the point-to-curve and point-to-surface constraints.
6. Additionofhandlingnearmass-less objects,such asropes. Theseobjects onlytransfer forces (the total force on them isalways zero), butdo notsigni cantlycontribute any mass and therefore can not be modeled by dynas, because the motion equations (see Section 3.2.1)can notbeapplied for(near)zeromasses.
7. A mechanism for creating controllers. Controllers are a means to actively steer the dynas(whereasconstraintsareonlyreactive), seeSection3.3. A fewcontrollers should be provided asexamples. This shouldbe donein such a way that itis easy to extend the librarywithmore controllers later.
8. Increase therobustnessof theinverse dynamicssubsystembyaddinghandlingof over-speci edand under-speci edconstraint combinations.
9. Sharetheexperienceinsoftwareengineeringtechniquesbyprovidingconsultancytothe Madymo groupat TNO.
Duringtheprojecttheserequirementsweremonitoredinregularmeetingstoensurethatthey weremet to satisfaction.
3.1 Methodology
3.1.1 Views and notation
Indesigningsoftwarethat dealswithsimulatingphysics, threeviews can be distinguished:
Physics view Inthisview thereal world physics is observed.
Modeling view Inthisview a(mathematical)modeloftherealworldphysicsis described. Some detailsfrom thereal world physicsare abstractedfromorsimpli ed.
Software view In thisviewthe modelis implementedinasoftware structure.
For example, a physical object has mass. This mass can be modeled in a variable m which denotes theamount ofmass oftheobject. Inthesoftware, thisamount willlikelybe admin-istrated inanattribute massof a classwhichimplementsthemodel ofthephysicalobject.
Thesethree viewsareverycloselyrelatedandarereferred toinaninterwoven mannerinthe text of this thesis. As shown above, the physicalview is denoted inregular font, the model viewinslantednotation,andthesoftwareviewintypewriterfont. However,strictlyadhering to this notation would give a very verbose text. For example: many modeling parameters map one-to-one to class attributes. Explicitlyintroducingthese would be very cumbersome, and therefore themodelingparametersare useddirectlyinsome of thepseudo codeand the text. So, for sake of brevity, it is sometimes left to the reader to distinguish between the di erent views, given thecontext inwhich they areused.
3.1.2 Software design phases
Three designphasesareidenti edinthesoftwaredesignof Dynamo:
Interfacedesign Theinterfaceshouldprovideauniformwayforanapplicationdesignerto specifytherequiredmotionsand relationsbetweenthe objectsthat areanimated.
Algorithm design The algorithms have to be selected and (if necessary) newlydesigned, to providethefunctionalityspeci edbythe interface.
Code design The code should be structured such that the di erent classes interface with eachotherinamannerthattheycollaboratetoimplementthealgorithms. Inthisphase the interfaces thatareinternal to thelibraryare designed.
Moreinformationonthesephases, andanexplanationofhowthesephases tinto theoverall software lifecyclemodelarepresentedinSection6.2.
3.1.3 Object orientation
Thedesign process and methodologyare discussedinmore detailinSection 6,butone very important designdecisionneeds tobe mentionedhere:
Design decision 1: The designof Dynamois doneinan objectorientedfashion.
Thereareseveral reasonsforthis, some of whichare:
Objectorientationallows foracommonlanguageindiscussingandreasoningaboutthe system in all of its views (even so much that the notational problemdiscussed at the endofSection3.1.1 arises). ThemechanicalsystemsthataremodeledbyDynamolend themselvesverywelltomodelinginanobjectorientedfashion,becausetheidenti cation oftheclassesofphysicalobjectssuchasdynasandconstraintscomesverynatural. The physical`interfaces'of these classesalso map very easilyto the software interfaces (for example: sincepushingand pullinga dynais theonlywayto make itmove,a method suchasapplyforcecomes naturally;seeSection3.2.2).
Objectorienteddesignallows fora clearseparationbetweeninterfaceand implementa-tion. This allows for a higher level of abstraction at the interfacelevel, which aids in clearer (and oftensmaller)interfaces. Theimplementationbehind thoseinterfaces can be variedand improved, withouta ecting other components (changingand improving algorithms is a constant focus of attention in research-oriented design such asthat of Dynamo).
Object orientationsupports several constructs(such asinheritance and object compo-sition)to specifyrelationsbetween subcomponents(classes). Thisallows a designerto providea clearstructure betweenthe componentsof a library,asopposed to just pro-viding alargecollection offunctionsandprocedures. Usingthisstructure,abstractions can be carried throughinthedesignin all itsphases, sothat forexamplethe abstract notion of a constraint can be foundfrom requirements to implementation. This helps in understanding,maintaining and extendingthe library. Ofspecialinterest is the in-heritance mechanism, which can be used to specifyabstract interfaces. In this way, a generalinterfacespeci cationsuchas\theconstraintinterface"canbede ned,enabling a uniforminterfaceforallconstraints(ease ofunderstandingand interfacing),and pro-vidinga\recipe"forextendingthelibrarywithnewconstrainttypes(onesimplyhasto implement the methods speci edby theconstraint interface to implement a new type of new constraint). Such useof inheritanceallows forre-use ofstructure.
Manyframeworksarealsobasedonobjectorienteddesign(mainlybecause ofthesame reasons that Dynamo is designed in an object oriented fashion), making it easier to incorporate Dynamo inthese frameworks. Seefor exampleSection 5.5 foran example that uses theadvantages of object orientation next to the functionalityof Dynamo to modela more complexsystem.
The following sections address these design issuesfor the di erent components of Dynamo. First the forward dynamics subsystem is discussed. The controller and inverse dynamics subsystems then buildon the forward dynamics functionality in calculating forces that can be inputforthe forwarddynamicscalculations.
3.2 Forward dynamics
This section dealswith the designand implementationof theforward dynamics subsystem. This subsystemshouldallowmotionsofthe movingobjectsinan animationto becomputed inaccordance tophysics,i.e. based oninertiaandappliedforces. Inorderto understandthe issues,we rstneedtotakealookatthephysicalconceptsthatarebeingmodeled. Thisis pre-sentedinSection3.2.1below. Fromthephysics,wecan derivethesoftwareinterfacethrough whichthis functionalitycan be made available. This ispresented inSection3.2.2. From the equations given by physics, we can also derive the mathematical schemes needed to imple-ment theinterface(Section 3.2.3). Once these schemes areknown,we can incorporate them in a re nement step of theclasses forwhich we earlier de nedthe interfaces (Section 3.2.4). Afterwards,a shortsummaryof thepreceding sectionsis given (Section 3.2.5).
So we rst take alook at thephysicsview to derivethemodel.
3.2.1 Physics
The motion equations below are the results of a more detailed derivation presented in [Witt77 ]. Foramore comprehensive treatment ofthephysics seeforexample[Alonz 70 ].
Geometries and dynas
Westartbyintroducingthemovingobjectswhicharepresentintheanimations/simulations. We will call these the geometries (also called rigid bodies sometimes). For dynamics it is importanttoknowthepositionandorientationofageometry,andtheirderivatives: linearand angularvelocity. Thisinformationiscalledthemotionstateofthegeometry(asmentionedin Section2.2). Thepositionandorientationofageometryde nealocalcoordinatesystem. This induces the distinctionbetween local coordinates and world coordinates. Local coordinates areveryconvenient to denotepositionswithina geometry(sincethegeometries areassumed to be rigid,these remain constant), but world coordinates are required when comparingfor examplepositionsbetween geometries (asis oftendoneinconstraints).
Thepositionof theorigin ofa geometry willbe denoted byvector z,and theorientation by orthonormal rotation matrixA. Conversion from local coordinates to world coordinates p isdoneusing:
p=z+A
The linear velocity of the geometry is denoted by v, and the angular velocity by !. The velocity inworld coordinatesp_ of a point of a dynawith localcoordinates can also easily
be expressed:
_
p=v+!A
Themotionofarigidbodyisgoverned byits massm (determiningitslinearinertia),and by its tensorof inertiaJ (determiningits rotational inertia,see alsoAppendixD).
Forthemotionsequationspresentedinthenextsection tobevalid,thereisanextra require-ment on the position z of a dyna: the local coordinate system of the dyna should have its origininits centerof gravity (i.e. z needs to denotethe centroidofthedyna).
Itisalsoveryconvenientiftheaxesofthelocalcoordinatesystemarealignedwiththemajor axesofinertiaofthedyna: thiscausesthetensorofinertiato bediagonalwhenexpressedin thelocalcoordinate system. We willrequirethisproperty frominitialorientation A.
Forrigidobjects,thetensorofinertiamatrixJ isaconstant(anddiagonalaswejustassumed) when it is expressed in the local coordinate system of the object. It is variable when it is expressed in world coordinates: J
w
=AJA T
. More on the tensor of inertia for the special case of one-dimensionalobjects can befoundinAppendixD.
Inertia
When there are no external in uences to a dyna, it moves under the e ects of inertiaonly. For the position this means that z_ = v. Naturally, the angular velocity ! and the change in orientation (
_
A ) are also related. This relation can easily be expressed when orientation is represented by Euler parameters (see section 2.1.3 of [Witt77 ]). An orientation which correspondsto arotationwithrespecttotheworldcoordinatesystemoveranangle around axis u (where jjujj = 1) is represented by a unit-length four-dimensional vector q, where q 0 = cos ( 2 ) and q 1::3 = usin( 2
). Orientation matrix A can be expressed in terms of q as follows: A=
2(q 2 0 +q 2 1 ) 1 2(q 1 q 2 +q 0 q 3 ) 2(q 1 q 3 q 0 q 2 ) 2(q 2 q 1 q 0 q 3 ) 2(q 2 0 +q 2 2 ) 1 2(q 2 q 3 +q 0 q 1 ) 2(q 3 q 1 +q 0 q 2 ) 2(q 3 q 2 q 0 q 1 ) 2(q 2 0 +q 2 3 ) 1
(3.1)So it is suÆcient to know q to be able to calculate A (and vice versa). The correspondence between! and the change inorientation is more easily given in terms of therepresentation usingEulerparameters(specifyingtherelationbetween!andq)._ Introducingthe?operator), thisrelationis given by:
_ q = 1 2
0 ! T ! !~
q = def !?q~ ! =
0 ! 3 ! 2 ! 3 0 ! 1 ! 2 ! 1 0
(3.2)Forces and torques
When forcesand/or torques are beingapplied to a dyna, thiscauses changes in its velocity (linear and/or angular). Central forces F (i.e. forces applied to the centroid of the dyna) a ect thelinearvelocity,according to
mv_ =F
TorquesM a ect theangularvelocity, accordingto
J w _ !+!J~ w !=M
A non-central force f (applied to point p) a ects both the linear and angular velocity. Its e ect on the linear velocity is equivalent to a central force F = f, while its e ect on the angular velocityis equivalentto that ofa torque
1
M =(p z)f =Af.
Togetherwith theequationsforinertiamotion,theequationsof motionof a dynathenare:
_ z = v (3.3) _ v = 1 m F (3.4) _ q = !?q (3.5) _ ! = J 1 w (M !J~ w !) (3.6)
Forwardintegrationoftheseequationsovertimeyieldsthemotionofadyna. Thisisdiscussed further inSection3.2.3.
Impulse changes
Forforces working over a very shorttime (like those occuring as a result of collisions), it is often easy to work directlywith changesinimpulse:
(t)=
Z
t0
f( )d
Using thisapproach,theexact valuesof the forcesduringthetime intervalfrom 0 to tneed not be known, just theiroverall e ect on the change of impulse . This is especially useful for very large forcesworking over a very shorttime interval. In the limitwhere theinterval duration goes to zero and the forces to in nity, the impulse change remains nite. In this
1
Notethatthetorqueisderivedfromthe(non-central)force: forcesarefundamental(asusedinforexample Newton'saction={reactionprinciple), whiletorques arenot. Wewillbase adesigndecision(designdecision 7)onthislater.
case,thee ectsofnon-impulsiveforcesand inertiavanish,sothatthemotionequationsthen become (for animpulsiveforce appliedto thepoint withlocal coordinates):
mv = J w ! = A or: v = 1 m (3.7) ! = J 1 w (A) (3.8)
3.2.2 The forward dynamics classes
The geometry class
In the previous section, the physics were presented that are to be modeled in the Dynamo software. Thissection willpresent the structureand theinterfaceof thesoftwarewhich will provideeasy accessto thisfunctionality.
As mentioned in Section 2.2, host systems for Dynamo will already have some notion of a geometry(forpurposesofrenderingetc.),whileDynamowillonlyadministratetheproperties it needs incompanionobjects:
Design decision 2: Dynamowillmanagecompanionobjects tothehostsystem's geome-tries.
Callbacks are de ned,which a user must implement in order to specifyhow the companion objects can interface with the geometries from the host system. Using these callbacks, a user can tellDynamo how to retrieve andupdatethemotion state ofa geometryof thehost system.
Design decision 3: Companion objects interface with their counterparts through call-backs.
Formore informationon thecallbacks,see Section2.2, and theuser manual[Baren 99 ].
Dynamo'scompanion objectsarethereforenothingmorethanameansfortheotherDynamo classes to interfacewith thegeometries of thehost system: thecompanion objects adminis-trates and providesaccess to themotionstate(s) of thehost system'sgeometries throughan interface which is known to the other Dynamo classes. The callbacks are used to keep the motion statesof thecompanion objects and thehostsystem's geometries incorrespondence.
Since time is discrete in animation systems (in an animation, the only points in time of interest to an end user are those for which a frame is rendered, since those are the only points in time he/she can observe), there are two motion states which are of interest when calculatingthemotionsfromoneframetothenext. The rstisthemotionstateforthelatest
frame(sometimesalsocalledthecurrentframe),whichistheresultofthemotioncalculations from last frame, and which functions as a starting point for the motions of the geometries betweenthisframeand thenext. Thesecond isthemotionstate forthenext frame,whichis beingcalculated. We willdenotethetimestampofthe formerwitht, andthedurationtime betweenframeswithh (sothetimestampof thenext motionstate ist+h). The companion objectsshouldadministrate boththese motionstates.
Thereisacleardistinctionbetweentwogroupsofgeometries. Ontheonehand,therearethe dynasfor which the motions are calculated by Dynamo. On the other hand, there are also othergeometrieswhicharenotunderDynamo'scontrol(thesecanbekinematicallycontrolled forexample).
Fromtheviewofthehostsystem,therearemanydi erentkindsofgeometries,eachofwhich have di erent motion behavior. Allgeometry types have incommon that they administrate and provideaccess to theirmotionstates. Theydi er inthe way thattheirmotion behavior is modeled, and therebyin the way new motion states are calculated. Only for dynas, new motion states are calculated using equations 3.3{3.8. Other types of geometry might use kinematics calculations for example. To be able to interface with every type of geometry present in the host system, the companion objects in Dynamo are companions to the host system's geometries: a geometry class is introduced in Dynamo (from now on, reference to \geometry"willrefertoDynamo'sgeometryclass,andnotthehostsystem'sgeometries,unless speci callystatedotherwise). Dynamo's companion objects to the host system'sgeometries belong to thisgeometryclass.
Thegeometryclassprovidesmethodsforinspectingandchangingbothmotionstates,andfor conversionsbetweenlocaland worldcoordinatesbased on bothmotion states.
The dyna class
Dynamo'sdynasaddthedynamics-speci cmotionbehaviortosomeofthegeometries. There-fore:
Design decision 4: Inheritence from the geometry class is used to introduce forward dynamics.
This allows for reuse of the motion states administrated by the geometry class, and the interface through which it can be accessed. It also allows reuse of the callbacks the user de ned to communicate the motion states. The specialization of the geometry class is the dynaclass.
Thedynaclassaddsseveral attributesandmethodswithrespectto itsparentclasstoprovide access to thedynamicmotionfunctionality.
First of all, it adds attributes for administrating the required mass and tensor of inertia of the dyna,sothedynaclasscan evaluatethemotion equationspresentedinSection 3.2.1.
Second,therearemethodsthatcan becalledto applyforces,torques andimpulsechangesto a dyna:
applyforce(p:point; frombasis:geometry; f:vector):void; // Force `f'isappliedto point`p',expressed inthecoordinate // systemof `frombasis'
applycenterforce(f:vector):void; // finworld coordinates applytorque(t:vector):void;
// t inworld coordinates
applyimpulse(p:point; frombasis:geometry; i:vector):void; // Impulsechange`i'is appliedto point `p',expressed // inthecoordinate systemof `frombasis'
We willcomeback to theexact semanticsof these methodsinSection3.2.3.
Third,it shouldbe possibleto trigger a dyna to advancea frame. This entails making sure that the next motion state is up to date with respect to theequations of motion presented in Section 3.2.1, and the forces, torques and impulse changes that have been applied to it. Then the frameadvancement (e ectivelyincreasing time) is performed byshifting the next motion state into the current motion state. The corresponding method of the dyna class is callednextframe.
Regular (once a frame) triggering, interleaved with calls to the four `apply'methods intro-duced above, mimics the motion under in uence of inertia and the applied forces, torques andimpulsechanges.
Oneotherextensiontotheinterfaceofthedynaclassismade: theintroductionofthenotionof velocitydamping,asa crude (butvery e ective)approximationof frictionwith themedium in which a dyna exists. Methods are introduced to set and retrieve the velocity damping factor (which liesbetween zero and one). This factor is applied to thevelocity(bothlinear and angular) when shifting thenext motion state to the current motion state at theend of each frame (taking thestep sizeinto account), therebyslowingthe dynadown.
The dyna systemclass
Several properties are global to all dynas. The step size h signifying thelength of the time intervalbetween one frame and the next is such a property. To improve the interface with an application designer, a class is introduced to administrate such global properties, and to provide a single entry-point for the dynamics calculations. The dynasystem class can maintainalistofalldynas,andtrigger allofofthem inresponseto asinglecallto amethod \dynamics()"ofitsown. Thisway,anapplicationdesignerdoesnothavetotriggereach dyna individually. After triggering the dynas,the dynasystemcan also callthe callbacks for each of thedynastocopytheirnewlycalculatedmotionstatesto thehostsystem, inthedynamics method (see theendof Section3.2.4).
For the dynasystem to be able to keep track of all dynas, a dyna registers itself with the dynasysteminits constructor, and `checks itself out' in its destructor. Similarly, each non-dynageometryregistersitself withthedynasystem,sothedynasystemcan makesure thatit callsthecallbackswhichensurethatthecompanion'smotionstate isupdatedwhenrequired (see thepseudocode inSection3.3),relievingtheuserofkeepingthecompanionsupto date.
Havingacentraldynasystemclass, enableseasy introductionofotherglobalproperties,such asgravity. Thedynasystemclasswillhaveanattributeforadministratingthisproperty(with methods for setting and retrieving it), and the dynasystem can apply gravity to all dynas (usingtheirapplycenterforcemethod)initsdynamicsmethod,sothatgravitywillalwaysbe applied withoutany userinteraction,oncea user hasspeci edthegravity vector.
Intermezzo
Methodssuch asthemethodsof thedynasystemthat adynacan useto register itself,and a method like the `trigger' method of thedyna classare part of theinterfaceof the class that theybelongto. However,they areonlyforusewithinDynamoitself. Thesearenotmethods thata user of Dynamo should use (or even be aware of). Many more such methods willbe introduced furtheron inthisthesis. One of theshort-comings of C++ (inwhichDynamo is implemented) isthelackofsupportforthis`library',or`package', or`module'notion, which wouldallowaneasywaytospecifywhichmethodsareexportedoutsidethelibrary,andwhich aremeantforinternaluseonly(thefriendnotioninC++istoocumbersomeforthispurpose since it needs explicit naming of the classes within the library, which limits extendibility). This is a speci c example of the problem identi ed in Section 4.3.1 of [Aksit92 ]. In the following text, we will denote this di erence as `exported methods' and `internal methods' (so these are two ` avors' of publicmethods, next to privatemethods). In the user manual of Dynamo (see[Baren 99 ]),onlytheexternalmethodsaredocumentedto deter auser from using the internal ones, but C++ does not provide a mechanism to prevent the user from actually usingany internalmethods.
3.2.3 Integration of the motion equations
In theprevious two sections, the interfaces of the forwarddynamics classes were presented. Tobeable toimplementthefunctionalityspeci edbythese interfaces,themotionequations presentedinSection3.2.1 needto be integrated.
Integrationstep size
The step size h of the simulation is administratedvia the dynasystemclass. This step size can be adjustedbythe user. Furthermore:
Design decision 5: The step size is determined by the user alone: Dynamo does not modifythe stepsize.
Thisdecisionismadebecausethetimestepisglobaltothewholeanimation. IfDynamowere tochangethetimestepituses,thedynasnewmotionstatewouldbecalculatedforadi erent timestamp than the motion states of the geometries that are notunder Dynamo's control. Changingthetimestep,and notifyingthehostsystemofthischange,wouldrequirethehost system to recalculate the motion states for the other geometries also, which it may not be
abletodo(itwouldrequireallother motioncalculationcomponentstosupportbacktracking aframe).
Also,inreal-timecomputeranimation,thestepsizeis oftendirectlyrelatedto thesimulated time between subsequent frames. If Dynamo were to change the step size on its own, this wouldresultinaperceivedslowingdownorspeedingupoftheanimation. Suchtimewarping isnotdesirable.
Since integration errors depend on the step size, this also provides users with a mechanism for error control: users should use a smaller time step if they want more accurate results. Unfortunately, it also burdens the user with the responsibility to provide a step size small enoughforthecalculationstobestable. Theuserisaidedinsolvingthisproblembyproviding methodsthatreportbackforexampletheconstrainterror: these canindicateifthestepsize istoo large.
Time warping can also result from changes in the amount of calculations that need to be performed to calculate the next motion state: if suddenly there is a large increase in the requiredcalculations(whenforexamplemanycollisionsoccur),there willalsobeaperceived slow-down. To try to counter thise ect,thefollowing designdecisionismade:
Design decision 6: Dynamodoesnotsubdividethetimestepspeci edbytheuser'sstep size
Inother words, Dynamo doesnot try to computea consistentmotion state for alldynasfor any moment between t and t+h. Again, such an intermediate motion state might require informationaboutthe non-dynamicallycontrolled geometries, which might notbe available fromthecomponentsthat calculatetheirmotions.
Thisstilldoesnotguarantyaconstant computationtimeperframe,butat leastitprecludes scenarioswherethesimulation comesto agrindinghaltbecause forexamplemanycollisions causethetimeintervalbetweentand t+hto be subdividedinmanysmallerintervals which each cost asmuch timeto integrate asnormallythewholeintervalwould.
Ausercancompensatefortheremainingtimewarpingbysettingthetimestepforeachframe basedonthedi erencebetweenthesimulatedtimeandthetimeprovidedbyareal-timeclock (similarto thewayJavahandles thecorrespondence betweenreal timeand simulatedtime).
Designdecision6posessomelimitationstothewaysDynamocanhandleforexamplecollisions (seeSection 4.1.14)and applyimpulses(seedesign decision10).
Constant forces and torques
To facilitate theintegration of themotion equationspresented inSection3.2.1, equation 3.6 isrewrittento _ ! =AJ 1 A T (M !AJ~ A T !) (3.9) (usingJ w =AJA T
),sothattheconstantJ can beused,and onlythemotionstate variables arenon-constant,next to thetotal force F and total torqueM.
ThedependencyofF andM ontimewithinthetimespanfromoneframetothenextdepends on the semantics of the applyforce, applycenterforce, and applytorque methods. Because of thediscrete modelingof time, we will take these semanticsto be that a constant force or torqueisappliedtothedynabetweenthecurrentframeandthenext. Butthetermconstant is stillopento interpretation,asis shownbelow.
For a non-central force, the applied torque M is speci ed by M = Af, as was seen in Section 3.2.1. This means that the torque exerted by f depends on the orientation of the object. Asthetorqueis being applied,its rotationalacceleration isaddedto already present inertialrotationofthedyna. Sotheorientationwilllikelychangebetweentandt+h,causing in turna change intheamount oftorque thatresultsfrom forcef.
There aretwo ways to interpretthe termconstant whenit comesto force vector f:
1. It can be kept constant in world coordinates. This means that the above mentioned e ect takesplace, and theresultingtorque willvary overthetimeinterval [t::t+h].
2. It can be kept constant in local coordinates. This means that the force will \rotate along" with the dyna. This yields a constant torque between t and t+h, but as the directionoftheforcenowvariesinworldcoordinates,thelineare ect oftheforce isno longerconstant.
Likewise, we can choose to keep the application point of the force constant in local or in world coordinates (but keeping it constant in world coordinates does not make much sense mechanically).
Many forces that are encountered are reaction forces, which come in pairs: a force which is applied to one dyna, and a counter force which is applied to another dyna (satisfying Newton's action= reaction principle). To ensurethat such reaction force pairsremain each other's opposites during the entire integration interval (see the footnote on page 15), we chooseto keep forcesconstant inworld coordinates.
Design decision 7: When applyingforces, theforces arekept constant inworld coordi-nates, and the application point is kept constant in local coordinates, over the course of theintegration interval.
Methods of integration
Design decision7resultsina constantcentralforceF betweentand t+h,whichmeansthat equations 3.3and 3.4can be integratedanalytically:
z t+h = z t +v t h+ 1 2 F m h 2 (3.10) v t+h = v t + F m h (3.11)
For the orientation and angular velocity, numerical integration is required however. For a numerical integrator @: <q t+h ;! t+h >=@(<q t ;! t >;M) (3.12)
Several integrators areavailable, each with their own precision and computationalcost (see AppendixBforsomeexamples). Toallowauser tobalancetheprecisionandcomputational costinhis/her application,we decide:
Design decision 8: We willleave thechoice ofthe motionintegrator open to theuser.
We willtake a lookat thee ects of thisdecisionon thesoftware inthenext section.
Themoreadvancedintegratorsevaluatethedi erentialequations(equations 3.5and3.9)for several timestampsbetweentandt+h. Thisiswherethee ect ofthevaryingtorquecomes in.
We might choose notto account for thise ect, by justcalculating the torquecorresponding to a given force once (using for example the orientation of the current motion state) and then using that torque for each evaluation of the di erential equations. This makes forces that rotate along with the dyna for the duration of the frame indistinguishablefrom forces that do not. In thiscase, the error which is made is of the order of the discretizationstep. Accounting for this e ect at each evaluation of the di erential equations means that the error is of theorder of theintegrator which isused: ifan integrator chooses to evaluatethe di erential equations for di erent time stamps in the integration interval, it will get more accurate results. Inthiswaya usercan getmore accurateresultsbychoosinga higher-order integrator. SuchscalabilitycanbeadvantageouswhenusingDynamoforsimulationpurposes, so thelittleoverhead causedbytheextra administrationistaken forgranted:
Design decision 9: We re-evaluate the torque caused by a force at each evaluation of equation 3.9 to account forthechangeinorientation,and therebygain moreaccuracy.
Beforelookingathowthesedecisionsareincorporatedinthedynaclass,westillhavetode ne the precisesemantics ofthe applyimpulsemethod.
Impulse changes
Intakingatimestepfromtimettotimet+h,impulse(ex)changescouldoccuratanymoment during the interval. However, allowing an impulse change to be applied at any moment in the interval would require sub-sampling the time interval (by maintaining a list of impulse-changes-to-be-applied,anddoingtheintegrationofregularforcesinapiecewisefashion). This would computationally be very expensive, which is unacceptablein the time-criticalmotion integration, as was decided in design decision 6. It would also mean that the user's choice inmotionintegratorwouldgive him/herno guarantyforthecomputationale ortusedinan integration step, which would mean a user would not be able to properly balance precision against computationalspeed. Therefore:
Design decision 10: Impulsechanges areappliedat thestartoftheintegration interval, not atarbitrary pointsintime duringtheinterval.
Theimpulsechangesarealsoassumedtobeinstantaneous,sothatindeedthee ectsofinertia andnon-impulsiveforcescanbeneglectedasthee ectsoftheimpulsechangesarecalculated, as denotedbyequations3.7 and 3.8.
3.2.4 Forward dynamics classes revisited
Now that we know how the motion equations will be integrated, we can revisit the for-ward dynamics classes and see how the notions discovered in the previous section in uence these classes, and can help implementing the methods speci ed in the interfaces given in Section 3.2.2.
The motion state after impulse changes
Due to the possible application of impulse changes at the start of the integration interval, there arenowthree motionstatesofinteresttoadynawhileitismakingthestepfromtimet tot+h. Thetwomotionstatesadministratedbythegeometryclass,andalsothemotionstate at time t after impulsesare applied. The two motion states at time t will be distinguished by t and t
+
. The former is the motion state before application of impulse changes (and is therefore thecommitted result of the motion calculationsof theprevious frame), and the latter is the motion state after application of impulse changes (which may change as more impulsechanges areapplied).
The dyna class is extended with an attribute for the third motion state, and the related methods, so that it can be queried and set just as the two other motion states using the methodsofthegeometryclass. Themethodsinheritedfromthegeometryclassforsettingthe motionstate at timet arere-implementedto setboththestates at t and t
+ .
Forthemotion equations,application ofan impulsechangeresultsin:
v t + = v t +v t ! t + = ! t +! t
and sincepositionsare nota ectedinan in nitelyshorttimeinterval:
z t + = z t =z t A t + = A t =A t q t + = q t =q t
Time integration for theregularforcesis now basedon z t ,v t +,A t ,q t and ! t +.
The motion integrator classes
To implement equation 3.12 we will have to use a numerical integrator to provide function @,and design decision8 speci esthat the choice of the motion integrator is up to theuser. The motion integration is therefore encapsulated in a separate class hierarchy: an abstract
mintegratorclass speci esthe interfacewhicha dyna can useto integrateits motion state, whileinheritanceisusedto implementthisinterfaceindi erentwaysspeci cto eachmotion integrator. Thiswillalso allowforeasy addition ofnew motion integratortypes.
Theintegratemethodisthemostimportantmethodofthemintegratorclass. Itreturnsthe rotationalpartofthemotionstateforthenextframe giventherotationalpartofthemotion state forthisframeand a reference to thedynawhose motionstate is being integrated.
The dyna classo ers supportfor the motionsintegrators throughthe addition of (internal) method\ode". Thismethodimplementstheright-handsidesofequations3.5and3.6yielding for(therotationalpartof)motionstatey: y_=ode(y)(theforcesandtorquesthatarerequired forthis evaluation are attributes of the dyna and do not have to be given as parameters of theodemethod).
The dual representation of the orientation (in both A and q) causes a small complication: after a new value for q has been calculated, A should be updated as well: the integrator shouldcalltheA2qmethodofthenewmotionstate. Toillustratetheuseofthismethod,here aretheintegratemethodsforEulerand forRungeKutta4(moreonthespeci cintegrators andthediscretization ofthemotion equationscan befoundinAppendixB):
euler::integrate(y,d): supvecf var ytemp: supvec;
ytemp:=y+hd.ode(y); ytemp.q2A(); return ytemp;
g
Note that thismethod, despite the syntax used here, will be implemented in C++ (as the methodbelowand mostother code presentedinthisreport).
rungekutta4::integrate(y,d): supvecf var k1, k2, k3, k4, ytemp : supvec;
k1:=d.ode(y ); ytemp:=y+ h 2 k1; ytemp.q2A(); k2:=d.ode(ytemp); ytemp:=y+ h 2 k2; ytemp.q2A();
k3:=d.ode(ytemp); ytemp:=y+hk3; ytemp.q2A(); k4:=d.ode(ytemp); ytemp:=y+ h 6 (k1+2k2+2k3+k4); ytemp.q2A(); return ytemp; g
Theadministrationoftheintegratorishandledbythedynasystemclasswhichwasintroduced earlier: itwillhaveareferencetotheintegratorwhichiscurrentlyused,andthereisamethod throughwhichthedynas can obtainthisreferencewhen they needa motionintegrator.
Dyna's attributes
Motion integration could be performed in each call to applyforce, applytorque, applycenterforce, and applyimpulse. However, since it is likely that more than one force,
torque and/or impulsechange is applied to a dyna before the next motion state is required, weadd some administrationto keeptrack oftheforcesand torqueswhichhave beenapplied so far, to minimizethe number of times that the (computationally expensive) integration is performed. Integration is thenonlyperformed when thenext motionstate isrequired. This optimisation is possiblebecause {under its interface{ thedyna class is solely responsiblefor how itperformsits functions.
Design decision 11: Motionintegrationisperformedon demand,ratherthaneachtime a new forceortorque isapplied.
Two private boolean attributes ShouldIntegrateTrans and ShouldIntegrateRot are intro-duced. The formeradminstrates theinvalidityofthenext motion state's positionand linear velocity. The latter the invalidityof the orientation and angular velocity. This way, the re-sultsof an integration step can be reuseduntilthey areinvalidated byapplication of forces, torques,impulsechanges,ora shiftto thenext frame. Amethod integrateisintroducedin thedyna class. This method is called each time when the next motion state it required,to make sure thatthat next motionstate is upto date.
Thetotalamountofcentralforce(duetocallstoapplycenterforceandapplyforce)andtotal amount of torque (dueto calls to applytorque,butnot to applyforce)are administratedby thedyna: ithasvector-typedattributesFandMto dothis. Forthetorquee ect offorcesthat areappliedthroughtheapplyforcemethod,alist fListofforce/vertexpairswillhaveto be maintainedsince the particular orientation at which the torque is required is not known in advance,sincethisdependson thetimestampforwhichtheintegrator needsto evaluatethe motionequations.
Dyna's methods
The pseudocode belowshows how the considerations presentedin theprevious sectionsare capturedinthedynaclass.
At the very end of each frame the newly calculated motion state is shifted into both cur-rent motion states (accounting for the velocity damping). Both motion states at t be-ing equal signi es that no impulses have yet been applied during the ensuing frame. The ShouldIntegrateTrans and ShouldIntegrateRot attributes are set to true to indicate that therehas beenno integration step performed yetforthenext frame:
dyna::nextframe()f integrate(); <z t ;v t ;q t ;A t ;! t >:=<z t+h ;vdv t+h ;q t+h ;A t+h ;vd! t+h >; <z t +;v t +;q t +;A t +;! t + >:=<z t ;v t ;q t ;A t ;! t >; F:= ~ 0; M:= ~ 0; fList.makeempty(); ShouldIntegrateTrans:=true; ShouldIntegrateRot:=true; g
Thefour methodswhichcan be used to applythe di erentkinds of forcesto a dyna simply modifytheadminstrationoftheforces,torquesand/orimpulseswhichhavebeenappliedthis frame: dyna::applycenterforce(f)f F+=f; ShouldIntegrateTrans:=true; g dyna::applytorque(M)f M+=M; ShouldIntegrateRot:=true; g dyna::applyforce(p,g,f)f var p local :point; p local
:=tolocal(g.toworld(p)); F+=f; fList.add(<p local ,f>); ShouldIntegrateTrans:=true; ShouldIntegrateRot:=true; g dyna::applyimpulse(p,g,)f var p local :point; p local
:=tolocal(g.toworld(p)); v t ++= 1 m ; ! t ++=A t J 1 A T t (A t p local ); ShouldIntegrateTrans:=true; ShouldIntegrateRot:=true; g
As seen in the last method, impulse changes applied during a frame are directly accounted forin the motion state for t
+
. This motion state then functions as a starting point for the algebraicaland numericalintegration to the motionstate at t+h.
Toreallyintegratethemotionvariables,thedynaclasshasthe(internal) integratemethod:
dyna::integrate()f if(ShouldIntegrateTrans)f z t+h :=z t + +v t +h+ 1 2 F m h 2 ; v t+h :=v t ++ F m h; ShouldIntegrateTrans:=false; g if(ShouldIntegrateRot)f <q t+h ;! t+h
>:=dynasystem.getintegrator().integrate( q t ,! t ,this); ShouldIntegrateRot:=false; g g
Forthemotionintegrator to beableto integratetherotationalpart ofthemotion state,the dyna provides theodemethod(only to be calledbythemotion integrators):
