applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     transformPoints
  26
  27 Description
  28     Transforms the mesh points in the polyMesh directory according to the
  29     translate, rotate and scale options.
  30
  31 Usage
  32     Options are:
  33
  34     -translate vector
  35         Translates the points by the given vector,
  36
  37     -rotate (vector vector)
  38         Rotates the points from the first vector to the second,
  39
  40      or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
  41      or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
  42      or -rotateAlongVector (vector and angle)
  43
  44     -scale vector
  45         Scales the points by the given vector.
  46
  47     -cylToCart (vector vector)
  48         Assumes that constant/points is defined in cylindrical coordinates:
  49         (radialPosition tangentialPosition axialPosition) for a coordinate
  50         system with origin: first vec, axis: second vec, direction: third vec
  51         Radial and axial positions should be in [m].
  52         Tangential positions should be in [radians].
  53         Transforms the points to Cartesian positions.
  54
  55     The any or all of the three options may be specified and are processed
  56     in the above order.
  57
  58     With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
  59     it will also read & transform vector & tensor fields.
  60
  61     Note:
  62     yaw (rotation about z)
  63     pitch (rotation about y)
  64     roll (rotation about x)
  65
  66 \*---------------------------------------------------------------------------*/
  67
  68 #include "argList.H"
  69 #include "objectRegistry.H"
  70 #include "foamTime.H"
  71 #include "fvMesh.H"
  72 #include "volFields.H"
  73 #include "surfaceFields.H"
  74 #include "ReadFields.H"
  75 #include "pointFields.H"
  76 #include "transformField.H"
  77 #include "transformGeometricField.H"
  78 #include "IStringStream.H"
  79 #include "RodriguesRotation.H"
  80 #include "cylindricalCS.H"
  81
  82 using namespace Foam;
  83 using namespace Foam::mathematicalConstant;
  84
  85 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  86
  87 template<class GeoField>
  88 void readAndRotateFields
  89 (
  90     PtrList<GeoField>& flds,
  91     const fvMesh& mesh,
  92     const tensor& T,
  93     const IOobjectList& objects
  94 )
  95 {
  96     ReadFields(mesh, objects, flds);
  97     forAll(flds, i)
  98     {
  99         Info<< "Transforming " << flds[i].name() << endl;
 100         dimensionedTensor dimT("t", flds[i].dimensions(), T);
 101         transform(flds[i], dimT, flds[i]);
 102     }
 103 }
 104
 105
 106 void rotateFields(const argList& args, const Time& runTime, const tensor& T)
 107 {
 108 #   include "createNamedMesh.H"
 109
 110     // Read objects in time directory
 111     IOobjectList objects(mesh, runTime.timeName());
 112
 113     // Read vol fields.
 114
 115     PtrList<volScalarField> vsFlds;
 116     readAndRotateFields(vsFlds, mesh, T, objects);
 117
 118     PtrList<volVectorField> vvFlds;
 119     readAndRotateFields(vvFlds, mesh, T, objects);
 120
 121     PtrList<volSphericalTensorField> vstFlds;
 122     readAndRotateFields(vstFlds, mesh, T, objects);
 123
 124     PtrList<volSymmTensorField> vsymtFlds;
 125     readAndRotateFields(vsymtFlds, mesh, T, objects);
 126
 127     PtrList<volTensorField> vtFlds;
 128     readAndRotateFields(vtFlds, mesh, T, objects);
 129
 130     // Read surface fields.
 131
 132     PtrList<surfaceScalarField> ssFlds;
 133     readAndRotateFields(ssFlds, mesh, T, objects);
 134
 135     PtrList<surfaceVectorField> svFlds;
 136     readAndRotateFields(svFlds, mesh, T, objects);
 137
 138     PtrList<surfaceSphericalTensorField> sstFlds;
 139     readAndRotateFields(sstFlds, mesh, T, objects);
 140
 141     PtrList<surfaceSymmTensorField> ssymtFlds;
 142     readAndRotateFields(ssymtFlds, mesh, T, objects);
 143
 144     PtrList<surfaceTensorField> stFlds;
 145     readAndRotateFields(stFlds, mesh, T, objects);
 146
 147     mesh.write();
 148 }
 149
 150
 151 //  Main program:
 152
 153 int main(int argc, char *argv[])
 154 {
 155 #   include "addRegionOption.H"
 156     argList::validOptions.insert("translate", "vector");
 157     argList::validOptions.insert("rotate", "(vector vector)");
 158     argList::validOptions.insert("rotateAlongVector", "(vector angleInDegree)");
 159     argList::validOptions.insert("rollPitchYaw", "(roll pitch yaw)");
 160     argList::validOptions.insert("yawPitchRoll", "(yaw pitch roll)");
 161     argList::validOptions.insert("rotateFields", "");
 162     argList::validOptions.insert("scale", "vector");
 163     argList::validOptions.insert("cylToCart", "(originVec axisVec directionVec)");
 164
 165 #   include "setRootCase.H"
 166 #   include "createTime.H"
 167
 168     word regionName = polyMesh::defaultRegion;
 169     fileName meshDir;
 170
 171     if (args.optionReadIfPresent("region", regionName))
 172     {
 173         meshDir = regionName/polyMesh::meshSubDir;
 174     }
 175     else
 176     {
 177         meshDir = polyMesh::meshSubDir;
 178     }
 179
 180     pointIOField points
 181     (
 182         IOobject
 183         (
 184             "points",
 185             runTime.findInstance(meshDir, "points"),
 186             meshDir,
 187             runTime,
 188             IOobject::MUST_READ,
 189             IOobject::NO_WRITE,
 190             false
 191         )
 192     );
 193
 194
 195     if (args.options().empty())
 196     {
 197         FatalErrorIn(args.executable())
 198             << "No options supplied, please use one or more of "
 199                "-translate, -rotate, -scale, or -cylToCart options."
 200             << exit(FatalError);
 201     }
 202
 203     if (args.optionFound("translate"))
 204     {
 205         vector transVector(args.optionLookup("translate")());
 206
 207         Info<< "Translating points by " << transVector << endl;
 208
 209         points += transVector;
 210     }
 211
 212     if (args.optionFound("rotate"))
 213     {
 214         Pair<vector> n1n2(args.optionLookup("rotate")());
 215         n1n2[0] /= mag(n1n2[0]);
 216         n1n2[1] /= mag(n1n2[1]);
 217         tensor T = rotationTensor(n1n2[0], n1n2[1]);
 218
 219         Info<< "Rotating points by " << T << endl;
 220
 221         points = transform(T, points);
 222
 223         if (args.optionFound("rotateFields"))
 224         {
 225             rotateFields(args, runTime, T);
 226         }
 227     }
 228     else if (args.optionFound("rollPitchYaw"))
 229     {
 230         vector v(args.optionLookup("rollPitchYaw")());
 231
 232         Info<< "Rotating points by" << nl
 233             << "    roll  " << v.x() << nl
 234             << "    pitch " << v.y() << nl
 235             << "    yaw   " << v.z() << endl;
 236
 237
 238         // Convert to radians
 239         v *= pi/180.0;
 240
 241         quaternion R(v.x(), v.y(), v.z());
 242
 243         Info<< "Rotating points by quaternion " << R << endl;
 244         points = transform(R, points);
 245
 246         if (args.optionFound("rotateFields"))
 247         {
 248             rotateFields(args, runTime, R.R());
 249         }
 250     }
 251     else if (args.optionFound("yawPitchRoll"))
 252     {
 253         vector v(args.optionLookup("yawPitchRoll")());
 254
 255         Info<< "Rotating points by" << nl
 256             << "    yaw   " << v.x() << nl
 257             << "    pitch " << v.y() << nl
 258             << "    roll  " << v.z() << endl;
 259
 260
 261         // Convert to radians
 262         v *= pi/180.0;
 263
 264         scalar yaw = v.x();
 265         scalar pitch = v.y();
 266         scalar roll = v.z();
 267
 268         quaternion R = quaternion(vector(0, 0, 1), yaw);
 269         R *= quaternion(vector(0, 1, 0), pitch);
 270         R *= quaternion(vector(1, 0, 0), roll);
 271
 272         Info<< "Rotating points by quaternion " << R << endl;
 273         points = transform(R, points);
 274
 275         if (args.optionFound("rotateFields"))
 276         {
 277             rotateFields(args, runTime, R.R());
 278         }
 279     }
 280     else if (args.optionFound("rotateAlongVector"))
 281     {
 282         vector rotationAxis;
 283         scalar rotationAngle;
 284
 285         args.optionLookup("rotateAlongVector")()
 286             >> rotationAxis
 287             >> rotationAngle;
 288
 289         tensor T = RodriguesRotation(rotationAxis, rotationAngle);
 290
 291         Info << "Rotating points by " << T << endl;
 292
 293         points = transform(T, points);
 294
 295         if (args.options().found("rotateFields"))
 296         {
 297             rotateFields(args, runTime, T);
 298         }
 299     }
 300
 301
 302
 303     if (args.optionFound("scale"))
 304     {
 305         vector scaleVector(args.optionLookup("scale")());
 306
 307         Info<< "Scaling points by " << scaleVector << endl;
 308
 309         points.replace(vector::X, scaleVector.x()*points.component(vector::X));
 310         points.replace(vector::Y, scaleVector.y()*points.component(vector::Y));
 311         points.replace(vector::Z, scaleVector.z()*points.component(vector::Z));
 312     }
 313
 314     if (args.optionFound("cylToCart"))
 315     //HN, 140523
 316     {
 317         vectorField n1n2(args.optionLookup("cylToCart")());
 318         n1n2[1] /= mag(n1n2[1]);
 319         n1n2[2] /= mag(n1n2[2]);
 320
 321         cylindricalCS ccs
 322         (
 323             "ccs",
 324             n1n2[0],
 325             n1n2[1],
 326             n1n2[2],
 327             false // Use radians
 328         );
 329
 330         points = ccs.globalPosition(points);
 331     }
 332
 333     // Set the precision of the points data to 10
 334     IOstream::defaultPrecision(10);
 335
 336     Info << "Writing points into directory " << points.path() << nl << endl;
 337     points.write();
 338
 339     return 0;
 340 }
 341
 342
 343 // ************************************************************************* //