transform: Handle codirectional and contradirectional transformation vectors

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=416

applications/utilities/mesh/manipulation/rotateMesh/rotateMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Application
     rotateMesh
 
 Description
-    Rotates the mesh and fields from the direcion n1 to the direction n2.
+    Rotates the mesh and fields from the direction n1 to direction n2.
 
 \*---------------------------------------------------------------------------*/
 
@@ -71,16 +71,16 @@ int main(int argc, char *argv[])
     argList::validArgs.append("n1");
     argList::validArgs.append("n2");
 
-#   include "setRootCase.H"
+    #include "setRootCase.H"
-#   include "createTime.H"
+    #include "createTime.H"
 
-    vector n1 = args.argRead<vector>(1);
+    vector n1(args.argRead<vector>(1));
     n1 /= mag(n1);
 
-    vector n2 = args.argRead<vector>(2);
+    vector n2(args.argRead<vector>(2));
     n2 /= mag(n2);
 
-    tensor T = rotationTensor(n1, n2);
+    tensor T(rotationTensor(n1, n2));
 
     {
         pointIOField points
@@ -109,8 +109,7 @@ int main(int argc, char *argv[])
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
-#   include "createMesh.H"
+    #include "createMesh.H"
-
 
     forAll(timeDirs, timeI)
     {

src/OpenFOAM/primitives/transform/transform.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,10 +50,29 @@ inline tensor rotationTensor
 {
     const scalar s = n1 & n2;
     const vector n3 = n1 ^ n2;
-    return
+    const scalar magSqrN3 = magSqr(n3);
-        s*I
+
-      + (1 - s)*sqr(n3)/(magSqr(n3) + VSMALL)
+    // n1 and n2 define a plane n3
-      + (n2*n1 - n1*n2);
+    if (magSqrN3 > SMALL)
+    {
+        // Return rotational transformation tensor in the n3-plane
+        return
+            s*I
+          + (1 - s)*sqr(n3)/magSqrN3
+          + (n2*n1 - n1*n2);
+    }
+    // n1 and n2 are contradirectional
+    else if (s < 0)
+    {
+        // Return mirror transformation tensor
+        return I + 2*n1*n2;
+    }
+    // n1 and n2 are codirectional
+    else
+    {
+        // Return null transformation tensor
+        return I;
+    }
 }
 
 
@@ -85,7 +104,6 @@ inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
 template<class Cmpt>
 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
 {
-    //return tt & t & tt.T();
     return Tensor<Cmpt>
     (
         (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.xx()
@@ -190,15 +208,16 @@ inline symmTensor transformMask<symmTensor>(const tensor& t)
     return symm(t);
 }
 
+
 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
 //  Is guaranteed to return increasing number but is not correct
 //  angle. Used for sorting angles.  All input vectors need to be normalized.
 //
-// Calculates scalar which increases with angle going from e0 to vec in
+//  Calculates scalar which increases with angle going from e0 to vec in
-// the coordinate system e0, e1, e0^e1
+//  the coordinate system e0, e1, e0^e1
 //
-// Jumps from 2PI -> 0 at -SMALL so parallel vectors with small rounding errors
+//  Jumps from 2*pi -> 0 at -SMALL so hopefully parallel vectors with small
-// should hopefully still get the same quadrant.
+//  rounding errors should still get the same quadrant.
 //
 inline scalar pseudoAngle
 (