ENC: sixDoFRigidBodyMotion. Adding fixedOrientation constraint.

Adding report Switch to turn on reporting.

Renaming linearSphericalAngularSpring to sphericalAngularSpring and
getting a moment contribution from each axis.  Can supply a reference
orientation.

src/postProcessing/functionObjects/forces/Make/files

@@ -16,7 +16,7 @@ sDoFRBMR = $(sDoFRBM)/sixDoFRigidBodyMotionRestraint
 $(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
 $(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C
 $(sDoFRBMR)/linearSpring/linearSpring.C
-$(sDoFRBMR)/linearSphericalAngularSpring/linearSphericalAngularSpring.C
+$(sDoFRBMR)/sphericalAngularSpring/sphericalAngularSpring.C
 
 sDoFRBMC = $(sDoFRBM)/sixDoFRigidBodyMotionConstraint
 
@@ -25,6 +25,7 @@ $(sDoFRBMC)/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C
 $(sDoFRBMC)/fixedPoint/fixedPoint.C
 $(sDoFRBMC)/fixedPlane/fixedPlane.C
 $(sDoFRBMC)/fixedLine/fixedLine.C
+$(sDoFRBMC)/fixedOrientation/fixedOrientation.C
 
 pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
 pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C

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

@@ -45,9 +45,12 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
 
             restraints_[rI].restrain(*this, rP, rF, rM);
 
-            // Info<< "Restraint " << rI << " force " << rF << nl
+            if (report_)
-            //     << "Restraint " << rI << " moment " << rM
+            {
-            //     << endl;
+                Info<< "Restraint " << restraints_[rI].name() << ": "
+                    << "force " << rF << " moment " << rM
+                    << endl;
+            }
 
             a() += rF/mass_;
 
@@ -121,16 +124,16 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
                 << "iterations (" << maxConstraintIters_ << ") exceeded." << nl
                 << exit(FatalError);
         }
-        else
+        else if (report_)
         {
             Info<< "sixDoFRigidBodyMotion constraints converged in "
                 << iter << " iterations"
-                // << nl << "Constraint force: " << cFA << nl
+                << nl << "Constraint force: " << cFA << nl
-                // << "Constraint moment: " << cMA
+                << "Constraint moment: " << cMA
                 << endl;
         }
 
-        // Add the constrain forces and moments to the motion state variables
+        // Add the constraint forces and moments to the motion state variables
         a() += cFA/mass_;
 
         // The moment of constraint forces has already been added
@@ -150,7 +153,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     maxConstraintIters_(0),
     refCentreOfMass_(vector::zero),
     momentOfInertia_(diagTensor::one*VSMALL),
-    mass_(VSMALL)
+    mass_(VSMALL),
+    report_(false)
 {}
 
 
@@ -164,7 +168,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     const vector& tau,
     scalar mass,
     const point& refCentreOfMass,
-    const diagTensor& momentOfInertia
+    const diagTensor& momentOfInertia,
+    bool report
 )
 :
     motionState_
@@ -181,7 +186,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     maxConstraintIters_(0),
     refCentreOfMass_(refCentreOfMass),
     momentOfInertia_(momentOfInertia),
-    mass_(mass)
+    mass_(mass),
+    report_(report)
 {}
 
 
@@ -193,7 +199,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
     maxConstraintIters_(0),
     refCentreOfMass_(dict.lookupOrDefault("refCentreOfMass", centreOfMass())),
     momentOfInertia_(dict.lookup("momentOfInertia")),
-    mass_(readScalar(dict.lookup("mass")))
+    mass_(readScalar(dict.lookup("mass"))),
+    report_(dict.lookupOrDefault<Switch>("report", false))
 {
     addRestraints(dict);
 
@@ -212,7 +219,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     maxConstraintIters_(sDoFRBM.maxConstraintIters()),
     refCentreOfMass_(sDoFRBM.refCentreOfMass()),
     momentOfInertia_(sDoFRBM.momentOfInertia()),
-    mass_(sDoFRBM.mass())
+    mass_(sDoFRBM.mass()),
+    report_(sDoFRBM.report())
 {}
 
 
@@ -347,6 +355,11 @@ void Foam::sixDoFRigidBodyMotion::updateForce
         v() += 0.5*deltaT*a();
 
         pi() += 0.5*deltaT*tau();
+
+        if(report_)
+        {
+            status();
+        }
     }
 
     Pstream::scatter(motionState_);
