Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

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
+            mag(angles - quaternion(order, angles).eulerAngles(order))
-            (
-                angles
-              - quaternion(quaternion::rotationSequence(rs), angles)
-               .eulerAngles(quaternion::rotationSequence(rs))
-            )
           > SMALL
         )
         {
             FatalErrorInFunction
-                << "Inconsistent conversion for rotation sequence "
+                << "Inconsistent conversion for euler rotation order "
-                << rs << exit(FatalError)
+                << orderName << nl << exit(FatalError);
-                << nl;
+        }
+
+
+        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

@@ -1,2 +1,5 @@
-/* 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

@@ -1,6 +1,8 @@
 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) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 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;
+        Info<< "Rotating points by " << rot << endl;
-
+        points = transform(rot, points);
-        points = transform(rotT, 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
+            << "    about " << axis << nl
-            << "    angle " << axisAngle.second() << nl;
+            << "    angle " << angle << nl;
 
-        const quaternion quat
+        const tensor rot(axisAngle::rotation(axis, angle, true));
-        (
-            axisAngle.first(),
-            degToRad(axisAngle.second())
-        );
 
-        Info<< "Rotating points by quaternion " << quat << endl;
+        Info<< "Rotating points by " << rot << endl;
-        points = transform(quat, points);
+        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 tensor rot(euler::rotation(euler::eulerOrder::XYZ, v, true));
 
-        const quaternion quat(quaternion::rotationSequence::XYZ, v);
+        Info<< "Rotating points by " << rot << endl;
-
+        points = transform(rot, points);
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, 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 tensor rot(euler::rotation(euler::eulerOrder::ZYX, v, true));
 
-        const quaternion quat(quaternion::rotationSequence::ZYX, v);
+        Info<< "Rotating points by " << rot << endl;
-
+        points = transform(rot, points);
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, points);
 
         if (doRotateFields)
         {
-            rotateFields(args, runTime, quat.R());
+            rotateFields(args, runTime, rot);
         }
     }
 

applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    surfaceMeshTriangulate
+    surfaceMeshExtract
 
 Group
     grpSurfaceUtilities
@@ -47,6 +47,7 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "polyMesh.H"
+#include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "ListListOps.H"
 #include "uindirectPrimitivePatch.H"
@@ -58,6 +59,57 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+labelList getSelectedPatches
+(
+    const polyBoundaryMesh& patches,
+    const wordRes& whitelist,
+    const wordRes& blacklist
+)
+{
+    DynamicList<label> patchIDs(patches.size());
+
+    for (const polyPatch& pp : patches)
+    {
+        if (isType<emptyPolyPatch>(pp))
+        {
+            continue;
+        }
+        else if (Pstream::parRun() && bool(isA<processorPolyPatch>(pp)))
+        {
+            break; // No processor patches for parallel output
+        }
+
+        const word& patchName = pp.name();
+
+        bool accept = false;
+
+        if (whitelist.size())
+        {
+            const auto matched = whitelist.matched(patchName);
+
+            accept =
+            (
+                matched == wordRe::LITERAL
+              ? true
+              : (matched == wordRe::REGEX && !blacklist.match(patchName))
+            );
+        }
+        else
+        {
+            accept = !blacklist.match(patchName);
+        }
+
+        if (accept)
+        {
+            patchIDs.append(pp.index());
+        }
+    }
+
+    return patchIDs.shrink();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
@@ -69,6 +121,10 @@ int main(int argc, char *argv[])
     );
     timeSelector::addOptions();
 
+    // Less frequently used - reduce some clutter
+    argList::setAdvanced("decomposeParDict");
+    argList::setAdvanced("noFunctionObjects");
+
     argList::addArgument("output", "The output surface file");
 
     #include "addRegionOption.H"
@@ -78,18 +134,26 @@ int main(int argc, char *argv[])
         "Exclude processor patches"
     );
     argList::addOption
+    (
+        "faceZones",
+        "wordRes",
+        "Specify single or multiple faceZones to extract\n"
+        "Eg, 'cells' or '( slice \"mfp-.*\" )'"
+    );
+    argList::addOption
     (
         "patches",
-        "wordRes"
+        "wordRes",
         "Specify single patch or multiple patches to extract.\n"
         "Eg, 'top' or '( front \".*back\" )'"
     );
     argList::addOption
     (
-        "faceZones",
+        "excludePatches",
         "wordRes",
-        "Specify single or multiple faceZones to extract\n"
+        "Specify single patch or multiple patches to exclude from writing."
-        "Eg, 'cells' or '( slice \"mfp-.*\" )'"
+        " Eg, 'outlet' or '( inlet \".*Wall\" )'",
+        true  // mark as an advanced option
     );
 
     #include "setRootCase.H"
@@ -121,6 +185,25 @@ int main(int argc, char *argv[])
         Info<< "Excluding all processor patches." << nl << endl;
     }
 
+    wordRes includePatches, excludePatches;
+    if (args.readListIfPresent<wordRe>("patches", includePatches))
+    {
+        Info<< "Including patches " << flatOutput(includePatches)
+            << nl << endl;
+    }
+    if (args.readListIfPresent<wordRe>("excludePatches", excludePatches))
+    {
+        Info<< "Excluding patches " << flatOutput(excludePatches)
+            << nl << endl;
+    }
+
+    // Non-mandatory
+    const wordRes selectedFaceZones(args.getList<wordRe>("faceZones", false));
+    if (selectedFaceZones.size())
+    {
+        Info<< "Including faceZones " << flatOutput(selectedFaceZones)
+            << nl << endl;
+    }
 
     Info<< "Reading mesh from time " << runTime.value() << endl;
 
@@ -171,53 +254,25 @@ int main(int argc, char *argv[])
         // Construct table of patches to include.
         const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
 
-        labelList includePatches;
+        labelList patchIds =
+        (
+            (includePatches.size() || excludePatches.size())
+          ? getSelectedPatches(bMesh, includePatches, excludePatches)
+          : includeProcPatches
+          ? identity(bMesh.size())
+          : identity(bMesh.nNonProcessor())
+        );
 
-        if (args.found("patches"))
+        labelList faceZoneIds;
-        {
-            includePatches =
-                bMesh.patchSet(args.getList<wordRe>("patches")).sortedToc();
-        }
-        else if (includeProcPatches)
-        {
-            includePatches = identity(bMesh.size());
-        }
-        else
-        {
-            includePatches = identity(bMesh.nNonProcessor());
-        }
-
-
-        labelList includeFaceZones;
 
         const faceZoneMesh& fzm = mesh.faceZones();
 
-        if (args.found("faceZones"))
+        if (selectedFaceZones.size())
         {
-            const wordList allZoneNames(fzm.names());
+            faceZoneIds = fzm.indices(selectedFaceZones);
-
-            const wordRes zoneNames(args.getList<wordRe>("faceZones"));
-
-            labelHashSet hashed(2*fzm.size());
-
-            for (const wordRe& zoneName : zoneNames)
-            {
-                labelList zoneIDs = findStrings(zoneName, allZoneNames);
-                hashed.insert(zoneIDs);
-
-                if (zoneIDs.empty())
-                {
-                    WarningInFunction
-                        << "Cannot find any faceZone name matching "
-                        << zoneName << endl;
-                }
-            }
-
-            includeFaceZones = hashed.sortedToc();
 
             Info<< "Additionally extracting faceZones "
-                << UIndirectList<word>(allZoneNames, includeFaceZones)
+                << fzm.names(selectedFaceZones) << nl;
-                << endl;
         }
 
 
@@ -232,7 +287,7 @@ int main(int argc, char *argv[])
             //  processor patches)
             HashTable<label> patchSize(1024);
             label nFaces = 0;
-            for (const label patchi : includePatches)
+            for (const label patchi : patchIds)
             {
                 const polyPatch& pp = bMesh[patchi];
                 patchSize.insert(pp.name(), pp.size());
@@ -240,7 +295,7 @@ int main(int argc, char *argv[])
             }
 
             HashTable<label> zoneSize(1024);
-            for (const label zonei : includeFaceZones)
+            for (const label zonei : faceZoneIds)
             {
                 const faceZone& pp = fzm[zonei];
                 zoneSize.insert(pp.name(), pp.size());
@@ -289,7 +344,7 @@ int main(int argc, char *argv[])
             compactZones.setCapacity(nFaces);
 
             // Collect faces on patches
-            for (const label patchi : includePatches)
+            for (const label patchi : patchIds)
             {
                 const polyPatch& pp = bMesh[patchi];
                 forAll(pp, i)
@@ -299,7 +354,7 @@ int main(int argc, char *argv[])
                 }
             }
             // Collect faces on faceZones
-            for (const label zonei : includeFaceZones)
+            for (const label zonei : faceZoneIds)
             {
                 const faceZone& pp = fzm[zonei];
                 forAll(pp, i)

applications/utilities/surface/surfaceTransformPoints/Make/options

@@ -1,5 +1,7 @@
 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):
+    The rollPitchYaw and yawPitchRoll options take three angles (degrees)
-    - roll (rotation about x) followed by
+    that describe the intrinsic Euler rotation.
-    - pitch (rotation about y) followed by
-    - yaw (rotation about z)
 
-    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 "axisAngleRotation.H"
-#include "unitConversion.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;
+        Info<< "Rotating points by " << rot << endl;
-
+        points = transform(rot, points);
-        points = transform(rotT, 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
+            << "    about " << axis << nl
-            << "    angle " << axisAngle.second() << nl;
+            << "    angle " << angle << nl;
 
-        const quaternion quat
+        const tensor rot(axisAngle::rotation(axis, angle, true));
-        (
-            axisAngle.first(),
-            degToRad(axisAngle.second())
-        );
 
-        Info<< "Rotating points by quaternion " << quat << endl;
+        Info<< "Rotating points by " << rot << endl;
-        points = transform(quat, points);
+        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 " << rot << endl;
-
+        points = transform(rot, points);
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, 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 tensor rot(euler::rotation(euler::eulerOrder::ZYX, v, true));
 
-        const quaternion quat(quaternion::rotationSequence::ZYX, v);
+        Info<< "Rotating points by " << rot << endl;
-
+        points = transform(rot, points);
-        Info<< "Rotating points by quaternion " << quat << endl;
-        points = transform(quat, 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
+    // Public Types
-    typedef scalar cmptType;
 
-    //- Euler-angle rotation sequence
+        //- Component type
-    enum rotationSequence
+        typedef scalar cmptType;
-    {
+
-        ZYX, ZYZ, ZXY, ZXZ, YXZ, YXY, YZX, YZY, XYZ, XYX, XZY, XZX
+        //- 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
+        //- Construct rotation quaternion given direction d and angle theta
-        //  and angle theta
         inline quaternion(const vector& d, const scalar theta);
 
-        //- Construct a rotation quaternion given the direction d
+        //- Construct a rotation quaternion given direction d
-        //- and cosine angle cosTheta and a if d is normalized
+        //- 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
+        //- Construct from three Euler rotation angles
-        inline quaternion(const rotationSequence rs, const vector& 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
+       //- Return the Euler rotation angles corresponding to the
-       //- specified rotation sequence
+       //- specified rotation order
-       inline vector eulerAngles(const rotationSequence rs) const;
+       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;
+        }
 
         case YZX:
+        {
             operator=(quaternion(vector(0, 1, 0), angles.x()));
             operator*=(quaternion(vector(0, 0, 1), angles.y()));
             operator*=(quaternion(vector(1, 0, 0), angles.z()));
             break;
+        }
 
         case YZY:
+        {
             operator=(quaternion(vector(0, 1, 0), angles.x()));
             operator*=(quaternion(vector(0, 0, 1), angles.y()));
             operator*=(quaternion(vector(0, 1, 0), angles.z()));
             break;
+        }
 
         case XYZ:
+        {
             operator=(quaternion(vector(1, 0, 0), angles.x()));
             operator*=(quaternion(vector(0, 1, 0), angles.y()));
             operator*=(quaternion(vector(0, 0, 1), angles.z()));
             break;
+        }
 
         case XYX:
+        {
             operator=(quaternion(vector(1, 0, 0), angles.x()));
             operator*=(quaternion(vector(0, 1, 0), angles.y()));
             operator*=(quaternion(vector(1, 0, 0), angles.z()));
             break;
+        }
 
         case XZY:
+        {
             operator=(quaternion(vector(1, 0, 0), angles.x()));
             operator*=(quaternion(vector(0, 0, 1), angles.y()));
             operator*=(quaternion(vector(0, 1, 0), angles.z()));
             break;
+        }
 
         case XZX:
+        {
             operator=(quaternion(vector(1, 0, 0), angles.x()));
             operator*=(quaternion(vector(0, 0, 1), angles.y()));
             operator*=(quaternion(vector(1, 0, 0), angles.z()));
             break;
+        }
 
         default:
             FatalErrorInFunction
-                << "Unknown rotation sequence " << rs << abort(FatalError);
+                << "Unknown euler rotation order "
+                << int(order) << abort(FatalError);
             break;
     }
 }
@@ -201,7 +233,7 @@ inline Foam::quaternion::quaternion
          && rotationTensor.xx() > rotationTensor.zz()
         )
         {
-            scalar s = 2.0*Foam::sqrt
+            const scalar s = 2.0*Foam::sqrt
             (
                 1.0
               + rotationTensor.xx()
@@ -219,7 +251,7 @@ inline Foam::quaternion::quaternion
             rotationTensor.yy() > rotationTensor.zz()
         )
         {
-            scalar s = 2.0*Foam::sqrt
+            const scalar s = 2.0*Foam::sqrt
             (
                 1.0
               + rotationTensor.yy()
@@ -234,7 +266,7 @@ inline Foam::quaternion::quaternion
         }
         else
         {
-            scalar s = 2.0*Foam::sqrt
+            const scalar s = 2.0*Foam::sqrt
             (
                 1.0
               + rotationTensor.zz()
@@ -373,7 +405,7 @@ inline Foam::vector Foam::quaternion::threeAxes
 
 inline Foam::vector Foam::quaternion::eulerAngles
 (
-    const rotationSequence rs
+    const eulerOrder order
 ) const
 {
     const scalar w2 = sqr(w());
@@ -381,9 +413,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
     const scalar y2 = sqr(v().y());
     const scalar z2 = sqr(v().z());
 
-    switch(rs)
+    switch (order)
     {
         case ZYX:
+        {
             return threeAxes
             (
                 2*(v().x()*v().y() + w()*v().z()),
@@ -393,8 +426,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 - x2 - y2 + z2
             );
             break;
+        }
 
         case ZYZ:
+        {
             return twoAxes
             (
                 2*(v().y()*v().z() - w()*v().x()),
@@ -404,8 +439,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(w()*v().y() - v().x()*v().z())
             );
             break;
+        }
 
         case ZXY:
+        {
             return threeAxes
             (
                 2*(w()*v().z() - v().x()*v().y()),
@@ -415,8 +452,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 - x2 - y2 + z2
             );
             break;
+        }
 
         case ZXZ:
+        {
             return twoAxes
             (
                 2*(v().x()*v().z() + w()*v().y()),
@@ -426,8 +465,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(v().y()*v().z() + w()*v().x())
             );
             break;
+        }
 
         case YXZ:
+        {
             return threeAxes
             (
                 2*(v().x()*v().z() + w()*v().y()),
@@ -437,8 +478,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 - x2 + y2 - z2
             );
             break;
+        }
 
         case YXY:
+        {
             return twoAxes
             (
                 2*(v().x()*v().y() - w()*v().z()),
@@ -448,8 +491,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(w()*v().x() - v().y()*v().z())
             );
             break;
+        }
 
         case YZX:
+        {
             return threeAxes
             (
                 2*(w()*v().y() - v().x()*v().z()),
@@ -459,8 +504,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 - x2 + y2 - z2
             );
             break;
+        }
 
         case YZY:
+        {
             return twoAxes
             (
                 2*(v().y()*v().z() + w()*v().x()),
@@ -470,8 +517,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(v().x()*v().y() + w()*v().z())
             );
             break;
+        }
 
         case XYZ:
+        {
             return threeAxes
             (
                 2*(w()*v().x() - v().y()*v().z()),
@@ -481,8 +530,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 + x2 - y2 - z2
             );
             break;
+        }
 
         case XYX:
+        {
             return twoAxes
             (
                 2*(v().x()*v().y() + w()*v().z()),
@@ -492,8 +543,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(v().x()*v().z() + w()*v().y())
             );
             break;
+        }
 
         case XZY:
+        {
             return threeAxes
             (
                 2*(v().y()*v().z() + w()*v().x()),
@@ -503,8 +556,10 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 w2 + x2 - y2 - z2
             );
             break;
+        }
 
         case XZX:
+        {
             return twoAxes
             (
                 2*(v().x()*v().z() - w()*v().y()),
@@ -514,12 +569,16 @@ inline Foam::vector Foam::quaternion::eulerAngles
                 2*(w()*v().z() - v().x()*v().y())
             );
             break;
+        }
+
         default:
             FatalErrorInFunction
-                << "Unknown rotation sequence " << rs << abort(FatalError);
+                << "Unknown euler rotation order "
-            return Zero;
+                << int(order) << abort(FatalError);
             break;
     }
+
+    return Zero;
 }
 
 

src/meshTools/coordinate/rotation/EulerCoordinateRotation.C

@@ -62,6 +62,7 @@ namespace Foam
 
 Foam::tensor Foam::coordinateRotations::euler::rotation
 (
+    const eulerOrder order,
     const vector& angles,
     bool degrees
 )
@@ -83,14 +84,163 @@ Foam::tensor Foam::coordinateRotations::euler::rotation
 
     // https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix
 
-    // Z1-X2-Z3 rotation
+    switch (order)
-    return
+    {
-        tensor
+        // Proper Euler angles
-        (
+
-            c1*c3 - c2*s1*s3, -c1*s3 - c2*c3*s1,  s1*s2,
+        case eulerOrder::XZX:  // X1-Z2-X3 rotation
-            c3*s1 + c1*c2*s3,  c1*c2*c3 - s1*s3, -c1*s2,
+        {
-            s2*s3,             c3*s2,             c2
+            return tensor
-        );
+            (
+                (   c2  ), (      -c3*s2      ), (      s2*s3        ),
+                ( c1*s2 ), ( c1*c2*c3 - s1*s3 ), ( -c3*s1 - c1*c2*s3 ),
+                ( s1*s2 ), ( c1*s3 + c2*c3*s1 ), (  c1*c3 - c2*s1*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::XYX:  // X1-Y2-X3 rotation
+        {
+            return tensor
+            (
+                (    c2  ), (      s2*s3       ), (      c3*s2 ),
+                (  s1*s2 ), ( c1*c3 - c2*s1*s3 ), ( -c1*s3 - c2*c3*s1 ),
+                ( -c1*s2 ), ( c3*s1 + c1*c2*s3 ), (  c1*c2*c3 - s1*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::YXY:  // Y1-X2-Y3 rotation
+        {
+            return tensor
+            (
+                ( c1*c3 - c2*s1*s3 ), ( s1*s2 ), ( c1*s3 + c2*c3*s1 ),
+                (     s2*s3        ), (   c2  ), (    -c3*s2        ),
+                ( -c3*s1 -c1*c2*s3 ), ( c1*s2 ), ( c1*c2*c3 - s1*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::YZY:  // Y1-Z2-Y3 rotation
+        {
+            return tensor
+            (
+                ( c1*c2*c3 - s1*s3 ), ( -c1*s2 ), ( c3*s1 + c1*c2*s3 ),
+                (        c3*s2     ), (   c2   ), (     s2*s3        ),
+                (-c1*s3 - c2*c3*s1 ), (  s1*s2 ), ( c1*c3 - c2*s1*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::ZYZ:  // Z1-Y2-Z3 rotation
+        {
+            return tensor
+            (
+                ( c1*c2*c3 - s1*s3 ), ( -c3*s1 - c1*c2*s3 ), ( c1*s2 ),
+                ( c1*s3 + c2*c3*s1 ), (  c1*c3 - c2*s1*s3 ), ( s1*s2 ),
+                (     -c3*s2       ), (      s2*s3 ),        (   c2  )
+            );
+            break;
+        }
+
+        case eulerOrder::ZXZ:  // Z1-X2-Z3 rotation
+        {
+            return tensor
+            (
+                ( c1*c3 - c2*s1*s3 ), ( -c1*s3 - c2*c3*s1 ), (  s1*s2 ),
+                ( c3*s1 + c1*c2*s3 ), ( c1*c2*c3 - s1*s3  ), ( -c1*s2 ),
+                (     s2*s3        ), (      c3*s2        ), (    c2  )
+            );
+            break;
+        }
+
+
+            // Tait-Bryan angles
+
+        case eulerOrder::XZY:  // X1-Z2-Y3 rotation
+        {
+            return tensor
+            (
+                (      c2*c3       ), (  -s2  ), (       c2*s3      ),
+                ( s1*s3 + c1*c3*s2 ), ( c1*c2 ), ( c1*s2*s3 - c3*s1 ),
+                ( c3*s1*s2 - c1*s3 ), ( c2*s1 ), ( c1*c3 + s1*s2*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::XYZ:  // X1-Y2-Z3 rotation
+        {
+            return tensor
+            (
+                (      c2*c3       ), (    -c2*s3        ), (    s2  ),
+                ( c1*s3 + c3*s1*s2 ), ( c1*c3 - s1*s2*s3 ), ( -c2*s1 ),
+                ( s1*s3 - c1*c3*s2 ), ( c3*s1 + c1*s2*s3 ), (  c1*c2 )
+            );
+            break;
+        }
+
+        case eulerOrder::YXZ:  // Y1-X2-Z3 rotation
+        {
+            return tensor
+            (
+                ( c1*c3 + s1*s2*s3 ), ( c3*s1*s2 - c1*s3 ), ( c2*s1 ),
+                (     c2*s3        ), (        c2*c3     ), (  -s2  ),
+                ( c1*s2*s3 - c3*s1 ), ( c1*c3*s2 + s1*s3 ), ( c1*c2 )
+            );
+            break;
+        }
+
+        case eulerOrder::YZX:  // Y1-Z2-X3 rotation
+        {
+            return tensor
+            (
+                (  c1*c2 ), ( s1*s3 - c1*c3*s2 ), ( c3*s1 + c1*s2*s3 ),
+                (  s2    ), ( c2*c3            ), ( -c2*s3           ),
+                ( -c2*s1 ), ( c1*s3 + c3*s1*s2 ), ( c1*c3 - s1*s2*s3 )
+            );
+            break;
+        }
+
+        case eulerOrder::ZYX:  // Z1-Y2-X3 rotation
+        {
+            return tensor
+            (
+                ( c1*c2 ), ( c1*s2*s3 - c3*s1 ), ( s1*s3 + c1*c3*s2 ),
+                ( c2*s1 ), ( c1*c3 + s1*s2*s3 ), ( c3*s1*s2 - c1*s3 ),
+                (  -s2  ), (      c2*s3       ), (        c2*c3     )
+            );
+            break;
+        }
+
+        case eulerOrder::ZXY:  // Z1-X2-Y3 rotation
+        {
+            return tensor
+            (
+                ( c1*c3 - s1*s2*s3 ), ( -c2*s1 ), ( c1*s3 + c3*s1*s2 ),
+                ( c3*s1 + c1*s2*s3 ), (  c1*c2 ), ( s1*s3 - c1*c3*s2 ),
+                (    -c2*s3 ),        (   s2   ), (     c2*c3        )
+            );
+            break;
+        }
+
+        default:
+            FatalErrorInFunction
+                << "Unknown euler rotation order "
+                << int(order) << abort(FatalError);
+            break;
+    }
+
+    return tensor::I;
+}
+
+
+Foam::tensor Foam::coordinateRotations::euler::rotation
+(
+    const vector& angles,
+    bool degrees
+)
+{
+    return rotation(eulerOrder::ZXZ, angles, degrees);
 }
 
 
@@ -100,7 +250,8 @@ Foam::coordinateRotations::euler::euler()
 :
     coordinateRotation(),
     angles_(Zero),
-    degrees_(true)
+    degrees_(true),
+    order_(eulerOrder::ZXZ)
 {}
 
 
@@ -108,7 +259,8 @@ Foam::coordinateRotations::euler::euler(const euler& crot)
 :
     coordinateRotation(crot),
     angles_(crot.angles_),
-    degrees_(crot.degrees_)
+    degrees_(crot.degrees_),
+    order_(crot.order_)
 {}
 
 
@@ -120,7 +272,8 @@ Foam::coordinateRotations::euler::euler
 :
     coordinateRotation(),
     angles_(angles),
-    degrees_(degrees)
+    degrees_(degrees),
+    order_(eulerOrder::ZXZ)
 {}
 
 
@@ -134,7 +287,8 @@ Foam::coordinateRotations::euler::euler
 :
     coordinateRotation(),
     angles_(angle1, angle2, angle3),
-    degrees_(degrees)
+    degrees_(degrees),
+    order_(eulerOrder::ZXZ)
 {}
 
 
@@ -142,7 +296,16 @@ Foam::coordinateRotations::euler::euler(const dictionary& dict)
 :
     coordinateRotation(),
     angles_(dict.getCompat<vector>("angles", {{"rotation", 1806}})),
-    degrees_(dict.lookupOrDefault("degrees", true))
+    degrees_(dict.lookupOrDefault("degrees", true)),
+    order_
+    (
+        quaternion::eulerOrderNames.lookupOrDefault
+        (
+            "order",
+            dict,
+            quaternion::eulerOrder::ZXZ
+        )
+    )
 {}
 
 
@@ -182,6 +345,12 @@ void Foam::coordinateRotations::euler::writeEntry
         os.writeEntry("degrees", "false");
     }
 
+    // writeEntryIfDifferent, but with enumerated name
+    if (order_ != eulerOrder::ZXZ)
+    {
+        os.writeEntry("order", quaternion::eulerOrderNames[order_]);
+    }
+
     os.endBlock();
 }
 

src/meshTools/coordinate/rotation/EulerCoordinateRotation.H

@@ -61,6 +61,7 @@ SourceFiles
 #define coordinateRotations_euler_H
 
 #include "coordinateRotation.H"
+#include "quaternion.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,6 +78,16 @@ class euler
 :
     public coordinateRotation
 {
+public:
+
+    // Public Types
+
+    //- Euler-angle rotation order
+    using eulerOrder = quaternion::eulerOrder;
+
+
+private:
+
     // Private Data
 
         //- The rotation angles
@@ -85,6 +96,9 @@ class euler
         //- Angles measured in degrees
         bool degrees_;
 
+        //- The Euler-angle rotation order (default: zxz)
+        eulerOrder order_;
+
 
 public:
 
@@ -126,7 +140,15 @@ public:
 
         //- The rotation tensor calculated for the intrinsic Euler
         //- angles in z-x-z order
-        static tensor rotation(const vector& angles, bool degrees);
+        static tensor rotation(const vector& angles, bool degrees=false);
+
+        //- The rotation tensor calculated for given angles and order
+        static tensor rotation
+        (
+            const eulerOrder order,
+            const vector& angles,
+            bool degrees=false
+        );
 
 
     // Member Functions

src/meshTools/coordinate/rotation/axisAngleRotation.C

@@ -57,6 +57,24 @@ void Foam::coordinateRotations::axisAngle::checkSpec()
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::tensor Foam::coordinateRotations::axisAngle::rotation
+(
+    const vector& axis,
+    const scalar angle,
+    bool degrees
+)
+{
+    if (mag(angle) < VSMALL || mag(axis) < SMALL)
+    {
+        return sphericalTensor::I;  // identity rotation
+    }
+
+    return quaternion(axis, (degrees ? degToRad(angle) : angle)).R();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::coordinateRotations::axisAngle::axisAngle()
@@ -144,12 +162,7 @@ void Foam::coordinateRotations::axisAngle::clear()
 
 Foam::tensor Foam::coordinateRotations::axisAngle::R() const
 {
-    if (mag(angle_) < VSMALL || mag(axis_) < SMALL)
+    return rotation(axis_, angle_, degrees_);
-    {
-        return sphericalTensor::I; // identity rotation
-    }
-
-    return quaternion(axis_, (degrees_ ? degToRad(angle_) : angle_)).R();
 }
 
 

src/meshTools/coordinate/rotation/axisAngleRotation.H

@@ -130,6 +130,17 @@ public:
     virtual ~axisAngle() = default;
 
 
+    // Static Member Functions
+
+        //- The rotation tensor for given axis/angle
+        static tensor rotation
+        (
+            const vector& axis,
+            const scalar angle,
+            bool degrees=false
+        );
+
+
     // Member Functions
 
         //- Reset specification

tutorials/incompressible/lumpedPointMotion/building/steady/system/controlDict

@@ -14,7 +14,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-libs            ("libmeshTools.so" "liblumpedPointMotion.so");
+libs
+(
+    "libmeshTools.so"
+    "liblumpedPointMotion.so"
+    "libfvMotionSolvers.so"
+);
 
 application     simpleFoam;     // Change to pimpleFoam for transient