ENH: support all 12 Euler rotation orders (#1292)

- adjust naming of quaternion 'rotationSequence' to be 'eulerOrder'
  to reflect its purpose.

- provide rotation matrices directly for these rotation orders in
  coordinateRotations::euler for case in which the rotation tensor
  is required but not a quaternion.

applications/test/quaternion/Test-quaternion.C

@@ -68,6 +68,21 @@ void printRotation(const quaternion& quat)
 }
 
 
+bool equalTensors(const tensor& rot1, const tensor& rot2)
+{
+    for (direction cmpt=0; cmpt < tensor::nComponents; ++cmpt)
+    {
+        // Cannot be really picky, but SMALL is reasonable
+        if (mag(rot1[cmpt] - rot2[cmpt]) > SMALL)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -108,8 +123,17 @@ int main(int argc, char *argv[])
         "(vector angle)",
         "Rotate about the <vector> by <angle> degrees - eg, '((1 0 0) 45)'"
     );
+    argList::addBoolOption
+    (
+        "verbose",
+        "Additional verbosity"
+    );
+
     argList args(argc, argv);
 
+    const bool verbose = args.found("verbose");
+
+
     vector rotVector;
 
     if (args.readIfPresent("rollPitchYaw", rotVector))
@@ -122,10 +146,18 @@ int main(int argc, char *argv[])
 
         rotVector *= degToRad();
 
-        const quaternion quat(quaternion::rotationSequence::XYZ, rotVector);
+        const quaternion quat(quaternion::eulerOrder::XYZ, rotVector);
 
         printRotation(quat);
+
+        // Euler
+        const tensor rot
+        (
+            euler::rotation(euler::eulerOrder::XYZ, rotVector, false)
+        );
+        printRotation(rot);
     }
+
     if (args.readIfPresent("yawPitchRoll", rotVector))
     {
         Info<< nl
@@ -136,9 +168,16 @@ int main(int argc, char *argv[])
 
         rotVector *= degToRad();
 
-        const quaternion quat(quaternion::rotationSequence::ZYX, rotVector);
+        const quaternion quat(quaternion::eulerOrder::ZYX, rotVector);
 
         printRotation(quat);
+
+        // Euler
+        const tensor rot
+        (
+            euler::rotation(euler::eulerOrder::ZYX, rotVector, false)
+        );
+        printRotation(rot);
     }
     if (args.readIfPresent("euler", rotVector))
     {
@@ -152,7 +191,7 @@ int main(int argc, char *argv[])
 
         rotVector *= degToRad();
 
-        const quaternion quat(quaternion::rotationSequence::ZXZ, rotVector);
+        const quaternion quat(quaternion::eulerOrder::ZXZ, rotVector);
 
         printRotation(quat);
     }
@@ -266,23 +305,44 @@ int main(int argc, char *argv[])
 
     Info<< "Test conversion from and to Euler-angles" << nl;
     vector angles(0.1, 0.2, 0.3);
-    for (int rs=quaternion::ZYX; rs<quaternion::XZX; ++rs)
+
+    for (int rs : quaternion::eulerOrderNames.values())
     {
+        const quaternion::eulerOrder order = quaternion::eulerOrder(rs);
+        const word& orderName = quaternion::eulerOrderNames[order];
+
         if
         (
-            mag
-            (
-                angles
-              - quaternion(quaternion::rotationSequence(rs), angles)
-               .eulerAngles(quaternion::rotationSequence(rs))
-            )
+            mag(angles - quaternion(order, angles).eulerAngles(order))
           > SMALL
         )
         {
             FatalErrorInFunction
-                << "Inconsistent conversion for rotation sequence "
-                << rs << exit(FatalError)
-                << nl;
+                << "Inconsistent conversion for euler rotation order "
+                << orderName << nl << exit(FatalError);
+        }
+
+
+        tensor rotQ(quaternion(order, angles).R());
+        tensor rotE(euler::rotation(order, angles, false));
+
+        if (verbose)
+        {
+            Info<< "euler " << orderName << angles << nl;
+            printRotation(rotE);
+        }
+
+        if (!equalTensors(rotQ, rotE))
+        {
+            WarningInFunction
+                << "Inconsistent quaternion/euler rotation matrices for "
+                << orderName << nl;
+
+            printRotation(rotQ);
+            printRotation(rotE);
+
+            FatalError
+                << nl << exit(FatalError);
         }
     }
 

applications/test/sizeof/Make/options

-/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
 /* EXE_LIBS = -lfiniteVolume */

