transformPoints.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2017-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
Application
    transformPoints
 
Group
    grpMeshManipulationUtilities
 
Description
    Transforms the mesh points in the polyMesh directory according to the
    translate, rotate and scale options.
 
Usage
    Options are:
 
    -time value
        Specify the time to search from and apply the transformation
        (default is latest)
 
    -recentre
        Recentre using the bounding box centre before other operations
 
    -translate vector
        Translate the points by the given vector before rotations
 
    -rotate (vector vector)
        Rotate 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)
 
    -scale scalar|vector
        Scale the points by the given scalar or vector on output.
 
    The any or all of the three options may be specified and are processed
    in the above order.
 
    With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
    it will also read & transform vector & tensor fields.
 
Note
    roll (rotation about x)
    pitch (rotation about y)
    yaw (rotation about z)
 
\*---------------------------------------------------------------------------*/
 
#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "pointFields.H"
#include "ReadFields.H"
#include "regionProperties.H"
#include "transformField.H"
#include "transformGeometricField.H"
#include "axisAngleRotation.H"
#include "EulerCoordinateRotation.H"
 
using namespace Foam;
using namespace Foam::coordinateRotations;
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
template<class GeoField>
void readAndRotateFields
(
    PtrList<GeoField>& flds,
    const fvMesh& mesh,
    const dimensionedTensor& rotT,
    const IOobjectList& objects
)
{
    ReadFields(mesh, objects, flds);
    for (GeoField& fld : flds)
    {
        Info<< "Transforming " << fld.name() << endl;
        transform(fld, rotT, fld);
    }
}
 
 
void rotateFields
(
    const word& regionName,
    const Time& runTime,
    const tensor& rotT
)
{
    // Need dimensionedTensor for geometric fields
    const dimensionedTensor T(rotT);
 
    #include "createRegionMesh.H"
 
    // Read objects in time directory
    IOobjectList objects(mesh, runTime.timeName());
 
    // Read vol fields.
 
    PtrList<volScalarField> vsFlds;
    readAndRotateFields(vsFlds, mesh, T, objects);
 
    PtrList<volVectorField> vvFlds;
    readAndRotateFields(vvFlds, mesh, T, objects);
 
    PtrList<volSphericalTensorField> vstFlds;
    readAndRotateFields(vstFlds, mesh, T, objects);
 
    PtrList<volSymmTensorField> vsymtFlds;
    readAndRotateFields(vsymtFlds, mesh, T, objects);
 
    PtrList<volTensorField> vtFlds;
    readAndRotateFields(vtFlds, mesh, T, objects);
 
    // Read surface fields.
 
    PtrList<surfaceScalarField> ssFlds;
    readAndRotateFields(ssFlds, mesh, T, objects);
 
    PtrList<surfaceVectorField> svFlds;
    readAndRotateFields(svFlds, mesh, T, objects);
 
    PtrList<surfaceSphericalTensorField> sstFlds;
    readAndRotateFields(sstFlds, mesh, T, objects);
 
    PtrList<surfaceSymmTensorField> ssymtFlds;
    readAndRotateFields(ssymtFlds, mesh, T, objects);
 
    PtrList<surfaceTensorField> stFlds;
    readAndRotateFields(stFlds, mesh, T, objects);
 
    mesh.write();
}
 
 
// Retrieve scaling option
// - size 0 : no scaling
// - size 1 : uniform scaling
// - size 3 : non-uniform scaling
List<scalar> getScalingOpt(const word& optName, const argList& args)
{
    // readListIfPresent handles single or multiple values
    // - accept 1 or 3 values
 
    List<scalar> scaling;
    args.readListIfPresent(optName, scaling);
 
    if (scaling.size() == 1)
    {
        // Uniform scaling
    }
    else if (scaling.size() == 3)
    {
        // Non-uniform, but may actually be uniform
        if
        (
            equal(scaling[0], scaling[1])
         && equal(scaling[0], scaling[2])
        )
        {
            scaling.resize(1);
        }
    }
    else if (!scaling.empty())
    {
        FatalError
            << "Incorrect number of components, must be 1 or 3." << nl
            << "    -" << optName << ' ' << args[optName].c_str() << endl
            << exit(FatalError);
    }
 
    if (scaling.size() == 1 && equal(scaling[0], 1))
    {
        // Scale factor 1 == no scaling
        scaling.clear();
    }
 
    // Zero and negative scaling are permitted
 
    return scaling;
}
 
 
void applyScaling(pointField& points, const List<scalar>& scaling)
{
    if (scaling.size() == 1)
    {
        Info<< "Scaling points uniformly by " << scaling[0] << nl;
        points *= scaling[0];
    }
    else if (scaling.size() == 3)
    {
        Info<< "Scaling points by ("
            << scaling[0] << ' '
            << scaling[1] << ' '
            << scaling[2] << ')' << nl;
 
        points.replace(vector::X, scaling[0]*points.component(vector::X));
        points.replace(vector::Y, scaling[1]*points.component(vector::Y));
        points.replace(vector::Z, scaling[2]*points.component(vector::Z));
    }
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transform (translate / rotate / scale) mesh points.\n"
        "Note: roll=rotate about x, pitch=rotate about y, yaw=rotate about z"
    );
    argList::addOption
    (
        "time",
        "time",
        "Specify the time to search from and apply the transformation"
        " (default is latest)"
    );
    argList::addBoolOption
    (
        "recentre",
        "Recentre the bounding box before other operations"
    );
    argList::addOption
    (
        "translate",
        "vector",
        "Translate by specified <vector> before rotations"
    );
    argList::addBoolOption
    (
        "auto-origin",
        "Use bounding box centre as origin for rotations"
    );
    argList::addOption
    (
        "origin",
        "point",
        "Use specified <point> as origin for rotations"
    );
    argList::addOption
    (
        "rotate",
        "(vectorA vectorB)",
        "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)' degrees"
    );
    argList::addOption
    (
        "yawPitchRoll",
        "vector",
        "Rotate by '(yaw pitch roll)' degrees"
    );
    argList::addBoolOption
    (
        "rotateFields",
        "Read and transform vector and tensor fields too"
    );
    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'"
    );
 
    // Compatibility with surfaceTransformPoints
    argList::addOptionCompat("scale", {"write-scale", 0});
 
    #include "addAllRegionOptions.H"
    #include "setRootCase.H"
 
    const bool doRotateFields = args.found("rotateFields");
 
    // Verify that an operation has been specified
    {
        const List<word> operationNames
        {
            "recentre",
            "translate",
            "rotate",
            "rotate-angle",
            "rollPitchYaw",
            "yawPitchRoll",
            "scale"
        };
 
        if (!args.count(operationNames))
        {
            FatalError
                << "No operation supplied, "
                << "use at least one of the following:" << nl
                << "   ";
 
            for (const auto& opName : operationNames)
            {
                FatalError
                    << " -" << opName;
            }
 
            FatalError
                << nl << exit(FatalError);
        }
    }
 
    // ------------------------------------------------------------------------
 
    #include "createTime.H"
 
    if (args.found("time"))
    {
        if (args["time"] == "constant")
        {
            runTime.setTime(instant(0, "constant"), 0);
        }
        else
        {
            const scalar timeValue = args.get<scalar>("time");
            runTime.setTime(instant(timeValue), 0);
        }
    }
 
    // Handle -allRegions, -regions, -region
    #include "getAllRegionOptions.H"
 
    // ------------------------------------------------------------------------
 
    forAll(regionNames, regioni)
    {
        const word& regionName = regionNames[regioni];
        const word& regionDir =
        (
            regionName == polyMesh::defaultRegion ? word::null : regionName
        );
        const fileName meshDir = regionDir/polyMesh::meshSubDir;
 
        if (regionNames.size() > 1)
        {
            Info<< "region=" << regionName << nl;
        }
 
        pointIOField points
        (
            IOobject
            (
                "points",
                runTime.findInstance(meshDir, "points"),
                meshDir,
                runTime,
                IOobject::MUST_READ,
                IOobject::NO_WRITE,
                false
            )
        );
 
 
        // Begin operations
 
        vector v;
        if (args.found("recentre"))
        {
            v = boundBox(points).centre();
            Info<< "Adjust centre " << v << " -> (0 0 0)" << endl;
            points -= v;
        }
 
        if (args.readIfPresent("translate", v))
        {
            Info<< "Translating points by " << v << endl;
            points += v;
        }
 
        vector origin;
        bool useOrigin = args.readIfPresent("origin", origin);
        if (args.found("auto-origin") && !useOrigin)
        {
            useOrigin = true;
            origin = boundBox(points).centre();
        }
 
        if (useOrigin)
        {
            Info<< "Set origin for rotations to " << origin << endl;
            points -= origin;
        }
 
        if (args.found("rotate"))
        {
            Pair<vector> n1n2
            (
                args.lookup("rotate")()
            );
            n1n2[0].normalise();
            n1n2[1].normalise();
 
            const tensor rot(rotationTensor(n1n2[0], n1n2[1]));
 
            Info<< "Rotating points by " << rot << endl;
            points = transform(rot, points);
 
            if (doRotateFields)
            {
                rotateFields(regionName, runTime, rot);
            }
        }
        else if (args.found("rotate-angle"))
        {
            const Tuple2<vector, scalar> rotAxisAngle
            (
                args.lookup("rotate-angle")()
            );
 
            const vector& axis = rotAxisAngle.first();
            const scalar angle = rotAxisAngle.second();
 
            Info<< "Rotating points " << nl
                << "    about " << axis << nl
                << "    angle " << angle << nl;
 
            const tensor rot(axisAngle::rotation(axis, angle, true));
 
            Info<< "Rotating points by " << rot << endl;
            points = transform(rot, points);
 
            if (doRotateFields)
            {
                rotateFields(regionName, runTime, rot);
            }
        }
        else if (args.readIfPresent("rollPitchYaw", v))
        {
            Info<< "Rotating points by" << nl
                << "    roll  " << v.x() << nl
                << "    pitch " << v.y() << nl
                << "    yaw   " << v.z() << nl;
 
            const tensor rot(euler::rotation(euler::eulerOrder::XYZ, v, true));
 
            Info<< "Rotating points by " << rot << endl;
            points = transform(rot, points);
 
            if (doRotateFields)
            {
                rotateFields(regionName, runTime, rot);
            }
        }
        else if (args.readIfPresent("yawPitchRoll", v))
        {
            Info<< "Rotating points by" << nl
                << "    yaw   " << v.x() << nl
                << "    pitch " << v.y() << nl
                << "    roll  " << v.z() << nl;
 
            const tensor rot(euler::rotation(euler::eulerOrder::ZYX, v, true));
 
            Info<< "Rotating points by " << rot << endl;
            points = transform(rot, points);
 
            if (doRotateFields)
            {
                rotateFields(regionName, runTime, rot);
            }
        }
 
        // Output scaling
        applyScaling(points, getScalingOpt("scale", args));
 
        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()));
 
        Info<< "Writing points into directory "
            << runTime.relativePath(points.path()) << nl
            << endl;
        points.write();
    }
 
    Info<< nl << "End" << nl << endl;
 
    return 0;
}
 
 
// ************************************************************************* //