diff --git a/applications/test/quaternion/Test-quaternion.C b/applications/test/quaternion/Test-quaternion.C
index 1317ad4bf330258a72b0a5f72a32fa231cee1cab..852a31fed17d3b6bb030f3f1acca8dcb75da3a24 100644
--- a/applications/test/quaternion/Test-quaternion.C
+++ b/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);
}
}
diff --git a/applications/test/sizeof/Make/options b/applications/test/sizeof/Make/options
index 4e772fdf9d7bc94221d127458f9d2ca32850fe69..410c857f852235bdeccb91e12e9e1b1b919ab15b 100644
--- a/applications/test/sizeof/Make/options
+++ b/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 */
diff --git a/applications/utilities/mesh/manipulation/transformPoints/Make/options b/applications/utilities/mesh/manipulation/transformPoints/Make/options
index 16eb624becb6894353d6c74bb8bab5de37551b56..969020c4afaf5d784299462b9e1af282040ba6b4 100644
--- a/applications/utilities/mesh/manipulation/transformPoints/Make/options
+++ b/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
diff --git a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
index 698ec7f601d811eaebea0a70502c88f1f1e3abfb..f19e726fc8d1b1f811b514b783262147cc660bc8 100644
--- a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
+++ b/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;
-
- 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);
}
}
diff --git a/applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C b/applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C
index 7232b683f45b68edae61ee78a89e4c09be9794d8..6e3c8f5aec52b80f17dd39409473484e4c4c035f 100644
--- a/applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C
+++ b/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"
- "Eg, 'cells' or '( slice \"mfp-.*\" )'"
+ "Specify single patch or multiple patches to exclude from writing."
+ " 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;
-
- if (args.found("patches"))
- {
- includePatches =
- bMesh.patchSet(args.getList<wordRe>("patches")).sortedToc();
- }
- else if (includeProcPatches)
- {
- includePatches = identity(bMesh.size());
- }
- else
- {
- includePatches = identity(bMesh.nNonProcessor());
- }
-
+ labelList patchIds =
+ (
+ (includePatches.size() || excludePatches.size())
+ ? getSelectedPatches(bMesh, includePatches, excludePatches)
+ : includeProcPatches
+ ? identity(bMesh.size())
+ : identity(bMesh.nNonProcessor())
+ );
- labelList includeFaceZones;
+ labelList faceZoneIds;
const faceZoneMesh& fzm = mesh.faceZones();
- if (args.found("faceZones"))
+ if (selectedFaceZones.size())
{
- const wordList allZoneNames(fzm.names());
-
- 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();
+ faceZoneIds = fzm.indices(selectedFaceZones);
Info<< "Additionally extracting faceZones "
- << UIndirectList<word>(allZoneNames, includeFaceZones)
- << endl;
+ << fzm.names(selectedFaceZones) << nl;
}
@@ -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)
diff --git a/applications/utilities/surface/surfaceTransformPoints/Make/options b/applications/utilities/surface/surfaceTransformPoints/Make/options
index a504dd8617bfa4ff4d8fc1075d5d68d6db66db89..5ae72ebd31bdf7500846ee0dca7a4c123ff4fb43 100644
--- a/applications/utilities/surface/surfaceTransformPoints/Make/options
+++ b/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
diff --git a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
index 229df4d1b9bb4feae949ca0c5a15f6536c59c7c0..552863b9b45340d4495cf4119597f8accfc8d861 100644
--- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
+++ b/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;
diff --git a/src/OpenFOAM/primitives/quaternion/quaternion.C b/src/OpenFOAM/primitives/quaternion/quaternion.C
index a2b53b8391593c4dca0ba55f464382f43f9663da..9ced3d5317fff96ecdcd3e09f63b8c426372cd5a 100644
--- a/src/OpenFOAM/primitives/quaternion/quaternion.C
+++ b/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)
diff --git a/src/OpenFOAM/primitives/quaternion/quaternion.H b/src/OpenFOAM/primitives/quaternion/quaternion.H
index 339253e34bc605c7bf9bb647f2374ba099015ffb..42040e0c252729774181b7ab52517ee1398351b0 100644
--- a/src/OpenFOAM/primitives/quaternion/quaternion.H
+++ b/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
diff --git a/src/OpenFOAM/primitives/quaternion/quaternionI.H b/src/OpenFOAM/primitives/quaternion/quaternionI.H
index 13b4f70d9423b794228045e3bfb8e27f182d7daa..e692707c3011ce66ef8721e43bc4b267f45b5302 100644
--- a/src/OpenFOAM/primitives/quaternion/quaternionI.H
+++ b/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);
- return Zero;
+ << "Unknown euler rotation order "
+ << int(order) << abort(FatalError);
break;
}
+
+ return Zero;
}
diff --git a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C
index 98104f9c03045ea95211850ecc91219a401183c3..0de692db9ccbdb2b53e9f15323333844474546c2 100644
--- a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C
+++ b/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
- 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
- );
+ switch (order)
+ {
+ // Proper Euler angles
+
+ case eulerOrder::XZX: // X1-Z2-X3 rotation
+ {
+ 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();
}
diff --git a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H
index 05ff2d82ac2ed599070ae0d917d8df284d21640b..7d15ee65fce103a222f3d08ab984cc713a62e0ae 100644
--- a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H
+++ b/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
diff --git a/src/meshTools/coordinate/rotation/axisAngleRotation.C b/src/meshTools/coordinate/rotation/axisAngleRotation.C
index 1e88f383f61a6260794f250b5c0bfe41864daf15..430efd5bad25d7161788a01e258515b574a4cd8f 100644
--- a/src/meshTools/coordinate/rotation/axisAngleRotation.C
+++ b/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 sphericalTensor::I; // identity rotation
- }
-
- return quaternion(axis_, (degrees_ ? degToRad(angle_) : angle_)).R();
+ return rotation(axis_, angle_, degrees_);
}
diff --git a/src/meshTools/coordinate/rotation/axisAngleRotation.H b/src/meshTools/coordinate/rotation/axisAngleRotation.H
index 0b5533331e5cb860bdd334fae95059d45496124b..a4ff281b123d5ef35498aa120ddf39b519cdb724 100644
--- a/src/meshTools/coordinate/rotation/axisAngleRotation.H
+++ b/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
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/system/controlDict b/tutorials/incompressible/lumpedPointMotion/building/steady/system/controlDict
index 0f76b2efedb42d0abb02e578695b965ea71ff703..cfa0acb757cb87184426e3c46479d6647ac639fe 100644
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/system/controlDict
+++ b/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