@@ -383,8 +396,8 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
 (
     const point& pt,
-    const vector deltaForce,
+    const vector& deltaForce,
-    const vector deltaMoment,
+    const vector& deltaMoment,
     scalar deltaT
 ) const
 {
@@ -402,4 +415,30 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
 }
 
 
+Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
+(
+    const vector& v,
+    const vector& deltaMoment,
+    scalar deltaT
+) const
+{
+    vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
+
+    tensor QTemp = Q();
+
+    rotate(QTemp, piTemp, deltaT);
+
+    return (QTemp & v);
+}
+
+
+void Foam::sixDoFRigidBodyMotion::status() const
+{
+    Info<< "Centre of mass: " << centreOfMass() << nl
+        << "Linear velocity: " << v() << nl
+        << "Angular velocity: " << omega()
+        << endl;
+}
+
+
 // ************************************************************************* //

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

@@ -107,6 +107,9 @@ class sixDoFRigidBodyMotion
         //- Mass of the body
         scalar mass_;
 
+        //- Switch to turn reporting of motion data on and off
+        Switch report_;
+
 
     // Private Member Functions
 
@@ -154,9 +157,6 @@ class sixDoFRigidBodyMotion
             //- Return access to the centre of mass
             inline const point& refCentreOfMass() const;
 
-            //- Return access to the inertia tensor
-            inline const diagTensor& momentOfInertia() const;
-
             //- Return access to the orientation
             inline const tensor& Q() const;
 
@@ -178,9 +178,6 @@ class sixDoFRigidBodyMotion
             //- Return access to the centre of mass
             inline point& refCentreOfMass();
 
-            //- Return non-const access to the inertia tensor
-            inline diagTensor& momentOfInertia();
-
             //- Return non-const access to the orientation
             inline tensor& Q();
 
@@ -215,7 +212,8 @@ public:
             const vector& tau,
             scalar mass,
             const point& refCentreOfMass,
-            const diagTensor& momentOfInertia
+            const diagTensor& momentOfInertia,
+            bool report = false
         );
 
         //- Construct from dictionary
@@ -282,8 +280,18 @@ public:
         point predictedPosition
         (
             const point& pt,
-            const vector deltaForce,
+            const vector& deltaForce,
-            const vector deltaMoment,
+            const vector& deltaMoment,
+            scalar deltaT
+        ) const;
+
+        //- Predict the orientation of the supplied vector after deltaT
+        //  given the current motion state and the additional supplied
+        //  moment
+        vector predictedOrientation
+        (
+            const vector& v,
+            const vector& deltaMoment,
             scalar deltaT
         ) const;
 
@@ -294,20 +302,33 @@ public:
         //  motion state
         inline point currentVelocity(const point& pt) const;
 
+        //- Report the status of the motion
+        void status() const;
+
+
         // Access
 
             //- Return const access to the centre of mass
             inline const point& centreOfMass() const;
 
+            //- Return access to the inertia tensor
+            inline const diagTensor& momentOfInertia() const;
+
             //- Return const access to the mass
             inline scalar mass() const;
 
+            //- Return the report Switch
+            inline bool report() const;
+
 
         // Edit
 
             //- Return non-const access to the centre of mass
             inline point& centreOfMass();
 
