From b9c738dd1039bf54f0418f141e885d69f07554bc Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Fri, 5 Oct 2018 00:49:55 +0200 Subject: [PATCH] ENH: expose calculation of coordinate rotations as static methods - improve handling of degenerate cases for the two-axes specification. --- .../rotation/EulerCoordinateRotation.C | 138 ++++----- .../rotation/EulerCoordinateRotation.H | 30 +- .../rotation/STARCDCoordinateRotation.C | 84 ++---- .../rotation/STARCDCoordinateRotation.H | 26 +- .../coordinate/rotation/axesRotation.C | 278 ++++++++++-------- .../coordinate/rotation/axesRotation.H | 26 +- .../coordinate/rotation/coordinateRotation.C | 70 +++-- .../coordinate/rotation/coordinateRotation.H | 11 +- 8 files changed, 339 insertions(+), 324 deletions(-) diff --git a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C index 5a7d86e828..196a25a923 100644 --- a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C +++ b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -24,8 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "EulerCoordinateRotation.H" - -#include "mathematicalConstants.H" +#include "unitConversion.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -47,6 +46,43 @@ namespace Foam ); } +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::tensor Foam::EulerCoordinateRotation::rotation +( + const vector& angles, + bool degrees +) +{ + scalar phi = angles.component(vector::X); // 1. Rotate about Z + scalar theta = angles.component(vector::Y); // 2. Rotate about X + scalar psi = angles.component(vector::Z); // 3. Rotate about Z + + if (degrees) + { + phi *= degToRad(); + theta *= degToRad(); + psi *= degToRad(); + } + + return + tensor + ( + cos(phi)*cos(psi) - sin(phi)*sin(psi)*cos(theta), + -sin(phi)*cos(psi)*cos(theta) - cos(phi)*sin(psi), + sin(phi)*sin(theta), + + cos(phi)*sin(psi)*cos(theta) + sin(phi)*cos(psi), + cos(phi)*cos(psi)*cos(theta) - sin(phi)*sin(psi), + -cos(phi)*sin(theta), + + sin(psi)*sin(theta), + cos(psi)*sin(theta), + cos(theta) + ); +} + + // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // Foam::vector Foam::EulerCoordinateRotation::transform(const vector& st) const @@ -147,49 +183,6 @@ Foam::symmTensor Foam::EulerCoordinateRotation::transformVector } -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::EulerCoordinateRotation::calcTransform -( - const scalar phiAngle, - const scalar thetaAngle, - const scalar psiAngle, - const bool inDegrees -) -{ - scalar phi = phiAngle; - scalar theta = thetaAngle; - scalar psi = psiAngle; - - if (inDegrees) - { - phi *= constant::mathematical::pi/180.0; - theta *= constant::mathematical::pi/180.0; - psi *= constant::mathematical::pi/180.0; - } - - R_ = - ( - tensor - ( - cos(phi)*cos(psi) - sin(phi)*sin(psi)*cos(theta), - -sin(phi)*cos(psi)*cos(theta) - cos(phi)*sin(psi), - sin(phi)*sin(theta), - - cos(phi)*sin(psi)*cos(theta) + sin(phi)*cos(psi), - cos(phi)*cos(psi)*cos(theta) - sin(phi)*sin(psi), - -cos(phi)*sin(theta), - - sin(psi)*sin(theta), - cos(psi)*sin(theta), - cos(theta) - ) - ); - - Rtr_ = R_.T(); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::EulerCoordinateRotation::EulerCoordinateRotation() @@ -212,35 +205,25 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation Foam::EulerCoordinateRotation::EulerCoordinateRotation ( const vector& phiThetaPsi, - const bool inDegrees + const bool degrees ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - calcTransform - ( - phiThetaPsi.component(vector::X), - phiThetaPsi.component(vector::Y), - phiThetaPsi.component(vector::Z), - inDegrees - ); -} + R_(rotation(phiThetaPsi, degrees)), + Rtr_(R_.T()) +{} Foam::EulerCoordinateRotation::EulerCoordinateRotation ( - const scalar phiAngle, - const scalar thetaAngle, - const scalar psiAngle, - const bool inDegrees + const scalar phi, + const scalar theta, + const scalar psi, + const bool degrees ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - calcTransform(phiAngle, thetaAngle, psiAngle, inDegrees); -} + R_(rotation(vector(phi, theta, psi), degrees)), + Rtr_(R_.T()) +{} Foam::EulerCoordinateRotation::EulerCoordinateRotation @@ -248,19 +231,16 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation const dictionary& dict ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - const vector rotation(dict.lookup("rotation")); - - calcTransform + R_ ( - rotation.component(vector::X), - rotation.component(vector::Y), - rotation.component(vector::Z), - dict.lookupOrDefault("degrees", true) - ); -} + rotation + ( + dict.get("rotation"), + dict.lookupOrDefault("degrees", true) + ) + ), + Rtr_(R_.T()) +{} Foam::EulerCoordinateRotation::EulerCoordinateRotation diff --git a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H index 6b2c44f809..57dfd31c46 100644 --- a/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H +++ b/src/meshTools/coordinate/rotation/EulerCoordinateRotation.H @@ -73,7 +73,6 @@ class EulerCoordinateRotation : public coordinateRotation { - // Private Member Data //- Local-to-global transformation tensor @@ -83,18 +82,6 @@ class EulerCoordinateRotation tensor Rtr_; - // Private Member Functions - - //- Calculate transformation tensor - void calcTransform - ( - const scalar phiAngle, - const scalar thetaAngle, - const scalar psiAngle, - const bool inDegrees - ); - - public: //- Runtime type information @@ -113,16 +100,16 @@ public: EulerCoordinateRotation ( const vector& phiThetaPsi, - const bool inDegrees + const bool degrees ); //- Construct from components of rotation vector EulerCoordinateRotation ( - const scalar phiAngle, - const scalar thetaAngle, - const scalar psiAngle, - const bool inDegrees + const scalar phi, + const scalar theta, + const scalar psi, + const bool degrees ); //- Construct from dictionary @@ -144,6 +131,13 @@ public: } + // Static Member Functions + + //- The rotation tensor calculated for the specified Euler angles + //- interpreted as phi/theta/psi (z-x-z order) + static tensor rotation(const vector& angles, bool degrees); + + // Member Functions //- Reset rotation to an identity rotation diff --git a/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.C b/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.C index 2cb61d9d1d..744bbf1877 100644 --- a/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.C +++ b/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -24,8 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "STARCDCoordinateRotation.H" - -#include "mathematicalConstants.H" +#include "unitConversion.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -148,29 +147,26 @@ Foam::symmTensor Foam::STARCDCoordinateRotation::transformVector } -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // -void Foam::STARCDCoordinateRotation::calcTransform +Foam::tensor Foam::STARCDCoordinateRotation::rotation ( - const scalar rotZ, - const scalar rotX, - const scalar rotY, - const bool inDegrees + const vector& angles, + const bool degrees ) { - scalar x = rotX; - scalar y = rotY; - scalar z = rotZ; + scalar z = angles.component(vector::X); // 1. Rotate about Z + scalar x = angles.component(vector::Y); // 2. Rotate about X + scalar y = angles.component(vector::Z); // 3. Rotate about Y - if (inDegrees) + if (degrees) { - x *= constant::mathematical::pi/180.0; - y *= constant::mathematical::pi/180.0; - z *= constant::mathematical::pi/180.0; + x *= degToRad(); + y *= degToRad(); + z *= degToRad(); } - R_ = - ( + return tensor ( cos(y)*cos(z) - sin(x)*sin(y)*sin(z), @@ -184,10 +180,7 @@ void Foam::STARCDCoordinateRotation::calcTransform -cos(x)*sin(y), sin(x), cos(x)*cos(y) - ) - ); - - Rtr_ = R_.T(); + ); } @@ -213,20 +206,12 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation Foam::STARCDCoordinateRotation::STARCDCoordinateRotation ( const vector& rotZrotXrotY, - const bool inDegrees + const bool degrees ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - calcTransform - ( - rotZrotXrotY.component(vector::X), - rotZrotXrotY.component(vector::Y), - rotZrotXrotY.component(vector::Z), - inDegrees - ); -} + R_(rotation(rotZrotXrotY, degrees)), + Rtr_(R_.T()) +{} Foam::STARCDCoordinateRotation::STARCDCoordinateRotation @@ -234,14 +219,12 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation const scalar rotZ, const scalar rotX, const scalar rotY, - const bool inDegrees + const bool degrees ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - calcTransform(rotZ, rotX, rotY, inDegrees); -} + R_(rotation(vector(rotZ, rotX, rotY), degrees)), + Rtr_(R_.T()) +{} Foam::STARCDCoordinateRotation::STARCDCoordinateRotation @@ -249,19 +232,16 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation const dictionary& dict ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - const vector rotation(dict.lookup("rotation")); - - calcTransform + R_ ( - rotation.component(vector::X), - rotation.component(vector::Y), - rotation.component(vector::Z), - dict.lookupOrDefault("degrees", true) - ); -} + rotation + ( + dict.get("rotation"), + dict.lookupOrDefault("degrees", true) + ) + ), + Rtr_(R_.T()) +{} Foam::STARCDCoordinateRotation::STARCDCoordinateRotation diff --git a/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.H b/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.H index 13625e34ba..a6094a8653 100644 --- a/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.H +++ b/src/meshTools/coordinate/rotation/STARCDCoordinateRotation.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -65,7 +65,6 @@ class STARCDCoordinateRotation : public coordinateRotation { - // Private Member Data //- Local-to-Global transformation tensor @@ -75,18 +74,6 @@ class STARCDCoordinateRotation tensor Rtr_; - // Private Member Functions - - //- Calculate transformation tensor - void calcTransform - ( - const scalar rotZ, - const scalar rotX, - const scalar rotY, - const bool inDegrees - ); - - public: //- Runtime type information @@ -105,7 +92,7 @@ public: STARCDCoordinateRotation ( const vector& rotZrotXrotY, - const bool inDegrees + const bool degrees ); //- Construct from components of rotation vector @@ -114,7 +101,7 @@ public: const scalar rotZ, const scalar rotX, const scalar rotY, - const bool inDegrees + const bool degrees ); //- Construct from dictionary @@ -137,6 +124,13 @@ public: } + // Static Member Functions + + //- The rotation tensor calculated for the specified STARCD angles + //- interpreted as rotate-Z, rotate-X, rotate-Y + static tensor rotation(const vector& angles, bool degrees); + + // Member Functions //- Reset rotation to an identity rotation diff --git a/src/meshTools/coordinate/rotation/axesRotation.C b/src/meshTools/coordinate/rotation/axesRotation.C index 1e8973d662..e40e11fe9e 100644 --- a/src/meshTools/coordinate/rotation/axesRotation.C +++ b/src/meshTools/coordinate/rotation/axesRotation.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -47,6 +47,151 @@ namespace Foam } +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::tensor Foam::axesRotation::rotation +( + const vector& axis1, + const vector& axis2, + axisOrder order +) +{ + const scalar magAxis1(mag(axis1)); + scalar magAxis2(mag(axis2)); + + if (magAxis1 < ROOTVSMALL) + { + FatalErrorInFunction + << "Dominant coordinate axis cannot have zero length" + << nl << endl + << abort(FatalError); + } + + const vector ax1(axis1 / magAxis1); // normalise + vector ax2(axis2); + + if (magAxis2 < ROOTVSMALL) + { + // axis2 == Zero : Use best-guess for the second axis. + ax2 = findOrthogonal(axis1); + } + + // Remove colinear component + ax2 -= ((ax1 & ax2) * ax1); + + magAxis2 = mag(ax2); + + if (magAxis2 < SMALL) + { + WarningInFunction + << "axis1, axis2 appear to be co-linear: " + << axis1 << ", " << axis2 << " Revert to guessing axis2" + << nl << endl; + + ax2 = findOrthogonal(axis1); + + // Remove colinear component + ax2 -= ((ax1 & ax2) * ax1); + + magAxis2 = mag(ax2); + + if (magAxis2 < SMALL) + { + FatalErrorInFunction + << "Could not find an appropriate second axis" + << nl << endl + << abort(FatalError); + } + } + + ax2 /= magAxis2; // normalise + + + // The local axes are columns of the rotation matrix + + tensor rotTensor; + + switch (order) + { + case E1_E2: + { + rotTensor.col<0>(ax1); + rotTensor.col<1>(ax2); + rotTensor.col<2>(ax1^ax2); + break; + } + case E2_E3: + { + rotTensor.col<0>(ax1^ax2); + rotTensor.col<1>(ax1); + rotTensor.col<2>(ax2); + break; + } + case E3_E1: + case E3_E1_COMPAT: + { + rotTensor.col<0>(ax2); + rotTensor.col<1>(ax1^ax2); + rotTensor.col<2>(ax1); + break; + } + } + + return rotTensor; +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::axesRotation::read(const dictionary& dict) +{ + vector axis1, axis2; + axisOrder order = E3_E1; + + if + ( + dict.readIfPresent("e1", axis1) + && dict.readIfPresent("e2", axis2) + ) + { + order = E1_E2; + } + else if + ( + dict.readIfPresent("e2", axis1) + && dict.readIfPresent("e3", axis2) + ) + { + order = E2_E3; + } + else if + ( + dict.readIfPresent("e3", axis1) + && dict.readIfPresent("e1", axis2) + ) + { + order = E3_E1; + } + else if + ( + dict.readIfPresent("axis", axis1) + && dict.readIfPresent("direction", axis2) + ) + { + order = E3_E1_COMPAT; + } + else + { + FatalErrorInFunction + << "No entries of the type (e1, e2) or (e2, e3) or (e3, e1) found" + << exit(FatalError); + } + + R_ = rotation(axis1, axis2, order); + Rtr_ = R_.T(); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::axesRotation::axesRotation() @@ -77,11 +222,9 @@ Foam::axesRotation::axesRotation const axisOrder& order ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - setTransform(axis, dir, order); -} + R_(rotation(axis, dir, order)), + Rtr_(R_.T()) +{} Foam::axesRotation::axesRotation @@ -89,37 +232,9 @@ Foam::axesRotation::axesRotation const vector& axis ) : - R_(sphericalTensor::I), - Rtr_(sphericalTensor::I) -{ - direction maxCmpt = 0, dirCmpt = 1; - - scalar maxVal = mag(axis[maxCmpt]); - bool negative = (axis[maxCmpt] < 0); - - for (direction cmpt = 1; cmpt < vector::nComponents; ++cmpt) - { - const scalar val = mag(axis[cmpt]); - - if (maxVal < val) - { - maxVal = val; - maxCmpt = cmpt; - dirCmpt = maxCmpt+1; - negative = (axis[cmpt] < 0); - - if (dirCmpt >= vector::nComponents) - { - dirCmpt = 0; - } - } - } - - vector dir = Zero; - dir.component(dirCmpt) = (negative ? -1 : 1); - - setTransform(axis, dir, E3_E1); -} + R_(rotation(axis, findOrthogonal(axis), E3_E1)), + Rtr_(R_.T()) +{} Foam::axesRotation::axesRotation @@ -130,7 +245,7 @@ Foam::axesRotation::axesRotation R_(sphericalTensor::I), Rtr_(sphericalTensor::I) { - operator=(dict); + read(dict); } @@ -146,62 +261,6 @@ Foam::axesRotation::axesRotation // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -void Foam::axesRotation::setTransform -( - const vector& axis1, - const vector& axis2, - const axisOrder& order -) -{ - const vector a = axis1/mag(axis1); - vector b = axis2; - - b = b - (b & a)*a; - - if (mag(b) < SMALL) - { - FatalErrorInFunction - << "axis1, axis2 appear to be co-linear: " - << axis1 << ", " << axis2 << endl - << abort(FatalError); - } - - b = b/mag(b); - const vector c = a^b; - - // Global->local transformation - switch (order) - { - case E1_E2: - { - Rtr_ = tensor(a, b, c); - break; - } - case E2_E3: - { - Rtr_ = tensor(c, a, b); - break; - } - case E3_E1: - { - Rtr_ = tensor(b, c, a); - break; - } - default: - { - FatalErrorInFunction - << "Unhandled axes specification" << endl - << abort(FatalError); - - break; - } - } - - // Local->global transformation - R_ = Rtr_.T(); -} - - const Foam::tensorField& Foam::axesRotation::Tr() const { NotImplemented; @@ -298,36 +357,7 @@ Foam::symmTensor Foam::axesRotation::transformVector void Foam::axesRotation::operator=(const dictionary& dict) { - vector axis1, axis2; - - if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2)) - { - setTransform(axis1, axis2, E1_E2); - } - else if (dict.readIfPresent("e2", axis1) && dict.readIfPresent("e3", axis2)) - { - setTransform(axis1, axis2, E2_E3); - } - else if (dict.readIfPresent("e3", axis1) && dict.readIfPresent("e1", axis2)) - { - setTransform(axis1, axis2, E3_E1); - } - else if (dict.found("axis") || dict.found("direction")) - { - // Both "axis" and "direction" are required - // If one is missing the appropriate error message will be generated - dict.readEntry("axis", axis1); - dict.readEntry("direction", axis2); - - setTransform(axis1, axis2, E3_E1); - } - else - { - FatalErrorInFunction - << "not entry of the type (e1, e2) or (e2, e3) or (e3, e1) " - << "found " - << exit(FatalError); - } + read(dict); } diff --git a/src/meshTools/coordinate/rotation/axesRotation.H b/src/meshTools/coordinate/rotation/axesRotation.H index 4ac0a04928..3a3583211b 100644 --- a/src/meshTools/coordinate/rotation/axesRotation.H +++ b/src/meshTools/coordinate/rotation/axesRotation.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -73,9 +73,10 @@ public: // Note that these follow the right-hand rule. enum axisOrder { - E1_E2, //!< The axis is X-dominant, the direction is Y-dominant - E2_E3, //!< The axis is Y-dominant, the direction is Z-dominant - E3_E1 //!< The axis is Z-dominant, the direction is X-dominant + E1_E2, //!< The axis1 (dominant) is local X, axis2 is local Y. + E2_E3, //!< The axis1 (dominant) is local Y, axis2 is local Z. + E3_E1, //!< The axis1 (dominant) is local Z, axis2 is local X. + E3_E1_COMPAT, //!< E3_E1 specified as axis/direction. }; @@ -90,6 +91,12 @@ private: tensor Rtr_; + // Protected Member Functions + + //- Read from dictionary + void read(const dictionary& dict); + + public: //- Runtime type information @@ -149,14 +156,19 @@ public: Rtr_ = sphericalTensor::I; } - //- Set the transformation tensors from two axes (axis and direction) - void setTransform + + // Static Member Functions + + //- The rotation tensor calculated from two axes and their order. + // The input axes will be normalised. + static tensor rotation ( const vector& axis1, const vector& axis2, - const axisOrder& order = E3_E1 + axisOrder order = E3_E1 ); + //- Update the rotation for a list of cells virtual void updateCells(const polyMesh&, const labelList&) {} diff --git a/src/meshTools/coordinate/rotation/coordinateRotation.C b/src/meshTools/coordinate/rotation/coordinateRotation.C index 57565dbc45..41892d013d 100644 --- a/src/meshTools/coordinate/rotation/coordinateRotation.C +++ b/src/meshTools/coordinate/rotation/coordinateRotation.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -37,41 +37,65 @@ namespace Foam } -// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::vector Foam::coordinateRotation::findOrthogonal(const vector& axis) +{ + direction maxCmpt = 0; + scalar maxVal = mag(axis[maxCmpt]); + + for (direction cmpt=1; cmpt < vector::nComponents; ++cmpt) + { + const scalar val = mag(axis[cmpt]); + + if (maxVal < val) + { + maxVal = val; + maxCmpt = cmpt; + } + } + + direction cmpt = ((maxCmpt == vector::nComponents-1) ? 0 : (maxCmpt+1)); + + vector dirn(Zero); + dirn.component(cmpt) = ((axis[maxCmpt] < 0) ? -1 : 1); + + return dirn; +} + Foam::symmTensor Foam::coordinateRotation::transformPrincipal ( const tensor& tt, - const vector& st -) const + const vector& v +) { return symmTensor ( - tt.xx()*st.x()*tt.xx() - + tt.xy()*st.y()*tt.xy() - + tt.xz()*st.z()*tt.xz(), + tt.xx()*v.x()*tt.xx() + + tt.xy()*v.y()*tt.xy() + + tt.xz()*v.z()*tt.xz(), - tt.xx()*st.x()*tt.yx() - + tt.xy()*st.y()*tt.yy() - + tt.xz()*st.z()*tt.yz(), + tt.xx()*v.x()*tt.yx() + + tt.xy()*v.y()*tt.yy() + + tt.xz()*v.z()*tt.yz(), - tt.xx()*st.x()*tt.zx() - + tt.xy()*st.y()*tt.zy() - + tt.xz()*st.z()*tt.zz(), + tt.xx()*v.x()*tt.zx() + + tt.xy()*v.y()*tt.zy() + + tt.xz()*v.z()*tt.zz(), - tt.yx()*st.x()*tt.yx() - + tt.yy()*st.y()*tt.yy() - + tt.yz()*st.z()*tt.yz(), + tt.yx()*v.x()*tt.yx() + + tt.yy()*v.y()*tt.yy() + + tt.yz()*v.z()*tt.yz(), - tt.yx()*st.x()*tt.zx() - + tt.yy()*st.y()*tt.zy() - + tt.yz()*st.z()*tt.zz(), + tt.yx()*v.x()*tt.zx() + + tt.yy()*v.y()*tt.zy() + + tt.yz()*v.z()*tt.zz(), - tt.zx()*st.x()*tt.zx() - + tt.zy()*st.y()*tt.zy() - + tt.zz()*st.z()*tt.zz() + tt.zx()*v.x()*tt.zx() + + tt.zy()*v.y()*tt.zy() + + tt.zz()*v.z()*tt.zz() ); - } diff --git a/src/meshTools/coordinate/rotation/coordinateRotation.H b/src/meshTools/coordinate/rotation/coordinateRotation.H index 3a8396b509..4c725756a1 100644 --- a/src/meshTools/coordinate/rotation/coordinateRotation.H +++ b/src/meshTools/coordinate/rotation/coordinateRotation.H @@ -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) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -75,12 +75,14 @@ protected: // Protected member functions //- Transform principal - symmTensor transformPrincipal(const tensor&, const vector&) const; + static symmTensor transformPrincipal(const tensor&, const vector&); + + //- Determine best-guess for an orthogonal axis + static vector findOrthogonal(const vector& axis); public: - //- Runtime type information TypeName("coordinateRotation"); @@ -135,8 +137,7 @@ public: //- Destructor - virtual ~coordinateRotation() - {} + virtual ~coordinateRotation() = default; // Member Functions