applications/utilities/mesh/manipulation/transformPoints/Make/options

 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010, 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -57,7 +57,7 @@ Usage
     With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
     it will also read & transform vector & tensor fields.
 
-    Note:
+Note
     roll (rotation about x)
     pitch (rotation about y)
     yaw (rotation about z)
@@ -73,9 +73,11 @@ Usage
 #include "pointFields.H"
 #include "transformField.H"
 #include "transformGeometricField.H"
-#include "unitConversion.H"
+#include "axisAngleRotation.H"
+#include "EulerCoordinateRotation.H"
 
 using namespace Foam;
+using namespace Foam::coordinateRotations;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -284,40 +286,38 @@ int main(int argc, char *argv[])
         n1n2[0].normalise();
         n1n2[1].normalise();
 
-        const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
+        const tensor rot(rotationTensor(n1n2[0], n1n2[1]));
 
-        Info<< "Rotating points by " << rotT << endl;
-
-        points = transform(rotT, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, rotT);
+            rotateFields(args, runTime, rot);
         }
     }
     else if (args.found("rotate-angle"))
     {
-        const Tuple2<vector, scalar> axisAngle
+        const Tuple2<vector, scalar> rotAxisAngle
         (
             args.lookup("rotate-angle")()
         );
 
+        const vector& axis  = rotAxisAngle.first();
+        const scalar& angle = rotAxisAngle.second();
+
         Info<< "Rotating points " << nl
-            << "    about " << axisAngle.first() << nl
-            << "    angle " << axisAngle.second() << nl;
+            << "    about " << axis << nl
+            << "    angle " << angle << nl;
 
-        const quaternion quat
-        (
-            axisAngle.first(),
-            degToRad(axisAngle.second())
-        );
+        const tensor rot(axisAngle::rotation(axis, angle, true));
 
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, quat.R());
+            rotateFields(args, runTime, rot);
         }
     }
     else if (args.readIfPresent("rollPitchYaw", v))
@@ -327,16 +327,14 @@ int main(int argc, char *argv[])
             << "    pitch " << v.y() << nl
             << "    yaw   " << v.z() << nl;
 
-        v *= degToRad();
-
-        const quaternion quat(quaternion::rotationSequence::XYZ, v);
+        const tensor rot(euler::rotation(euler::eulerOrder::XYZ, v, true));
 
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, quat.R());
+            rotateFields(args, runTime, rot);
         }
     }
     else if (args.readIfPresent("yawPitchRoll", v))
@@ -346,16 +344,14 @@ int main(int argc, char *argv[])
             << "    pitch " << v.y() << nl
             << "    roll  " << v.z() << nl;
 
-        v *= degToRad();
-
-        const quaternion quat(quaternion::rotationSequence::ZYX, v);
+        const tensor rot(euler::rotation(euler::eulerOrder::ZYX, v, true));
 
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, quat.R());
+            rotateFields(args, runTime, rot);
         }
     }
 

applications/utilities/surface/surfaceTransformPoints/Make/options

 EXE_INC = \
-    -I$(LIB_SRC)/surfMesh/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lsurfMesh
+    -lsurfMesh \
+    -lmeshTools

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -33,12 +33,18 @@ Description
     Transform (scale/rotate) a surface.
     Like transformPoints but for surfaces.
 
-    The rollPitchYaw option takes three angles (degrees):
-    - roll (rotation about x) followed by
-    - pitch (rotation about y) followed by
-    - yaw (rotation about z)
+    The rollPitchYaw and yawPitchRoll options take three angles (degrees)
+    that describe the intrinsic Euler rotation.
 
-    The yawPitchRoll does yaw followed by pitch followed by roll.
+    rollPitchYaw
+    - roll (rotation about X) followed by
+    - pitch (rotation about Y) followed by
+    - yaw (rotation about Z)
+
+    yawPitchRoll
+    - yaw (rotation about Z) followed by
+    - pitch (rotation about Y) followed by
+    - roll (rotation about X)
 
 \*---------------------------------------------------------------------------*/
 
@@ -48,12 +54,12 @@ Description
 #include "transformField.H"
 #include "Pair.H"
 #include "Tuple2.H"
-#include "quaternion.H"
-#include "unitConversion.H"
-
+#include "axisAngleRotation.H"
+#include "EulerCoordinateRotation.H"
 #include "MeshedSurfaces.H"
 
 using namespace Foam;
+using namespace Foam::coordinateRotations;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -179,31 +185,29 @@ int main(int argc, char *argv[])
         n1n2[0].normalise();
         n1n2[1].normalise();
 
-        const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
+        const tensor rot(rotationTensor(n1n2[0], n1n2[1]));
 