+            //- Return non-const access to the inertia tensor
+            inline diagTensor& momentOfInertia();
+
             //- Return non-const access to the mass
             inline scalar& mass();
 

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedOrientation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "sixDoFRigidBodyMotion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFRigidBodyMotionConstraints
+{
+    defineTypeNameAndDebug(fixedOrientation, 0);
+    addToRunTimeSelectionTable
+    (
+        sixDoFRigidBodyMotionConstraint,
+        fixedOrientation,
+        dictionary
+    );
+};
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation
+(
+    const dictionary& sDoFRBMCDict
+)
+:
+    sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
+    fixedOrientation_()
+{
+    read(sDoFRBMCDict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::~fixedOrientation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain
+(
+    const sixDoFRigidBodyMotion& motion,
+    const vector& existingConstraintForce,
+    const vector& existingConstraintMoment,
+    scalar deltaT,
+    vector& constraintPosition,
+    vector& constraintForceIncrement,
+    vector& constraintMomentIncrement
+) const
+{
+    constraintMomentIncrement = vector::zero;
+
+    scalar maxTheta = -SMALL;
+
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+    {
+        vector axis = vector::zero;
+
+        axis[cmpt] = 1;
+
+        vector refDir = vector::zero;
+
+        refDir[(cmpt + 1) % 3] = 1;
+
+        vector predictedDir = motion.predictedOrientation
+        (
+            refDir,
+            existingConstraintMoment,
+            deltaT
+        );
+
+        axis = (fixedOrientation_ & axis);
+
+        refDir = (fixedOrientation_ & refDir);
+
+        // Removing any axis component from predictedDir
+        predictedDir -= (axis & predictedDir)*axis;
+
+        scalar theta = GREAT;
+
+        scalar magPredictedDir = mag(predictedDir);
+
+        if (magPredictedDir > VSMALL)
+        {
+            predictedDir /= magPredictedDir;
+
+            theta = acos(predictedDir & refDir);
+
+            // Recalculating axis to give correct sign to angle
+            axis = (refDir ^ predictedDir);
+
+            scalar magAxis = mag(axis);
+
+            if (magAxis > VSMALL)
+            {
+                axis /= magAxis;
+            }
+            else
+            {
+                axis = vector::zero;
+            }
+        }
+
+        if (theta > maxTheta)
+        {
+            maxTheta = theta;
+        }
+
+        constraintMomentIncrement +=
+           -relaxationFactor_
+           *theta*axis
+            *motion.momentOfInertia()[cmpt]/sqr(deltaT);
+    }
+
+    constraintPosition = motion.centreOfMass();
+
+    constraintForceIncrement = vector::zero;
+
+    return (mag(maxTheta) < tolerance_);
+}
+
+
+bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read
+(
+    const dictionary& sDoFRBMCDict
+)
+{
+    sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
+
+    fixedOrientation_ =
+        sDoFRBMCCoeffs_.lookupOrDefault<tensor>("fixedOrientation", I);
+
+    if (mag(mag(fixedOrientation_) - sqrt(3.0)) > 1e-9)
+    {
+        FatalErrorIn
+        (
+            "Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read"
+            "("
+                "const dictionary& sDoFRBMCDict"
+            ")"
+        )
+            << "fixedOrientation " << fixedOrientation_
+            << " is not a rotation tensor."
+            << exit(FatalError);
+    }
+
+    return true;
+}
+
+// ************************************************************************* //

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

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation
+
+Description
+    sixDoFRigidBodyMotionConstraint.  Orientation of body fixed global
+    space. Only valid where the predicted deviation from alignment is
+    < 90 degrees.
+
+SourceFiles
+    fixedOrientation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedOrientation_H
+#define fixedOrientation_H
+
+#include "sixDoFRigidBodyMotionConstraint.H"
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace sixDoFRigidBodyMotionConstraints
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class fixedOrientation Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedOrientation
+:
+    public sixDoFRigidBodyMotionConstraint
+{
+
+    // Private data
+
+        //- Reference orientation where there is no moment
+        tensor fixedOrientation_;
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedOrientation");
+
+
+    // Constructors
+
+        //- Construct from components
+        fixedOrientation
+        (
+            const dictionary& sDoFRBMCDict
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
+        {
+            return autoPtr<sixDoFRigidBodyMotionConstraint>
+            (
+                new fixedOrientation(*this)
+            );
+        }
+
+
+    // Destructor
+
+        virtual ~fixedOrientation();
+
+
+    // Member Functions
+
+        //- Calculate the constraint position, force and moment.
+        //  Global reference frame vectors.  Returns boolean stating
+        //  whether the constraint been converged to tolerance.
+        virtual bool constrain
+        (
+            const sixDoFRigidBodyMotion& motion,
+            const vector& existingConstraintForce,
+            const vector& existingConstraintMoment,
+            scalar deltaT,
+            vector& constraintPosition,
+            vector& constraintForceIncrement,
+            vector& constraintMomentIncrement
+        ) const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& sDoFRBMCCoeff);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C

@@ -53,7 +53,7 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::fixedPoint
 )
 :
     sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
-    fixedPoint_(sDoFRBMCCoeffs_.lookup("fixedPoint"))
+    fixedPoint_()
 {
     read(sDoFRBMCDict);
 }

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

@@ -126,13 +126,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass() const
 }
 
 
-inline const Foam::diagTensor&
-Foam::sixDoFRigidBodyMotion::momentOfInertia() const
-{
-    return momentOfInertia_;
-}
-
-
 inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q() const
 {
     return motionState_.Q();
@@ -169,12 +162,6 @@ inline Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass()
 }
 
 
-inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia()
-{
-    return momentOfInertia_;
-}
-
-
 inline Foam::tensor& Foam::sixDoFRigidBodyMotion::Q()
 {
     return motionState_.Q();
@@ -253,18 +240,37 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
 }
 
 
+inline const Foam::diagTensor&
+Foam::sixDoFRigidBodyMotion::momentOfInertia() const
+{
+    return momentOfInertia_;
+}
+
+
 inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const
 {
     return mass_;
 }
 
 
+inline bool Foam::sixDoFRigidBodyMotion::report() const
+{
+    return report_;
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
 {
     return motionState_.centreOfMass();
 }
 
 
+inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia()
+{
+    return momentOfInertia_;
+}
+
+
 inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass()
 {
     return mass_;

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C

@@ -39,6 +39,8 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("report")
+        << report_ << token::END_STATEMENT << nl;
 
     if (restraints_.size())
     {

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "linearSphericalAngularSpring.H"
+#include "sphericalAngularSpring.H"
 #include "addToRunTimeSelectionTable.H"
 #include "sixDoFRigidBodyMotion.H"
 
@@ -34,11 +34,11 @@ namespace Foam
 {
 namespace sixDoFRigidBodyMotionRestraints
 {
-    defineTypeNameAndDebug(linearSphericalAngularSpring, 0);
+    defineTypeNameAndDebug(sphericalAngularSpring, 0);
     addToRunTimeSelectionTable
     (
         sixDoFRigidBodyMotionRestraint,
-        linearSphericalAngularSpring,
+        sphericalAngularSpring,
         dictionary
     );
 };
@@ -47,14 +47,14 @@ namespace sixDoFRigidBodyMotionRestraints
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::
+Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
-linearSphericalAngularSpring
+sphericalAngularSpring
 (
     const dictionary& sDoFRBMRDict
 )
 :
     sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
-    refDir_(),
+    refQ_(),
     stiffness_(),
     damping_()
 {
@@ -64,15 +64,15 @@ linearSphericalAngularSpring
 
 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
 
-Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::
+Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
-~linearSphericalAngularSpring()
+~sphericalAngularSpring()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void
-Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::restrain
+Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::restrain
 (
     const sixDoFRigidBodyMotion& motion,
     vector& restraintPosition,
@@ -80,62 +80,62 @@ Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::restrain
     vector& restraintMoment
 ) const
 {
-    vector curDir = motion.currentOrientation(refDir_);
+    restraintMoment = vector::zero;
 
-    vector a = refDir_ ^ curDir;
+    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-
-    scalar magA = mag(a);
-
-    scalar theta = 0;
-
-    if (magA > VSMALL)
     {
-        a /= magA;
+        vector axis = vector::zero;
 
-        theta = acos(curDir & refDir_);
+        axis[cmpt] = 1;
-    }
+
-    else
+        vector refDir = vector::zero;
-    {
+
-        a = vector::zero;
+        refDir[(cmpt + 1) % 3] = 1;
+
+        vector newDir = motion.currentOrientation(refDir);
+
+        axis = (refQ_ & axis);
+
+        refDir = (refQ_ & refDir);
+
+        newDir -= (axis & newDir)*axis;
+
+        restraintMoment += -stiffness_*(refDir ^ newDir);
     }
 
-    restraintMoment = -stiffness_*theta*a - damping_*motion.omega();
+    restraintMoment += -damping_*motion.omega();
 
     restraintForce = vector::zero;
 
     // Not needed to be altered as restraintForce is zero, but set to
-    // centre of mass to be sure of no spurious moment
+    // centreOfMass to be sure of no spurious moment
     restraintPosition = motion.centreOfMass();
 }
 
 
-bool Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring ::read
+bool Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::read
 (
     const dictionary& sDoFRBMRDict
 )
 {
     sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
 
-    sDoFRBMRCoeffs_.lookup("referenceDirection") >> refDir_;
+    refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
 
-    scalar magRefDir(mag(refDir_));
+    if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
-
-    if (magRefDir > VSMALL)
-    {
-        refDir_ /= magRefDir;
-    }
-    else
     {
         FatalErrorIn
         (
-            "bool Foam::sixDoFRigidBodyMotionRestraints::"
+            "Foam::sixDoFRigidBodyMotionConstraints::sphericalAngularSpring::"
-            "linearSphericalAngularSpring ::read"
+            "read"
             "("
                 "const dictionary& sDoFRBMRDict"
             ")"
         )
-            << "referenceDirection has zero length"
+            << "referenceOrientation " << refQ_ << " is not a rotation tensor. "
-            << abort(FatalError);
+            << "mag(referenceOrientation) - sqrt(3) = "
+            << mag(refQ_) - sqrt(3.0) << nl
+            << exit(FatalError);
     }
 
     sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;

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

@@ -23,21 +23,22 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring
+    Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring
 
 Description
-    sixDoFRigidBodyMotionRestraints model.  Linear spherical angular spring.
+    sixDoFRigidBodyMotionRestraints model.  Spherical angular spring.
 
 SourceFiles
-    linearSphericalAngularSpring.C
+    sphericalAngularSpring.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef linearSphericalAngularSpring_H
+#ifndef sphericalAngularSpring_H
-#define linearSphericalAngularSpring_H
+#define sphericalAngularSpring_H
 
 #include "sixDoFRigidBodyMotionRestraint.H"
 #include "point.H"
+#include "tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,17 +49,17 @@ namespace sixDoFRigidBodyMotionRestraints
 {
 
 /*---------------------------------------------------------------------------*\
-                 Class linearSphericalAngularSpring Declaration
+                   Class sphericalAngularSpring Declaration
 \*---------------------------------------------------------------------------*/
 
-class linearSphericalAngularSpring
+class sphericalAngularSpring
 :
     public sixDoFRigidBodyMotionRestraint
 {
     // Private data
 
-        // Reference direction there is no spring reaction torque
+        //- Reference orientation where there is no moment
-        vector refDir_;
+        tensor refQ_;
 
         //- Spring stiffness coefficient (Nm/rad)
         scalar stiffness_;
@@ -70,13 +71,13 @@ class linearSphericalAngularSpring
 public:
 
     //- Runtime type information
-    TypeName("linearSphericalAngularSpring");
+    TypeName("sphericalAngularSpring");
 
 
     // Constructors
 
         //- Construct from components
-        linearSphericalAngularSpring
+        sphericalAngularSpring
         (
             const dictionary& sDoFRBMRDict
         );
@@ -86,14 +87,14 @@ public:
         {
             return autoPtr<sixDoFRigidBodyMotionRestraint>
             (
-                new linearSphericalAngularSpring(*this)
+                new sphericalAngularSpring(*this)
             );
         }
 
 
     // Destructor
 
-        virtual ~linearSphericalAngularSpring();
+        virtual ~sphericalAngularSpring();
 
 
     // Member Functions