src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Class
  26     Foam::sixDoFRigidBodyMotion
  27
  28 Description
  29     Six degree of freedom motion for a rigid body.  Angular momentum
  30     stored in body fixed reference frame.  Reference orientation of
  31     the body (where Q = I) must align with the cartesian axes such
  32     that the Inertia tensor is in principle component form.
  33
  34     Symplectic motion as per:
  35
  36     title = {Symplectic splitting methods for rigid body molecular dynamics},
  37     publisher = {AIP},
  38     year = {1997},
  39     journal = {The Journal of Chemical Physics},
  40     volume = {107},
  41     number = {15},
  42     pages = {5840-5851},
  43     url = {http://link.aip.org/link/?JCP/107/5840/1},
  44     doi = {10.1063/1.474310}
  45
  46     Can add restraints (i.e. a spring) and constraints (i.e. motion
  47     may only be on a plane).
  48
  49 SourceFiles
  50     sixDoFRigidBodyMotionI.H
  51     sixDoFRigidBodyMotion.C
  52     sixDoFRigidBodyMotionIO.C
  53
  54 \*---------------------------------------------------------------------------*/
  55
  56 #ifndef sixDoFRigidBodyMotion_H
  57 #define sixDoFRigidBodyMotion_H
  58
  59 #include "sixDoFRigidBodyMotionState.H"
  60 #include "pointField.H"
  61 #include "sixDoFRigidBodyMotionRestraint.H"
  62 #include "sixDoFRigidBodyMotionConstraint.H"
  63
  64
  65 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  66
  67 namespace Foam
  68 {
  69
  70 // Forward declaration of classes
  71 class Istream;
  72 class Ostream;
  73
  74 // Forward declaration of friend functions and operators
  75 class sixDoFRigidBodyMotion;
  76 Istream& operator>>(Istream&, sixDoFRigidBodyMotion&);
  77 Ostream& operator<<(Ostream&, const sixDoFRigidBodyMotion&);
  78
  79
  80 /*---------------------------------------------------------------------------*\
  81                       Class sixDoFRigidBodyMotion Declaration
  82 \*---------------------------------------------------------------------------*/
  83
  84 class sixDoFRigidBodyMotion
  85 {
  86     // Private data
  87
  88         //- Motion state data object
  89         sixDoFRigidBodyMotionState motionState_;
  90
  91         //- Restraints on the motion
  92         PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
  93
  94         //- Names of the restraints
  95         wordList restraintNames_;
  96
  97         //- Constaints on the motion
  98         PtrList<sixDoFRigidBodyMotionConstraint> constraints_;
  99
 100         //- Names of the constraints
 101         wordList constraintNames_;
 102
 103         //- Maximum number of iterations allowed to attempt to obey
 104         //  constraints
 105         label maxConstraintIterations_;
 106
 107         //- Centre of mass of initial state
 108         point initialCentreOfMass_;
 109
 110         //- Orientation of initial state
 111         tensor initialQ_;
 112
 113         //- Moment of inertia of the body in reference configuration
 114         //  (Q = I)
 115         diagTensor momentOfInertia_;
 116
 117         //- Mass of the body
 118         scalar mass_;
 119
 120         //- Switch to turn reporting of motion data on and off
 121         Switch report_;
 122
 123
 124     // Private Member Functions
 125
 126         //- Calculate the rotation tensor around the body reference
 127         //  frame x-axis by the given angle
 128         inline tensor rotationTensorX(scalar deltaT) const;
 129
 130         //- Calculate the rotation tensor around the body reference
 131         //  frame y-axis by the given angle
 132         inline tensor rotationTensorY(scalar deltaT) const;
 133
 134         //- Calculate the rotation tensor around the body reference
 135         //  frame z-axis by the given angle
 136         inline tensor rotationTensorZ(scalar deltaT) const;
 137
 138         //- Apply rotation tensors to Q for the given torque (pi) and deltaT
 139         inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
 140
 141         //- Apply the restraints to the object
 142         void applyRestraints();
 143
 144         //- Apply the constraints to the object
 145         void applyConstraints(scalar deltaT);
 146
 147         // Access functions retained as private because of the risk of
 148         // confusion over what is a body local frame vector and what is global
 149
 150         // Access
 151
 152             //- Return access to the motion state
 153             inline const sixDoFRigidBodyMotionState& motionState() const;
 154
 155             //- Return access to the restraints
 156             inline const PtrList<sixDoFRigidBodyMotionRestraint>&
 157                 restraints() const;
 158
 159             //- Return access to the restraintNames
 160             inline const wordList& restraintNames() const;
 161
 162             //- Return access to the constraints
 163             inline const PtrList<sixDoFRigidBodyMotionConstraint>&
 164                 constraints() const;
 165
 166             //- Return access to the constraintNames
 167             inline const wordList& constraintNames() const;
 168
 169             //- Return access to the maximum allowed number of
 170             //  constraint iterations
 171             inline label maxConstraintIterations() const;
 172
 173             //- Return access to the initial centre of mass
 174             inline const point& initialCentreOfMass() const;
 175
 176             //- Return access to the initial orientation
 177             inline const tensor& initialQ() const;
 178
 179             //- Return access to the orientation
 180             inline const tensor& Q() const;
 181
 182             //- Return access to velocity
 183             inline const vector& v() const;
 184
 185             //- Return access to acceleration
 186             inline const vector& a() const;
 187
 188             //- Return access to angular momentum
 189             inline const vector& pi() const;
 190
 191             //- Return access to torque
 192             inline const vector& tau() const;
 193
 194
 195         // Edit
 196
 197             //- Return access to the centre of mass
 198             inline point& initialCentreOfMass();
 199
 200             //- Return access to the centre of mass
 201             inline tensor& initialQ();
 202
 203             //- Return non-const access to the orientation
 204             inline tensor& Q();
 205
 206             //- Return non-const access to vector
 207             inline vector& v();
 208
 209             //- Return non-const access to acceleration
 210             inline vector& a();
 211
 212             //- Return non-const access to angular momentum
 213             inline vector& pi();
 214
 215             //- Return non-const access to torque
 216             inline vector& tau();
 217
 218
 219 public:
 220
 221     // Constructors
 222
 223         //- Construct null
 224         sixDoFRigidBodyMotion();
 225
 226         //- Construct from components
 227         sixDoFRigidBodyMotion
 228         (
 229             const point& centreOfMass,
 230             const tensor& Q,
 231             const vector& v,
 232             const vector& a,
 233             const vector& pi,
 234             const vector& tau,
 235             scalar mass,
 236             const point& initialCentreOfMass,
 237             const tensor& initialQ,
 238             const diagTensor& momentOfInertia,
 239             bool report = false
 240         );
 241
 242         //- Construct from dictionary
 243         sixDoFRigidBodyMotion(const dictionary& dict);
 244
 245         //- Construct as copy
 246         sixDoFRigidBodyMotion(const sixDoFRigidBodyMotion&);
 247
 248
 249     //- Destructor
 250     ~sixDoFRigidBodyMotion();
 251
 252
 253     // Member Functions
 254
 255         //- Add restraints to the motion, public to allow external
 256         //  addition of restraints after construction
 257         void addRestraints(const dictionary& dict);
 258
 259         //- Add restraints to the motion, public to allow external
 260         //  addition of restraints after construction
 261         void addConstraints(const dictionary& dict);
 262
 263         //- First leapfrog velocity adjust and motion part, required
 264         // before force calculation
 265         void updatePosition
 266         (
 267             scalar deltaT
 268         );
 269
 270         //- Second leapfrog velocity adjust part, required after motion and
 271         // force calculation
 272         void updateForce
 273         (
 274             const vector& fGlobal,
 275             const vector& tauGlobal,
 276             scalar deltaT
 277         );
 278
 279         //- Global forces supplied at locations, calculating net force
 280         //  and moment
 281         void updateForce
 282         (
 283             const pointField& positions,
 284             const vectorField& forces,
 285             scalar deltaT
 286         );
 287
 288         //- Transform the given initial state pointField by the current
 289         //  motion state
 290         inline tmp<pointField> currentPosition
 291         (
 292             const pointField& pInitial
 293         ) const;
 294
 295         //- Transform the given initial state point by the current motion
 296         //  state
 297         inline point currentPosition(const point& pInitial) const;
 298
 299         //- Transform the given initial state direction by the current
 300         //  motion state
 301         inline vector currentOrientation(const vector& vInitial) const;
 302
 303         //- Access the orientation tensor, Q.
 304         //  globalVector = Q & bodyLocalVector
 305         //  bodyLocalVector = Q.T() & globalVector
 306         inline const tensor& orientation() const;
 307
 308         //- Predict the position of the supplied initial state point
 309         //  after deltaT given the current motion state and the
 310         //  additional supplied force and moment
 311         point predictedPosition
 312         (
 313             const point& pInitial,
 314             const vector& deltaForce,
 315             const vector& deltaMoment,
 316             scalar deltaT
 317         ) const;
 318
 319         //- Predict the orientation of the supplied initial state
 320         //  vector after deltaT given the current motion state and the
 321         //  additional supplied moment
 322         vector predictedOrientation
 323         (
 324             const vector& vInitial,
 325             const vector& deltaMoment,
 326             scalar deltaT
 327         ) const;
 328
 329         //- Return the angular velocity in the global frame
 330         inline vector omega() const;
 331
 332         //- Return the velocity of a position given by the current
 333         //  motion state
 334         inline point currentVelocity(const point& pt) const;
 335
 336         //- Report the status of the motion
 337         void status() const;
 338
 339
 340         // Access
 341
 342             //- Return const access to the centre of mass
 343             inline const point& centreOfMass() const;
 344
 345             //- Return access to the inertia tensor
 346             inline const diagTensor& momentOfInertia() const;
 347
 348             //- Return const access to the mass
 349             inline scalar mass() const;
 350
 351             //- Return the report Switch
 352             inline bool report() const;
 353
 354
 355         // Edit
 356
 357             //- Return non-const access to the centre of mass
 358             inline point& centreOfMass();
 359
 360             //- Return non-const access to the inertia tensor
 361             inline diagTensor& momentOfInertia();
 362
 363             //- Return non-const access to the mass
 364             inline scalar& mass();
 365
 366
 367         //- Write
 368         void write(Ostream&) const;
 369
 370
 371     // IOstream Operators
 372
 373         friend Istream& operator>>(Istream&, sixDoFRigidBodyMotion&);
 374         friend Ostream& operator<<(Ostream&, const sixDoFRigidBodyMotion&);
 375 };
 376
 377
 378 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 379
 380 } // End namespace Foam
 381
 382 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 383
 384 #include "sixDoFRigidBodyMotionI.H"
 385
 386 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 387
 388 #endif
 389
 390 // ************************************************************************* //