-        Info<< "Rotating points by " << rotT << endl;
-
-        points = transform(rotT, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
     }
     else if (args.found("rotate-angle"))
     {
-        const Tuple2<vector, scalar> axisAngle
+        const Tuple2<vector, scalar> rotAxisAngle
         (
             args.lookup("rotate-angle")()
         );
 
+        const vector& axis  = rotAxisAngle.first();
+        const scalar& angle = rotAxisAngle.second();
+
         Info<< "Rotating points " << nl
-            << "    about " << axisAngle.first() << nl
-            << "    angle " << axisAngle.second() << nl;
+            << "    about " << axis << nl
+            << "    angle " << angle << nl;
 
-        const quaternion quat
-        (
-            axisAngle.first(),
-            degToRad(axisAngle.second())
-        );
+        const tensor rot(axisAngle::rotation(axis, angle, true));
 
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
     }
     else if (args.readIfPresent("rollPitchYaw", v))
     {
@@ -212,12 +216,10 @@ int main(int argc, char *argv[])
             << "    pitch " << v.y() << nl
             << "    yaw   " << v.z() << nl;
 
-        v *= degToRad();
+        const tensor rot(euler::rotation(euler::eulerOrder::XYZ, v, true));
 
-        const quaternion quat(quaternion::rotationSequence::XYZ, v);
-
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
     }
     else if (args.readIfPresent("yawPitchRoll", v))
     {
@@ -226,12 +228,10 @@ int main(int argc, char *argv[])
             << "    pitch " << v.y() << nl
             << "    roll  " << v.z() << nl;
 
-        v *= degToRad();
-
-        const quaternion quat(quaternion::rotationSequence::ZYX, v);
+        const tensor rot(euler::rotation(euler::eulerOrder::ZYX, v, true));
 
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
+        Info<< "Rotating points by " << rot << endl;
+        points = transform(rot, points);
     }
 
     List<scalar> scaling;

src/OpenFOAM/primitives/quaternion/quaternion.C

@@ -31,10 +31,30 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const char* const Foam::quaternion::typeName = "quaternion";
 const Foam::quaternion Foam::quaternion::zero(0, vector(0, 0, 0));
 const Foam::quaternion Foam::quaternion::I(1, vector(0, 0, 0));
 
+const Foam::Enum<Foam::quaternion::eulerOrder>
+Foam::quaternion::eulerOrderNames
+({
+    // Proper Euler angles
+    { eulerOrder::XZX, "xzx" },
+    { eulerOrder::XYX, "xyx" },
+    { eulerOrder::YXY, "yxy" },
+    { eulerOrder::YZY, "yzy" },
+    { eulerOrder::ZYZ, "zyz" },
+    { eulerOrder::ZXZ, "zxz" },
+
+    // Tait-Bryan angles
+    { eulerOrder::XZY, "xzy" },
+    { eulerOrder::XYZ, "xyz" },
+    { eulerOrder::YXZ, "yxz" },
+    { eulerOrder::YZX, "yzx" },
+    { eulerOrder::ZYX, "zyx" },
+    { eulerOrder::ZXY, "zxy" },
+});
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::quaternion::quaternion(Istream& is)

src/OpenFOAM/primitives/quaternion/quaternion.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "tensor.H"
 #include "word.H"
 #include "contiguous.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -90,14 +91,23 @@ class quaternion
 
 public:
 
-    //- Component type
-    typedef scalar cmptType;
+    // Public Types
 
-    //- Euler-angle rotation sequence
-    enum rotationSequence
-    {
-        ZYX, ZYZ, ZXY, ZXZ, YXZ, YXY, YZX, YZY, XYZ, XYX, XZY, XZX
-    };
+        //- Component type
+        typedef scalar cmptType;
+
+        //- Euler-angle rotation order
+        enum eulerOrder : unsigned char
+        {
+            // Proper Euler angles
+            XZX, XYX, YXY, YZY, ZYZ, ZXZ,
+
+            // Tait-Bryan angles
+            XZY, XYZ, YXZ, YZX, ZYX, ZXY
+        };
+
+        //- The names for Euler-angle rotation order
+        static const Enum<eulerOrder> eulerOrderNames;
 
 
     // Member Constants
@@ -108,7 +118,7 @@ public:
 
     // Static Data Members
 
-        static const char* const typeName;
+        static constexpr const char* const typeName = "quaternion";
 
         static const quaternion zero;
         static const quaternion I;
@@ -119,15 +129,17 @@ public:
         //- Construct null
         inline quaternion();
 
