Commit 00325db3 authored by Mark OLESEN's avatar Mark OLESEN
Browse files

ENH: additional options for transforming points (closes #660)

- The -rotate-angle option allows convenient specification of a
  rotation about an arbitrary axis. Eg, -rotate-angle '((1 1 1) 45)'

- The -origin option can be used to temporarily shift the origin
  for the rotation operations. For example,

      -origin '(0 0 1)' -rotate-angle '((1 0 0) 180)'

  for mirroring.
parent 95e357ed
......@@ -29,9 +29,13 @@ Description
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "quaternion.H"
#include "septernion.H"
#include "mathematicalConstants.H"
#include "Tuple2.H"
#include "IOstreams.H"
#include "transform.H"
using namespace Foam;
......@@ -39,6 +43,95 @@ using namespace Foam;
int main(int argc, char *argv[])
{
argList::addOption
(
"rollPitchYaw",
"vector",
"Rotate by '(roll pitch yaw)' in degrees"
);
argList::addOption
(
"yawPitchRoll",
"vector",
"Rotate by '(yaw pitch roll)' in degrees"
);
argList::addOption
(
"rotate-angle",
"(vector angle)",
"Rotate about the <vector> by <angle> degrees - eg, '((1 0 0) 45)'"
);
argList args(argc, argv);
vector rotVector;
if (args.optionReadIfPresent("rollPitchYaw", rotVector))
{
Info<< "Rotate by" << nl
<< " roll " << rotVector.x() << nl
<< " pitch " << rotVector.y() << nl
<< " yaw " << rotVector.z() << nl;
// degToRad
rotVector *= constant::mathematical::pi/180.0;
const quaternion quat(quaternion::rotationSequence::XYZ, rotVector);
Info<< "quaternion " << quat << endl;
Info<< "rotation = " << quat.R() << endl;
}
if (args.optionReadIfPresent("yawPitchRoll", rotVector))
{
Info<< "Rotate by" << nl
<< " yaw " << rotVector.x() << nl
<< " pitch " << rotVector.y() << nl
<< " roll " << rotVector.z() << nl;
// degToRad
rotVector *= constant::mathematical::pi/180.0;
const quaternion quat(quaternion::rotationSequence::ZYX, rotVector);
Info<< "quaternion " << quat << endl;
Info<< "rotation = " << quat.R() << endl;
}
if (args.optionFound("rotate-angle"))
{
const Tuple2<vector, scalar> axisAngle
(
args.optionLookup("rotate-angle")()
);
Info<< "Rotate" << nl
<< " about " << axisAngle.first() << nl
<< " angle " << axisAngle.second() << nl;
const quaternion quat
(
axisAngle.first(),
axisAngle.second() * constant::mathematical::pi/180.0 // degToRad
);
Info<< "quaternion " << quat << endl;
Info<< "rotation = " << quat.R() << endl;
Info<< "transform Ra = "
<< Ra
(
axisAngle.first() / mag(axisAngle.first()),
axisAngle.second() * constant::mathematical::pi/180.0
) << endl;
Info<< "-ve Ra = "
<< Ra
(
axisAngle.first() / mag(axisAngle.first()),
-axisAngle.second() * constant::mathematical::pi/180.0
) << endl;
}
Info<< nl << nl;
quaternion q(vector(1, 2, 3), 0.7853981);
Info<< "q " << q << endl;
......
......@@ -48,19 +48,19 @@ void RotateFields
(
const fvMesh& mesh,
const IOobjectList& objects,
const tensor& T
const tensor& rotT
)
{
// Search list of objects for volScalarFields
// Objects of field type
IOobjectList fields(objects.lookupClass(GeometricField::typeName));
forAllIter(IOobjectList, fields, fieldIter)
{
Info<< " Rotating " << fieldIter()->name() << endl;
GeometricField theta(*fieldIter(), mesh);
transform(theta, dimensionedTensor(T), theta);
theta.write();
GeometricField fld(*fieldIter(), mesh);
transform(fld, dimensionedTensor(rotT), fld);
fld.write();
}
}
......@@ -69,6 +69,12 @@ void RotateFields
int main(int argc, char *argv[])
{
argList::addNote
(
"Rotate mesh points and vector/tensor fields\n"
"Rotation from the <n1> vector to the <n2> vector"
);
timeSelector::addOptions();
argList::addArgument("n1");
......@@ -83,7 +89,7 @@ int main(int argc, char *argv[])
vector n2(args.argRead<vector>(2));
n2 /= mag(n2);
tensor T(rotationTensor(n1, n2));
const tensor rotT(rotationTensor(n1, n2));
{
pointIOField points
......@@ -100,7 +106,7 @@ int main(int argc, char *argv[])
)
);
points = transform(T, points);
points = transform(rotT, points);
// Set the precision of the points data to 10
IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
......@@ -123,15 +129,15 @@ int main(int argc, char *argv[])
// Search for list of objects for this time
IOobjectList objects(mesh, runTime.timeName());
RotateFields<volVectorField>(mesh, objects, T);
RotateFields<volSphericalTensorField>(mesh, objects, T);
RotateFields<volSymmTensorField>(mesh, objects, T);
RotateFields<volTensorField>(mesh, objects, T);
RotateFields<volVectorField>(mesh, objects, rotT);
RotateFields<volSphericalTensorField>(mesh, objects, rotT);
RotateFields<volSymmTensorField>(mesh, objects, rotT);
RotateFields<volTensorField>(mesh, objects, rotT);
RotateFields<surfaceVectorField>(mesh, objects, T);
RotateFields<surfaceSphericalTensorField>(mesh, objects, T);
RotateFields<surfaceSymmTensorField>(mesh, objects, T);
RotateFields<surfaceTensorField>(mesh, objects, T);
RotateFields<surfaceVectorField>(mesh, objects, rotT);
RotateFields<surfaceSphericalTensorField>(mesh, objects, rotT);
RotateFields<surfaceSymmTensorField>(mesh, objects, rotT);
RotateFields<surfaceTensorField>(mesh, objects, rotT);
}
Info<< "End\n" << endl;
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -40,6 +40,9 @@ Usage
-rotate (vector vector)
Rotates the points from the first vector to the second,
-rotate-angle (vector angle)
Rotate angle degrees about vector axis.
or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
......@@ -68,7 +71,6 @@ Usage
#include "pointFields.H"
#include "transformField.H"
#include "transformGeometricField.H"
#include "StringStream.H"
#include "mathematicalConstants.H"
using namespace Foam;
......@@ -81,7 +83,7 @@ void readAndRotateFields
(
PtrList<GeoField>& flds,
const fvMesh& mesh,
const tensor& T,
const tensor& rotT,
const IOobjectList& objects
)
{
......@@ -89,7 +91,7 @@ void readAndRotateFields
forAll(flds, i)
{
Info<< "Transforming " << flds[i].name() << endl;
dimensionedTensor dimT("t", flds[i].dimensions(), T);
const dimensionedTensor dimT("t", flds[i].dimensions(), rotT);
transform(flds[i], dimT, flds[i]);
}
}
......@@ -155,23 +157,34 @@ int main(int argc, char *argv[])
"Translate by specified <vector> - eg, '(1 0 0)' before rotations"
);
argList::addOption
(
"origin",
"point",
"Use specified <point> as origin for rotations"
);
argList::addOption
(
"rotate",
"(vectorA vectorB)",
"Transform as a rotation between <vectorA> and <vectorB> "
"- eg, '( (1 0 0) (0 0 1) )'"
"Rotate from <vectorA> to <vectorB> - eg, '((1 0 0) (0 0 1))'"
);
argList::addOption
(
"rotate-angle",
"(vector scalar)",
"Rotate <angle> degrees about <vector> - eg, '((1 0 0) 45)'"
);
argList::addOption
(
"rollPitchYaw",
"vector",
"Rotate by '(roll pitch yaw)' in degrees"
"Rotate by '(roll pitch yaw)' degrees"
);
argList::addOption
(
"yawPitchRoll",
"vector",
"Rotate by '(yaw pitch roll)' in degrees"
"Rotate by '(yaw pitch roll)' degrees"
);
argList::addBoolOption
(
......@@ -182,25 +195,54 @@ int main(int argc, char *argv[])
(
"scale",
"scalar | vector",
"Scale by the specified amount - eg, for a uniform [mm] to [m] scaling "
"use either (0.001 0.001 0.001)' or simply '0.001'"
"Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
"use either '(0.001 0.001 0.001)' or simply '0.001'"
);
#include "addRegionOption.H"
#include "setRootCase.H"
const bool doRotateFields = args.optionFound("rotateFields");
// Verify that an operation has been specified
{
const List<word> operationNames
{
"translate",
"rotate",
"rotate-angle",
"rollPitchYaw",
"yawPitchRoll",
"scale"
};
if (!args.optionCount(operationNames))
{
FatalError
<< "No operation supplied, "
<< "use least one of the following:" << nl
<< " ";
for (const auto& opName : operationNames)
{
FatalError
<< " -" << opName;
}
FatalError
<< nl << exit(FatalError);
}
}
#include "createTime.H"
word regionName = polyMesh::defaultRegion;
fileName meshDir;
fileName meshDir = polyMesh::meshSubDir;
if (args.optionReadIfPresent("region", regionName))
{
meshDir = regionName/polyMesh::meshSubDir;
}
else
{
meshDir = polyMesh::meshSubDir;
}
pointIOField points
(
......@@ -216,17 +258,6 @@ int main(int argc, char *argv[])
)
);
const bool doRotateFields = args.optionFound("rotateFields");
// this is not actually stringent enough:
if (args.options().empty())
{
FatalErrorInFunction
<< "No options supplied, please use one or more of "
"-translate, -rotate or -scale options."
<< exit(FatalError);
}
vector v;
if (args.optionReadIfPresent("translate", v))
{
......@@ -235,6 +266,14 @@ int main(int argc, char *argv[])
points += v;
}
vector origin;
const bool useOrigin = args.optionReadIfPresent("origin", origin);
if (useOrigin)
{
Info<< "Set origin for rotations to " << origin << endl;
points -= origin;
}
if (args.optionFound("rotate"))
{
Pair<vector> n1n2
......@@ -243,15 +282,41 @@ int main(int argc, char *argv[])
);
n1n2[0] /= mag(n1n2[0]);
n1n2[1] /= mag(n1n2[1]);
tensor T = rotationTensor(n1n2[0], n1n2[1]);
Info<< "Rotating points by " << T << endl;
const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
Info<< "Rotating points by " << rotT << endl;
points = transform(rotT, points);
if (doRotateFields)
{
rotateFields(args, runTime, rotT);
}
}
else if (args.optionFound("rotate-angle"))
{
const Tuple2<vector, scalar> axisAngle
(
args.optionLookup("rotate-angle")()
);
Info<< "Rotating points " << nl
<< " about " << axisAngle.first() << nl
<< " angle " << axisAngle.second() << nl;
const quaternion quat
(
axisAngle.first(),
axisAngle.second() * pi/180.0 // degToRad
);
points = transform(T, points);
Info<< "Rotating points by quaternion " << quat << endl;
points = transform(quat, points);
if (doRotateFields)
{
rotateFields(args, runTime, T);
rotateFields(args, runTime, quat.R());
}
}
else if (args.optionReadIfPresent("rollPitchYaw", v))
......@@ -264,14 +329,14 @@ int main(int argc, char *argv[])
// degToRad
v *= pi/180.0;
const quaternion R(quaternion::rotationSequence::XYZ, v);
const quaternion quat(quaternion::rotationSequence::XYZ, v);
Info<< "Rotating points by quaternion " << R << endl;
points = transform(R, points);
Info<< "Rotating points by quaternion " << quat << endl;
points = transform(quat, points);
if (doRotateFields)
{
rotateFields(args, runTime, R.R());
rotateFields(args, runTime, quat.R());
}
}
else if (args.optionReadIfPresent("yawPitchRoll", v))
......@@ -284,14 +349,14 @@ int main(int argc, char *argv[])
// degToRad
v *= pi/180.0;
const quaternion R(quaternion::rotationSequence::ZYX, v);
const quaternion quat(quaternion::rotationSequence::ZYX, v);
Info<< "Rotating points by quaternion " << R << endl;
points = transform(R, points);
Info<< "Rotating points by quaternion " << quat << endl;
points = transform(quat, points);
if (doRotateFields)
{
rotateFields(args, runTime, R.R());
rotateFields(args, runTime, quat.R());
}
}
......@@ -325,6 +390,13 @@ int main(int argc, char *argv[])
}
}
if (useOrigin)
{
Info<< "Unset origin for rotations from " << origin << endl;
points += origin;
}
// Set the precision of the points data to 10
IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -45,6 +45,7 @@ Description
#include "boundBox.H"
#include "transformField.H"
#include "Pair.H"
#include "Tuple2.H"
#include "quaternion.H"
#include "mathematicalConstants.H"
......@@ -74,47 +75,80 @@ int main(int argc, char *argv[])
"Translate by specified <vector> - eg, '(1 0 0)' before rotations"
);
argList::addOption
(
"origin",
"point",
"Use specified <point> as origin for rotations"
);
argList::addOption
(
"rotate",
"(vectorA vectorB)",
"Transform as a rotation between <vectorA> and <vectorB> "
"- eg, '( (1 0 0) (0 0 1) )'"
"Rotate from <vectorA> to <vectorB> - eg, '((1 0 0) (0 0 1))'"
);
argList::addOption
(
"scale",
"scalar | vector",
"Scale by the specified amount - eg, for a uniform [mm] to [m] scaling "
"use either (0.001 0.001 0.001)' or simply '0.001'"
"rotate-angle",
"(vector scalar)",
"Rotate <angle> degrees about <vector> - eg, '((1 0 0) 45)'"
);
argList::addOption
(
"rollPitchYaw",
"vector",
"Rotate by '(roll pitch yaw)' in degrees"
"Rotate by '(roll pitch yaw)' degrees"
);
argList::addOption
(
"yawPitchRoll",
"vector",
"Rotate by '(yaw pitch roll)' in degrees"
"Rotate by '(yaw pitch roll)' degrees"
);
argList::addOption
(
"scale",
"scalar | vector",
"Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
"use either '(0.001 0.001 0.001)' or simply '0.001'"
);
argList args(argc, argv);
// Verify that an operation has been specified
{
const List<word> operationNames
{
"translate",
"rotate",
"rotate-angle",
"rollPitchYaw",
"yawPitchRoll",
"scale"
};
if (!args.optionCount(operationNames))
{
FatalError
<< "No operation supplied, "
<< "use least one of the following:" << nl
<< " ";
for (const auto& opName : operationNames)
{
FatalError
<< " -" << opName;
}
FatalError
<< nl << exit(FatalError);
}
}
const fileName surfFileName = args[1];
const fileName outFileName = args[2];
Info<< "Reading surf from " << surfFileName << " ..." << nl
<< "Writing surf to " << outFileName << " ..." << endl;
if (args.options().empty())
{
FatalErrorInFunction
<< "No options supplied, please use one or more of "
"-translate, -rotate or -scale options."
<< exit(FatalError);
}
meshedSurface surf1(surfFileName);
pointField points(surf1.points());
......@@ -127,6 +161,14 @@ int main(int argc, char *argv[])
points += v;
}
vector origin;
const bool useOrigin = args.optionReadIfPresent("origin", origin);
if (useOrigin)
{
Info<< "Set origin for rotations to " << origin << endl;
points -= origin;
}
if (args.optionFound("rotate"))
{
Pair<vector> n1n2
......@@ -136,11 +178,31 @@ int main(int argc, char *argv[])
n1n2[0] /= mag(n1n2[0]);
n1n2[1] /= mag(n1n2[1]);
const tensor T = rotationTensor(n1n2[0], n1n2[1]);
const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
Info<< "Rotating points by " << T << endl;
Info<< "Rotating points by " << rotT << endl;
points = transform(T, points);
points = transform(rotT, points);
}
else if (