diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
index 86f1647d7..f23e00603 100644
--- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
+++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
@@ -76,10 +76,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::A
 {
     // Fix the global force for global rotation constraints
     dimensionedVector fAbs = force();
-    vector& fAbsVal = fAbs.value();
 
-    // Constrain rotation
-    fAbsVal = constrainTranslation(fAbsVal);
+    // Constrain translation
+    constrainTranslation(fAbs.value());
 
     return
     (
@@ -100,9 +99,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot
 {
     // Fix the global moment for global rotation constraints
     dimensionedVector mAbs = moment();
-    vector& mAbsVal = mAbs.value();
 
-    mAbsVal = constrainRotation(mAbsVal);
+    // Constrain rotation
+    constrainRotation(mAbs.value());
 
     return
         inv(momentOfInertia_)
@@ -124,7 +123,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::E
 }
 
 
-Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const
+void Foam::sixDOFqODE::constrainRotation(vector& vec) const
 {
     vector consVec(vector::zero);
 
@@ -165,12 +164,10 @@ Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const
             consVec.z() = 0;
         }
     }
-
-    return consVec;
 }
 
 
-Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const
+void Foam::sixDOFqODE::constrainTranslation(vector& vec) const
 {
     vector consVec(vector::zero);
 
@@ -211,8 +208,6 @@ Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const
             consVec.z() = 0;
         }
     }
-
-    return consVec;
 }
 
 
@@ -263,7 +258,7 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
             false
         )
     ),
-    referentRotation_(vector::zero, 0)
+    referentRotation_(vector(1, 0, 0), 0)
 {
     setCoeffs();
 }
@@ -502,7 +497,7 @@ void Foam::sixDOFqODE::update(const scalar delta)
     Uval.z() = coeffs_[5];
 
     // Constrain velocity
-    Uval = constrainTranslation(Uval);
+    constrainTranslation(Uval);
     coeffs_[3] = Uval.x();
     coeffs_[4] = Uval.y();
     coeffs_[5] = Uval.z();
@@ -515,7 +510,7 @@ void Foam::sixDOFqODE::update(const scalar delta)
     omegaVal.z() = coeffs_[8];
 
     // Constrain omega
-    omegaVal = constrainRotation(omegaVal);
+    constrainRotation(omegaVal);
     coeffs_[6] = omegaVal.x();
     coeffs_[7] = omegaVal.y();
     coeffs_[8] = omegaVal.z();
diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
index a6caf36eb..e61f54f87 100644
--- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
+++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
@@ -197,6 +197,14 @@ class sixDOFqODE
                 const dimensionedVector& omega
             ) const;
 
+            //- Constrain rotation vector in referent or global coordinate
+            //  system
+            void constrainRotation(vector& vec) const;
+
+            //- Constrain translation vector in referent or global coordinate
+            //  system
+            void constrainTranslation(vector& vec) const;
+
 
 public:
 
@@ -292,6 +300,9 @@ public:
             //  coordinate system
             inline void setOmega(const vector& omega);
 
+            //- Set referent coordinate system to apply constraints
+            inline void setReferentRotation(const HamiltonRodriguezRot& rot);
+
 
         // Average motion per time-step