Skip to content
Snippets Groups Projects
  • Mark OLESEN's avatar
    c75d5118
    ENH: improve, simplify, rationalize coordinate system handling (issue #863) · c75d5118
    Mark OLESEN authored
    Previously the coordinate system functionality was split between
    coordinateSystem and coordinateRotation. The coordinateRotation stored
    the rotation tensor and handled all tensor transformations.
    
    The functionality has now been revised and consolidated into the
    coordinateSystem classes. The sole purpose of coordinateRotation
    is now just to provide a selectable mechanism of how to define the
    rotation tensor (eg, axis-angle, euler angles, local axes) for user
    input, but after providing the appropriate rotation tensor it has
    no further influence on the transformations.
    
    --
    
    The coordinateSystem class now contains an origin and a base rotation
    tensor directly and various transformation methods.
    
      - The origin represents the "shift" for a local coordinate system.
    
      - The base rotation tensor represents the "tilt" or orientation
        of the local coordinate system in general (eg, for mapping
        positions), but may require position-dependent tensors when
        transforming vectors and tensors.
    
    For some coordinate systems (currently the cylindrical coordinate system),
    the rotation tensor required for rotating a vector or tensor is
    position-dependent.
    
    The new coordinateSystem and its derivates (cartesian, cylindrical,
    indirect) now provide a uniform() method to define if the rotation
    tensor is position dependent/independent.
    
    The coordinateSystem transform and invTransform methods are now
    available in two-parameter forms for obtaining position-dependent
    rotation tensors. Eg,
    
          ... = cs.transform(globalPt, someVector);
    
    In some cases it can be useful to use query uniform() to avoid
    storage of redundant values.
    
          if (cs.uniform())
          {
              vector xx = cs.transform(someVector);
          }
          else
          {
              List<vector> xx = cs.transform(manyPoints, someVector);
          }
    
    Support transform/invTransform for common data types:
       (scalar, vector, sphericalTensor, symmTensor, tensor).
    
    ====================
      Breaking Changes
    ====================
    
    - These changes to coordinate systems and rotations may represent
      a breaking change for existing user coding.
    
    - Relocating the rotation tensor into coordinateSystem itself means
      that the coordinate system 'R()' method now returns the rotation
      directly instead of the coordinateRotation. The method name 'R()'
      was chosen for consistency with other low-level entities (eg,
      quaternion).
    
      The following changes will be needed in coding:
    
          Old:  tensor rot = cs.R().R();
          New:  tensor rot = cs.R();
    
          Old:  cs.R().transform(...);
          New:  cs.transform(...);
    
      Accessing the runTime selectable coordinateRotation
      has moved to the rotation() method:
    
          Old:  Info<< "Rotation input: " << cs.R() << nl;
          New:  Info<< "Rotation input: " << cs.rotation() << nl;
    
    - Naming consistency changes may also cause code to break.
    
          Old:  transformVector()
          New:  transformPrincipal()
    
      The old method name transformTensor() now simply becomes transform().
    
    ====================
      New methods
    ====================
    
    For operations requiring caching of the coordinate rotations, the
    'R()' method can be used with multiple input points:
    
           tensorField rots(cs.R(somePoints));
    
       and later
    
           Foam::transformList(rots, someVectors);
    
    The rotation() method can also be used to change the rotation tensor
    via a new coordinateRotation definition (issue #879).
    
    The new methods transformPoint/invTransformPoint provide
    transformations with an origin offset using Cartesian for both local
    and global points. These can be used to determine the local position
    based on the origin/rotation without interpreting it as a r-theta-z
    value, for example.
    
    ================
      Input format
    ================
    
    - Streamline dictionary input requirements
    
      * The default type is cartesian.
      * The default rotation type is the commonly used axes rotation
        specification (with e1/e2/3), which is assumed if the 'rotation'
        sub-dictionary does not exist.
    
        Example,
    
        Compact specification:
    
            coordinateSystem
            {
                origin  (0 0 0);
                e2      (0 1 0);
                e3      (0.5 0 0.866025);
            }
    
        Full specification (also accepts the longer 'coordinateRotation'
        sub-dictionary name):
    
            coordinateSystem
            {
                type    cartesian;
                origin  (0 0 0);
    
                rotation
                {
                    type    axes;
                    e2      (0 1 0);
                    e3      (0.5 0 0.866025);
                }
            }
    
       This simplifies the input for many cases.
    
    - Additional rotation specification 'none' (an identity rotation):
    
          coordinateSystem
          {
              origin  (0 0 0);
              rotation { type none; }
          }
    
    - Additional rotation specification 'axisAngle', which is similar
      to the -rotate-angle option for transforming points (issue #660).
      For some cases this can be more intuitive.
    
      For example,
    
          rotation
          {
              type    axisAngle;
              axis    (0 1 0);
              angle   30;
          }
      vs.
          rotation
          {
              type    axes;
              e2      (0 1 0);
              e3      (0.5 0 0.866025);
          }
    
    - shorter names (or older longer names) for the coordinate rotation
      specification.
    
         euler         EulerRotation
         starcd        STARCDRotation
         axes          axesRotation
    
    ================
      Coding Style
    ================
    - use Foam::coordSystem namespace for categories of coordinate systems
      (cartesian, cylindrical, indirect). This reduces potential name
      clashes and makes a clearer declaration. Eg,
    
          coordSystem::cartesian csys_;
    
      The older names (eg, cartesianCS, etc) remain available via typedefs.
    
    - added coordinateRotations namespace for better organization and
      reduce potential name clashes.
    c75d5118
    History
    ENH: improve, simplify, rationalize coordinate system handling (issue #863)
    Mark OLESEN authored
    Previously the coordinate system functionality was split between
    coordinateSystem and coordinateRotation. The coordinateRotation stored
    the rotation tensor and handled all tensor transformations.
    
    The functionality has now been revised and consolidated into the
    coordinateSystem classes. The sole purpose of coordinateRotation
    is now just to provide a selectable mechanism of how to define the
    rotation tensor (eg, axis-angle, euler angles, local axes) for user
    input, but after providing the appropriate rotation tensor it has
    no further influence on the transformations.
    
    --
    
    The coordinateSystem class now contains an origin and a base rotation
    tensor directly and various transformation methods.
    
      - The origin represents the "shift" for a local coordinate system.
    
      - The base rotation tensor represents the "tilt" or orientation
        of the local coordinate system in general (eg, for mapping
        positions), but may require position-dependent tensors when
        transforming vectors and tensors.
    
    For some coordinate systems (currently the cylindrical coordinate system),
    the rotation tensor required for rotating a vector or tensor is
    position-dependent.
    
    The new coordinateSystem and its derivates (cartesian, cylindrical,
    indirect) now provide a uniform() method to define if the rotation
    tensor is position dependent/independent.
    
    The coordinateSystem transform and invTransform methods are now
    available in two-parameter forms for obtaining position-dependent
    rotation tensors. Eg,
    
          ... = cs.transform(globalPt, someVector);
    
    In some cases it can be useful to use query uniform() to avoid
    storage of redundant values.
    
          if (cs.uniform())
          {
              vector xx = cs.transform(someVector);
          }
          else
          {
              List<vector> xx = cs.transform(manyPoints, someVector);
          }
    
    Support transform/invTransform for common data types:
       (scalar, vector, sphericalTensor, symmTensor, tensor).
    
    ====================
      Breaking Changes
    ====================
    
    - These changes to coordinate systems and rotations may represent
      a breaking change for existing user coding.
    
    - Relocating the rotation tensor into coordinateSystem itself means
      that the coordinate system 'R()' method now returns the rotation
      directly instead of the coordinateRotation. The method name 'R()'
      was chosen for consistency with other low-level entities (eg,
      quaternion).
    
      The following changes will be needed in coding:
    
          Old:  tensor rot = cs.R().R();
          New:  tensor rot = cs.R();
    
          Old:  cs.R().transform(...);
          New:  cs.transform(...);
    
      Accessing the runTime selectable coordinateRotation
      has moved to the rotation() method:
    
          Old:  Info<< "Rotation input: " << cs.R() << nl;
          New:  Info<< "Rotation input: " << cs.rotation() << nl;
    
    - Naming consistency changes may also cause code to break.
    
          Old:  transformVector()
          New:  transformPrincipal()
    
      The old method name transformTensor() now simply becomes transform().
    
    ====================
      New methods
    ====================
    
    For operations requiring caching of the coordinate rotations, the
    'R()' method can be used with multiple input points:
    
           tensorField rots(cs.R(somePoints));
    
       and later
    
           Foam::transformList(rots, someVectors);
    
    The rotation() method can also be used to change the rotation tensor
    via a new coordinateRotation definition (issue #879).
    
    The new methods transformPoint/invTransformPoint provide
    transformations with an origin offset using Cartesian for both local
    and global points. These can be used to determine the local position
    based on the origin/rotation without interpreting it as a r-theta-z
    value, for example.
    
    ================
      Input format
    ================
    
    - Streamline dictionary input requirements
    
      * The default type is cartesian.
      * The default rotation type is the commonly used axes rotation
        specification (with e1/e2/3), which is assumed if the 'rotation'
        sub-dictionary does not exist.
    
        Example,
    
        Compact specification:
    
            coordinateSystem
            {
                origin  (0 0 0);
                e2      (0 1 0);
                e3      (0.5 0 0.866025);
            }
    
        Full specification (also accepts the longer 'coordinateRotation'
        sub-dictionary name):
    
            coordinateSystem
            {
                type    cartesian;
                origin  (0 0 0);
    
                rotation
                {
                    type    axes;
                    e2      (0 1 0);
                    e3      (0.5 0 0.866025);
                }
            }
    
       This simplifies the input for many cases.
    
    - Additional rotation specification 'none' (an identity rotation):
    
          coordinateSystem
          {
              origin  (0 0 0);
              rotation { type none; }
          }
    
    - Additional rotation specification 'axisAngle', which is similar
      to the -rotate-angle option for transforming points (issue #660).
      For some cases this can be more intuitive.
    
      For example,
    
          rotation
          {
              type    axisAngle;
              axis    (0 1 0);
              angle   30;
          }
      vs.
          rotation
          {
              type    axes;
              e2      (0 1 0);
              e3      (0.5 0 0.866025);
          }
    
    - shorter names (or older longer names) for the coordinate rotation
      specification.
    
         euler         EulerRotation
         starcd        STARCDRotation
         axes          axesRotation
    
    ================
      Coding Style
    ================
    - use Foam::coordSystem namespace for categories of coordinate systems
      (cartesian, cylindrical, indirect). This reduces potential name
      clashes and makes a clearer declaration. Eg,
    
          coordSystem::cartesian csys_;
    
      The older names (eg, cartesianCS, etc) remain available via typedefs.
    
    - added coordinateRotations namespace for better organization and
      reduce potential name clashes.
Test-coordinateSystem.C 5.61 KiB
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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
    Test-coordinateSystem

Description
    Expand coordinate system definitions

\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "Time.H"
#include "coordinateSystems.H"
#include "identityRotation.H"
#include "indirectCS.H"
#include "Fstream.H"
#include "IOstreams.H"
#include "transform.H"

using namespace Foam;


template<class T>
void testTransform(const coordinateSystem& cs, const point& p, const T& val)
{
    Info<< "    " << pTraits<T>::typeName << ": " << val
        << " transform: " << cs.transform(p, val)
        << " invTransform: " << cs.invTransform(p, val) << nl;

    // Info<< " both: " << cs.invTransform(p, cs.transform(p, val)) << nl;
}


void basicTests(const coordinateSystem& cs)
{
    cs.writeEntry(cs.name(), Info);

    if (isA<coordSystem::indirect>(cs))
    {
        Info<< "indirect from:" << nl;
        dynamicCast<const coordSystem::indirect>(cs).cs()
            .writeEntry(cs.name(), Info);
    }


    Info<< "rotation: " << cs.R() << nl;

    List<point> testPoints
    ({
        {1,0,0}, {0,1,0}, {0,0,1}, {1,1,1},
    });


    for (const point& p : testPoints)
    {
        Info<< nl
            << "  test point: " << p
            << " = local point " << cs.transformPoint(p)
            << " = local coord " << cs.localPosition(p) << nl;

        const vector v1(1, 1, 1);
        const tensor t1(tensor::I);
        const tensor t2(1, 2, 3, 4, 5, 6, 7, 8, 9);

        testTransform(cs, p, v1);
        testTransform(cs, p, t1);
        testTransform(cs, p, t2);
    }

    Info<< nl;
}


void doTest(const dictionary& dict)
{
    Info<< dict.dictName() << dict << nl;

    // Could fail?
    const bool throwingIOError = FatalIOError.throwExceptions();
    const bool throwingError = FatalError.throwExceptions();

    try
    {
        auto cs1ptr = coordinateSystem::New(dict, "");
        coordinateSystem& cs1 = *cs1ptr;
        cs1.rename(dict.dictName());

        basicTests(cs1);
    }
    catch (Foam::IOerror& err)
    {
        Info<< "Caught FatalIOError " << err << nl << endl;
    }
    catch (Foam::error& err)
    {
        Info<< "Caught FatalError " << err << nl << endl;
    }
    FatalError.throwExceptions(throwingError);
    FatalIOError.throwExceptions(throwingIOError);
}


void doTest(const objectRegistry& obr, const dictionary& dict)
{
    Info<< dict.dictName() << dict << nl;

    // Could fail?
    const bool throwingIOError = FatalIOError.throwExceptions();
    const bool throwingError = FatalError.throwExceptions();

    try
    {
        auto cs1ptr = coordinateSystem::New(obr, dict, word::null);
        coordinateSystem& cs1 = *cs1ptr;

        basicTests(cs1);
    }
    catch (Foam::IOerror& err)
    {
        Info<< "Caught FatalIOError " << err << nl << endl;
    }
    catch (Foam::error& err)
    {
        Info<< "Caught FatalError " << err << nl << endl;
    }
    FatalError.throwExceptions(throwingError);
    FatalIOError.throwExceptions(throwingIOError);
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//  Main program:

int main(int argc, char *argv[])
{
    argList::noParallel();
    argList::addArgument("dict .. dictN");
    argList args(argc, argv, false, true);

    if (args.found("case"))
    {
        Info<<"using case for tests" << nl;

        #include "createTime.H"

        const coordinateSystems& systems = coordinateSystems::New(runTime);

        Info<< systems.size() << " global systems" << nl;

        for (const coordinateSystem& cs : systems)
        {
            basicTests(cs);
        }

        // systems.write();

        for (label argi=1; argi < args.size(); ++argi)
        {
            const string& dictFile = args[argi];
            IFstream is(dictFile);

            dictionary inputDict(is);

            forAllConstIters(inputDict, iter)
            {
                if (iter().isDict())
                {
                    doTest(runTime, iter().dict());
                }
            }
        }
    }
    else if (args.size() <= 1)
    {
        Info<<"no coordinateSystem dictionaries to expand" << nl;
    }
    else
    {
        for (label argi=1; argi < args.size(); ++argi)
        {
            const string& dictFile = args[argi];
            IFstream is(dictFile);

            dictionary inputDict(is);

            forAllConstIters(inputDict, iter)
            {
                if (iter().isDict())
                {
                    doTest(iter().dict());
                }
            }
        }
    }

    return 0;
}


// ************************************************************************* //