From a61f07451faa306a231a27c9904f0b967920dbe4 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 22 Dec 2014 19:54:58 +0000
Subject: [PATCH] rotorDiskSource: Corrected coordinate transformations Now the
 rotor axis may be specified in other than the z-direction Resolves bug-report
 http://www.openfoam.org/mantisbt/view.php?id=1290

---
 .../derived/rotorDiskSource/rotorDiskSource.C | 42 +++++++++----------
 .../rotorDiskSourceTemplates.C                | 32 +++++++-------
 .../coordinateRotation/axesRotation.C         |  8 ++--
 .../coordinateRotation/axesRotation.H         |  7 ++--
 4 files changed, 46 insertions(+), 43 deletions(-)

diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
index 433407c84..f30f88cfd 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
@@ -69,14 +69,14 @@ namespace Foam
 
 void Foam::fv::rotorDiskSource::checkData()
 {
-    // set inflow type
+    // Set inflow type
     switch (selectionMode())
     {
         case smCellSet:
         case smCellZone:
         case smAll:
         {
-            // set the profile ID for each blade section
+            // Set the profile ID for each blade section
             profiles_.connectBlades(blade_.profileName(), blade_.profileID());
             switch (inletFlow_)
             {
@@ -96,7 +96,7 @@ void Foam::fv::rotorDiskSource::checkData()
                 }
                 case ifLocal:
                 {
-                    // do nothing
+                    // Do nothing
                     break;
                 }
                 default:
@@ -137,7 +137,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
 
     vector n = vector::zero;
 
-    // calculate cell addressing for selected cells
+    // Calculate cell addressing for selected cells
     labelList cellAddr(mesh_.nCells(), -1);
     UIndirectList<label>(cellAddr, cells_) = identity(cells_.size());
     labelList nbrFaceCellAddr(mesh_.nFaces() - nInternalFaces, -1);
@@ -157,10 +157,10 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
         }
     }
 
-    // correct for parallel running
+    // Correct for parallel running
     syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
 
-    // add internal field contributions
+    // Add internal field contributions
     for (label faceI = 0; faceI < nInternalFaces; faceI++)
     {
         const label own = cellAddr[mesh_.faceOwner()[faceI]];
@@ -189,7 +189,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
     }
 
 
-    // add boundary contributions
+    // Add boundary contributions
     forAll(pbm, patchI)
     {
         const polyPatch& pp = pbm[patchI];
@@ -262,7 +262,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
 
 void Foam::fv::rotorDiskSource::createCoordinateSystem()
 {
-    // construct the local rotor co-prdinate system
+    // Construct the local rotor co-prdinate system
     vector origin(vector::zero);
     vector axis(vector::zero);
     vector refDir(vector::zero);
@@ -274,7 +274,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
     {
         case gmAuto:
         {
-            // determine rotation origin (cell volume weighted)
+            // Determine rotation origin (cell volume weighted)
             scalar sumV = 0.0;
             const scalarField& V = mesh_.V();
             const vectorField& C = mesh_.C();
@@ -288,7 +288,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
             reduce(sumV, sumOp<scalar>());
             origin /= sumV;
 
-            // determine first radial vector
+            // Determine first radial vector
             vector dx1(vector::zero);
             scalar magR = -GREAT;
             forAll(cells_, i)
@@ -304,7 +304,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
             reduce(dx1, maxMagSqrOp<vector>());
             magR = mag(dx1);
 
-            // determine second radial vector and cross to determine axis
+            // Determine second radial vector and cross to determine axis
             forAll(cells_, i)
             {
                 const label cellI = cells_[i];
@@ -321,7 +321,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
             reduce(axis, maxMagSqrOp<vector>());
             axis /= mag(axis);
 
-            // correct the axis direction using a point above the rotor
+            // Correct the axis direction using a point above the rotor
             {
                 vector pointAbove(coeffs_.lookup("pointAbove"));
                 vector dir = pointAbove - origin;
@@ -345,7 +345,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
                 )
             );
 
-            // set the face areas and apply correction to calculated axis
+            // Set the face areas and apply correction to calculated axis
             // e.g. if cellZone is more than a single layer in thickness
             setFaceArea(axis, true);
 
@@ -405,20 +405,20 @@ void Foam::fv::rotorDiskSource::constructGeometry()
         {
             const label cellI = cells_[i];
 
-            // position in (planar) rotor co-ordinate system
+            // Position in (planar) rotor co-ordinate system
             x_[i] = coordSys_.localPosition(C[cellI]);
 
-            // cache max radius
+            // Cache max radius
             rMax_ = max(rMax_, x_[i].x());
 
-            // swept angle relative to rDir axis [radians] in range 0 -> 2*pi
+            // Swept angle relative to rDir axis [radians] in range 0 -> 2*pi
             scalar psi = x_[i].y();
 
-            // blade flap angle [radians]
+            // Blade flap angle [radians]
             scalar beta =
                 flap_.beta0 - flap_.beta1c*cos(psi) - flap_.beta2s*sin(psi);
 
-            // determine rotation tensor to convert from planar system into the
+            // Determine rotation tensor to convert from planar system into the
             // rotor cone system
             scalar c = cos(beta);
             scalar s = sin(beta);
@@ -601,7 +601,7 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict)
         coeffs_.lookup("fieldNames") >> fieldNames_;
         applied_.setSize(fieldNames_.size(), false);
 
-        // read co-ordinate system/geometry invariant properties
+        // Read co-ordinate system/geometry invariant properties
         scalar rpm(readScalar(coeffs_.lookup("rpm")));
         omega_ = rpm/60.0*mathematical::twoPi;
 
@@ -620,10 +620,10 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict)
         flap_.beta2s = degToRad(flap_.beta2s);
 
 
-        // create co-ordinate system
+        // Create co-ordinate system
         createCoordinateSystem();
 
-        // read co-odinate system dependent properties
+        // Read co-odinate system dependent properties
         checkData();
 
         constructGeometry();
diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
index 722afcdb0..77fb7adc8 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
@@ -43,7 +43,7 @@ void Foam::fv::rotorDiskSource::calculate
 {
     const scalarField& V = mesh_.V();
 
-    // logging info
+    // Logging info
     scalar dragEff = 0.0;
     scalar liftEff = 0.0;
     scalar AOAmin = GREAT;
@@ -57,19 +57,19 @@ void Foam::fv::rotorDiskSource::calculate
 
             const scalar radius = x_[i].x();
 
-            // velocity in local cylindrical reference frame
-            vector Uc = localAxesRotation_->transform(U[cellI], i);
+            // Transform velocity into local cylindrical reference frame
+            vector Uc = localAxesRotation_->invTransform(U[cellI], i);
 
-            // transform from rotor cylindrical into local coning system
+            // Transform velocity into local coning system
             Uc = R_[i] & Uc;
 
-            // set radial component of velocity to zero
+            // Set radial component of velocity to zero
             Uc.x() = 0.0;
 
-            // set blade normal component of velocity
+            // Set blade normal component of velocity
             Uc.y() = radius*omega_ - Uc.y();
 
-            // determine blade data for this radius
+            // Determine blade data for this radius
             // i2 = index of upper radius bound data point in blade list
             scalar twist = 0.0;
             scalar chord = 0.0;
@@ -78,14 +78,14 @@ void Foam::fv::rotorDiskSource::calculate
             scalar invDr = 0.0;
             blade_.interpolate(radius, twist, chord, i1, i2, invDr);
 
-            // flip geometric angle if blade is spinning in reverse (clockwise)
+            // Flip geometric angle if blade is spinning in reverse (clockwise)
             scalar alphaGeom = thetag[i] + twist;
             if (omega_ < 0)
             {
                 alphaGeom = mathematical::pi - alphaGeom;
             }
 
-            // effective angle of attack
+            // Effective angle of attack
             scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
             if (alphaEff > mathematical::pi)
             {
@@ -99,7 +99,7 @@ void Foam::fv::rotorDiskSource::calculate
             AOAmin = min(AOAmin, alphaEff);
             AOAmax = max(AOAmax, alphaEff);
 
-            // determine profile data for this radius and angle of attack
+            // Determine profile data for this radius and angle of attack
             const label profile1 = blade_.profileID()[i1];
             const label profile2 = blade_.profileID()[i2];
 
@@ -114,24 +114,24 @@ void Foam::fv::rotorDiskSource::calculate
             scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
             scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
 
-            // apply tip effect for blade lift
+            // Apply tip effect for blade lift
             scalar tipFactor = neg(radius/rMax_ - tipEffect_);
 
-            // calculate forces perpendicular to blade
+            // Calculate forces perpendicular to blade
             scalar pDyn = 0.5*rho[cellI]*magSqr(Uc);
 
             scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
             vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
 
-            // accumulate forces
+            // Accumulate forces
             dragEff += rhoRef_*localForce.y();
             liftEff += rhoRef_*localForce.z();
 
-            // convert force from local coning system into rotor cylindrical
+            // Transform force from local coning system into rotor cylindrical
             localForce = invR_[i] & localForce;
 
-            // convert force to global cartesian co-ordinate system
-            force[cellI] = localAxesRotation_->invTransform(localForce, i);
+            // Transform force into global Cartesian co-ordinate system
+            force[cellI] = localAxesRotation_->transform(localForce, i);
 
             if (divideVolume)
             {
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C
index c7a59d35c..03479be54 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C
@@ -103,9 +103,10 @@ void Foam::axesRotation::calcTransform
         }
     }
 
-    // the global->local transformation
+    // Global->local transformation
     Rtr_ = Rtr;
-    // the local->global transformation
+
+    // Local->global transformation
     R_ = Rtr.T();
 }
 
@@ -298,7 +299,8 @@ void Foam::axesRotation::operator=(const dictionary& dict)
     }
     else if (dict.found("axis") || dict.found("direction"))
     {
-        // let it bomb if only one of axis/direction is defined
+        // Both "axis" and "direction" are required
+        // If one is missing the appropriate error message will be generated
         order = e3e1;
         dict.lookup("axis") >> axis1;
         dict.lookup("direction") >> axis2;
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H
index 1478e4459..a0bc1b2c1 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H
+++ b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H
@@ -88,6 +88,7 @@ class axesRotation
             const axisOrder& order = e3e1
         );
 
+
 public:
 
     //- Runtime type information
@@ -134,7 +135,7 @@ public:
         //- Update the rotation for a list of cells
         virtual void updateCells(const polyMesh&, const labelList&)
         {
-            // do nothing
+            // Do nothing
         }
 
         //- Return local-to-global transformation tensor
@@ -196,14 +197,14 @@ public:
         ) const;
 
         //- Transform vectorField using transformation tensorField and return
-        // symmetrical tensorField
+        //  symmetric tensorField
         virtual tmp<symmTensorField> transformVector
         (
             const vectorField& st
         ) const;
 
         //- Transform vector using transformation tensor and return
-        // symmetrical tensor
+        //  symmetric tensor
         virtual symmTensor transformVector(const vector& st) const;