+        //- Construct zero initialized
+        inline quaternion(const Foam::zero);
+
         //- Construct given scalar and vector parts
         inline quaternion(const scalar w, const vector& v);
 
-        //- Construct a rotation quaternion given the direction d
-        //  and angle theta
+        //- Construct rotation quaternion given direction d and angle theta
         inline quaternion(const vector& d, const scalar theta);
 
-        //- Construct a rotation quaternion given the direction d
-        //- and cosine angle cosTheta and a if d is normalized
+        //- Construct a rotation quaternion given direction d
+        //- and cosine angle cosTheta and flag if d is normalized
         inline quaternion
         (
             const vector& d,
@@ -135,7 +147,8 @@ public:
             const bool normalized
         );
 
-        //- Construct a real from the given scalar part, the vector part = zero
+        //- Construct a real quaternion from the given scalar part,
+        //- the vector part = zero
         inline explicit quaternion(const scalar w);
 
         //- Construct a pure imaginary quaternion given the vector part,
@@ -146,8 +159,8 @@ public:
         //- (w = sqrt(1 - |sqr(v)|))
         static inline quaternion unit(const vector& v);
 
-        //- Construct from three Euler angles
-        inline quaternion(const rotationSequence rs, const vector& angles);
+        //- Construct from three Euler rotation angles
+        inline quaternion(const eulerOrder order, const vector& angles);
 
         //- Construct from a rotation tensor
         inline explicit quaternion(const tensor& rotationTensor);
@@ -169,9 +182,9 @@ public:
        //- The rotation tensor corresponding the quaternion
        inline tensor R() const;
 
-       //- Return a vector of euler angles corresponding to the
-       //- specified rotation sequence
-       inline vector eulerAngles(const rotationSequence rs) const;
+       //- Return the Euler rotation angles corresponding to the
+       //- specified rotation order
+       inline vector eulerAngles(const eulerOrder order) const;
 
        //- Return the quaternion normalized by its magnitude
        inline quaternion normalized() const;
@@ -238,7 +251,7 @@ inline quaternion normalize(const quaternion& q);
 inline quaternion inv(const quaternion& q);
 
 //- Return a string representation of a quaternion
-word name(const quaternion&);
+word name(const quaternion& q);
 
 //- Spherical linear interpolation of quaternions
 quaternion slerp

src/OpenFOAM/primitives/quaternion/quaternionI.H

@@ -31,6 +31,13 @@ inline Foam::quaternion::quaternion()
 {}
 
 
+inline Foam::quaternion::quaternion(const Foam::zero)
+:
+    w_(Zero),
+    v_(Zero)
+{}
+
+
 inline Foam::quaternion::quaternion(const scalar w, const vector& v)
 :
     w_(w),
@@ -88,87 +95,112 @@ inline Foam::quaternion Foam::quaternion::unit(const vector& v)
 
 inline Foam::quaternion::quaternion
 (
-    const rotationSequence rs,
+    const eulerOrder order,
     const vector& angles
 )
 {
-    switch(rs)
+    switch (order)
     {
         case ZYX:
+        {
             operator=(quaternion(vector(0, 0, 1), angles.x()));
             operator*=(quaternion(vector(0, 1, 0), angles.y()));
             operator*=(quaternion(vector(1, 0, 0), angles.z()));
             break;
+        }
 
         case ZYZ:
+        {
             operator=(quaternion(vector(0, 0, 1), angles.x()));
             operator*=(quaternion(vector(0, 1, 0), angles.y()));
             operator*=(quaternion(vector(0, 0, 1), angles.z()));
             break;
+        }
 
         case ZXY:
+        {
             operator=(quaternion(vector(0, 0, 1), angles.x()));
             operator*=(quaternion(vector(1, 0, 0), angles.y()));
             operator*=(quaternion(vector(0, 1, 0), angles.z()));
             break;
+        }
 
         case ZXZ:
+        {
             operator=(quaternion(vector(0, 0, 1), angles.x()));
             operator*=(quaternion(vector(1, 0, 0), angles.y()));
             operator*=(quaternion(vector(0, 0, 1), angles.z()));
             break;
+        }
 
         case YXZ:
+        {
             operator=(quaternion(vector(0, 1, 0), angles.x()));
             operator*=(quaternion(vector(1, 0, 0), angles.y()));
             operator*=(quaternion(vector(0, 0, 1), angles.z()));
             break;
+        }
 
         case YXY:
+        {
             operator=(quaternion(vector(0, 1, 0), angles.x()));
             operator*=(quaternion(vector(1, 0, 0), angles.y()));
             operator*=(quaternion(vector(0, 1, 0), angles.z()));
             break;