Merge branch 'master' of ssh://opencfd:8007/home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H

@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(rho, U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPIMPLE() <= 1)
 {

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H

@@ -4,12 +4,12 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        fvOptions(rho, U)
     );
 
     UEqn().relax();
 
-    mrfZones.addCoriolis(rho, UEqn());
-
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
     // including the spherical part of the resistance in the momentum diagonal
@@ -30,7 +30,7 @@
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & ((UEqn() == fvOptions(rho, U))().H() - gradp);
+            U = trTU() & (UEqn().H() - gradp);
         }
         U.correctBoundaryConditions();
 
@@ -42,7 +42,7 @@
 
         fvOptions.constrain(UEqn());
 
-        solve(UEqn() == -fvc::grad(p) + fvOptions(rho, U));
+        solve(UEqn() == -fvc::grad(p));
 
         fvOptions.correct(U);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createZones.H

@@ -1,6 +1,3 @@
-    IOMRFZoneList mrfZones(mesh);
-    mrfZones.correctBoundaryVelocity(U);
-
     IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H

@@ -10,11 +10,11 @@
 
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & (UEqn() == fvOptions(rho, U))().H();
+        HbyA = trTU() & UEqn().H();
     }
     else
     {
-        HbyA = trAU()*(UEqn() == fvOptions(rho, U))().H();
+        HbyA = trAU()*UEqn().H();
     }
 
     UEqn.clear();
@@ -27,7 +27,7 @@
         fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
-    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -27,14 +27,13 @@ Application
 Description
     Steady-state solver for turbulent flow of compressible fluids with
     RANS turbulence modelling, implicit or explicit porosity treatment
-    and MRF for HVAC and similar applications.
+    and run-time selectable finite volume sources.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "RASModel.H"
-#include "IOMRFZoneList.H"
 #include "fvIOoptionList.H"
 #include "IOporosityModelList.H"
 #include "simpleControl.H"

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(rho, U))().H();
+HbyA = rAU*UEqn().H();
 
 UEqn.clear();
 

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options

@@ -4,11 +4,13 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lradiationModels

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H

@@ -9,10 +9,14 @@
         fvm::ddt(T)
       + fvm::div(phi, T)
       - fvm::laplacian(kappaEff, T)
+     ==
+        radiation->ST(rhoCpRef, T)
     );
 
     TEqn.relax();
     TEqn.solve();
 
+    radiation->correct();
+
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -49,6 +49,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "pimpleControl.H"
+#include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +60,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createIncompressibleRadiationModel.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H

@@ -17,7 +17,8 @@
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        fvOptions(rho, he)
+        radiation->Sh(thermo)
+      + fvOptions(rho, he)
     );
 
     EEqn.relax();
@@ -29,4 +30,5 @@
     fvOptions.correct(he);
 
     thermo.correct();
+    radiation->correct();
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -40,13 +40,14 @@ Description
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "zeroGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
 #include "solidRegionDiffNo.H"
 #include "solidThermo.H"
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
-
+#include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -37,6 +37,7 @@ Description
 #include "solidThermo.H"
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
+#include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H

@@ -3,7 +3,9 @@
     {
         fvScalarMatrix hEqn
         (
-           - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+            thermo.isotropic()
+            ? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+            : -fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)")
            + fvOptions(rho, h)
         );
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -1,4 +1,5 @@
     // Initialise solid field pointer lists
+    PtrList<coordinateSystem> coordinates(solidRegions.size());
     PtrList<solidThermo> thermos(solidRegions.size());
     PtrList<radiation::radiationModel> radiations(solidRegions.size());
     PtrList<fv::IOoptionList> solidHeatSources(solidRegions.size());
@@ -23,6 +24,16 @@
             new fv::IOoptionList(solidRegions[i])
         );
 
+        if (!thermos[i].isotropic())
+        {
+            Info<< "    Adding coordinateSystems\n" << endl;
+            coordinates.set
+            (
+                i,
+                coordinateSystem::New(solidRegions[i], thermos[i])
+            );
+        }
+
         IOobject betavSolidIO
         (
             "betavSolid",

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -8,10 +8,36 @@
     tmp<volScalarField> tcp = thermo.Cp();
     const volScalarField& cp = tcp();
 
-    tmp<volScalarField> talpha = thermo.alpha();
+    tmp<volSymmTensorField> tAnialpha;
-    const volScalarField& alpha = talpha();
+
-    tmp<volScalarField> tkappa = thermo.kappa();
+    if (!thermo.isotropic())
-    const volScalarField& kappa = tkappa();
+    {
+        tmp<volVectorField> tkappaCp = thermo.Kappa()/cp;
+
+        const coordinateSystem& coodSys = coordinates[i];
+        tAnialpha =
+            tmp<volSymmTensorField>
+            (
+                new volSymmTensorField
+                (
+                    IOobject
+                    (
+                        "Anialpha",
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh,
+                    tkappaCp().dimensions(),
+                    zeroGradientFvPatchVectorField::typeName
+                )
+            );
+
+        volSymmTensorField& Anialpha = tAnialpha();
+        Anialpha.internalField() = coodSys.R().transformVector(tkappaCp());
+        Anialpha.correctBoundaryConditions();
+    }
 
     volScalarField& h = thermo.he();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C

@@ -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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,38 +56,4 @@ Foam::scalar Foam::solidRegionDiffNo
     return DiNum;
 }
 
-
-Foam::scalar Foam::solidRegionDiffNo
-(
-    const fvMesh& mesh,
-    const Time& runTime,
-    const volScalarField& Cprho,
-    const volSymmTensorField& kappadirectional
-)
-{
-    scalar DiNum = 0.0;
-    scalar meanDiNum = 0.0;
-
-    volScalarField kappa(mag(kappadirectional));
-
-    //- Take care: can have fluid domains with 0 cells so do not test for
-    //  zero internal faces.
-    surfaceScalarField kapparhoCpbyDelta
-    (
-        mesh.surfaceInterpolation::deltaCoeffs()
-      * fvc::interpolate(kappa)
-      / fvc::interpolate(Cprho)
-    );
-
-    DiNum = gMax(kapparhoCpbyDelta.internalField())*runTime.deltaT().value();
-
-    meanDiNum = (average(kapparhoCpbyDelta)).value()*runTime.deltaT().value();
-
-    Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
-        << " max: " << DiNum << endl;
-
-    return DiNum;
-}
-
-
 // ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,15 +41,6 @@ namespace Foam
         const volScalarField& Cprho,
         const volScalarField& kappa
     );
-
-    scalar solidRegionDiffNo
-    (
-        const fvMesh& mesh,
-        const Time& runTime,
-        const volScalarField& Cprho,
-        const volSymmTensorField& kappa
-    );
-
 }
 
 #endif

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -4,6 +4,16 @@
 {
     #include "setRegionSolidFields.H"
 
+    tmp<volScalarField> magKappa;
+    if (thermo.isotropic())
+    {
+        magKappa = thermo.kappa();
+    }
+    else
+    {
+        magKappa = mag(thermo.Kappa());
+    }
+
     DiNum = max
     (
         solidRegionDiffNo
@@ -11,8 +21,9 @@
             solidRegions[i],
             runTime,
             rho*cp,
-                kappa
+            magKappa()
         ),
         DiNum
     );
+
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -9,7 +9,11 @@ if (finalIter)
         tmp<fvScalarMatrix> hEqn
         (
             fvm::ddt(betav*rho, h)
-          - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+          - (
+               thermo.isotropic()
+             ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+             : fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)")
+            )
           ==
             fvOptions(rho, h)
         );

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H

@@ -5,12 +5,14 @@
       + fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
       + SRF->Su()
+     ==
+        fvOptions(Urel)
     );
 
     UrelEqn().relax();
 
     fvOptions.constrain(UrelEqn());
 
-    solve(UrelEqn() == -fvc::grad(p) + fvOptions(Urel));
+    solve(UrelEqn() == -fvc::grad(p));
 
     fvOptions.correct(Urel);

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H

@@ -1,6 +1,6 @@
 volScalarField rAUrel(1.0/UrelEqn().A());
 volVectorField HbyA("HbyA", Urel);
-HbyA = rAUrel*(UrelEqn() == fvOptions(Urel))().H();
+HbyA = rAUrel*UrelEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H

@@ -5,6 +5,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
+  ==
+    fvOptions(U)
 );
 
 UEqn().relax();

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -1,5 +1,5 @@
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H

@@ -3,6 +3,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
 );
 
 
@@ -12,7 +14,7 @@ fvOptions.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p_gh) + fvOptions(U));
+    solve(UEqn() == -fvc::grad(p_gh));
 
     fvOptions.correct(U);
 }

applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H

@@ -2,7 +2,7 @@ volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -18,6 +18,7 @@ surfaceScalarField phiHbyA
 
 adjustPhi(phiHbyA, U, p_gh);
 
+fvOptions.relativeFlux(phiHbyA);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -36,7 +36,6 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "IOporosityModelList.H"
 #include "fvIOoptionList.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"

applications/solvers/multiphase/multiphaseEulerFoam/Make/options

@@ -1,4 +1,4 @@
-EXE_INC = -g \
+EXE_INC = \
     -IphaseModel/lnInclude \
     -ImultiphaseSystem/lnInclude \
     -ImultiphaseFixedFluxPressure \

applications/test/PatchTools/Make/files

@@ -0,0 +1,3 @@
+Test-PatchTools.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-PatchTools

applications/test/PatchTools/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume

applications/test/PatchTools/Test-PatchTools.C

@@ -0,0 +1,297 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     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
+    testPatchTools
+
+Description
+    Test app for PatchTools functionality
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchTools.H"
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "Time.H"
+#include "OBJstream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//template<class PatchType>
+//Foam::tmp<Foam::pointField>
+//areaPointNormals
+//(
+//    const polyMesh& mesh,
+//    const PatchType& p,
+//    const labelList& meshFaces
+//)
+//{
+//    // Assume patch is smaller than the globalData().coupledPatch() (?) so
+//    // loop over patch meshPoints.
+//
+//    const labelList& meshPoints = p.meshPoints();
+//
+//    const globalMeshData& globalData = mesh.globalData();
+//    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+//    const Map<label>& coupledPatchMP = coupledPatch.meshPointMap();
+//    const mapDistribute& map = globalData.globalPointSlavesMap();
+//    const globalIndexAndTransform& transforms =
+//        globalData.globalTransforms();
+//
+//
+//    // 1. Start off with local (area-weighted) normals
+//    //    (note:without calculating pointNormals
+//    //     to avoid them being stored)
+//
+//    tmp<pointField> textrudeN(new pointField(p.nPoints(), vector::zero));
+//    pointField& extrudeN = textrudeN();
+//    {
+//        const faceList& localFaces = p.localFaces();
+//        const vectorField& faceAreas = mesh.faceAreas();
+//
+//        forAll(localFaces, faceI)
+//        {
+//            const face& f = localFaces[faceI];
+//            const vector& n = faceAreas[meshFaces[faceI]];
+//            forAll(f, fp)
+//            {
+//                extrudeN[f[fp]] += n;
+//            }
+//        }
+//    }
+//
+//
+//    // Collect local pointFaces
+//    List<List<point> > pointFaceNormals(map.constructSize());
+//    {
+//        const vectorField& faceAreas = mesh.faceAreas();
+//
+//        forAll(meshPoints, patchPointI)
+//        {
+//            label meshPointI = meshPoints[patchPointI];
+//            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//            if (fnd != coupledPatchMP.end())
+//            {
+//                label coupledPointI = fnd();
+//
+//                List<point>& pNormals = pointFaceNormals[coupledPointI];
+//                const labelList& pFaces = p.pointFaces()[patchPointI];
+//                pNormals.setSize(pFaces.size());
+//                forAll(pFaces, i)
+//                {
+//                    pNormals[i] =  faceAreas[meshFaces[pFaces[i]]];
+//                }
+//            }
+//        }
+//    }
+//
+//    // Pull remote data into local slots
+//    map.distribute
+//    (
+//        transforms,
+//        pointFaceNormals,
+//        listTransform()
+//    );
+//
+//
+//    // Combine normals
+//    const labelListList& slaves = globalData.globalPointSlaves();
+//    const labelListList& transformedSlaves =
+//        globalData.globalPointTransformedSlaves();
+//
+//
+//    pointField coupledPointNormals(map.constructSize(), vector::zero);
+//
+//    forAll(meshPoints, patchPointI)
+//    {
+//        label meshPointI = meshPoints[patchPointI];
+//        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//        if (fnd != coupledPatchMP.end())
+//        {
+//            label coupledPointI = fnd();
+//            const labelList& slaveSlots = slaves[coupledPointI];
+//            const labelList& transformedSlaveSlots =
+//                transformedSlaves[coupledPointI];
+//
+//            label nFaces = slaveSlots.size()+transformedSlaveSlots.size();
+//            if (nFaces > 0)
+//            {
+//                // Combine
+//                point& n = coupledPointNormals[coupledPointI];
+//
+//                n += sum(pointFaceNormals[coupledPointI]);
+//
+//                forAll(slaveSlots, i)
+//                {
+//                    n += sum(pointFaceNormals[slaveSlots[i]]);
+//                }
+//                forAll(transformedSlaveSlots, i)
+//                {
+//                    n += sum(pointFaceNormals[transformedSlaveSlots[i]]);
+//                }
+//
+//                // Put back into slave slots
+//                forAll(slaveSlots, i)
+//                {
+//                    coupledPointNormals[slaveSlots[i]] = n;
+//                }
+//                forAll(transformedSlaveSlots, i)
+//                {
+//                    coupledPointNormals[transformedSlaveSlots[i]] = n;
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Send back
+//    map.reverseDistribute
+//    (
+//        transforms,
+//        coupledPointNormals.size(),
+//        coupledPointNormals,
+//        mapDistribute::transform()
+//    );
+//
+//
+//    // Override patch normals
+//    forAll(meshPoints, patchPointI)
+//    {
+//        label meshPointI = meshPoints[patchPointI];
+//        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//        if (fnd != coupledPatchMP.end())
+//        {
+//            label coupledPointI = fnd();
+//            extrudeN[patchPointI] = coupledPointNormals[coupledPointI];
+//        }
+//    }
+//
+//    extrudeN /= mag(extrudeN)+VSMALL;
+//
+//    return textrudeN;
+//}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+    argList::validArgs.append("patch");
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+#   include "createMesh.H"
+
+    const word patchName = args[1];
+    label patchI = mesh.boundaryMesh().findPatchID(patchName);
+    const polyPatch& pp = mesh.boundaryMesh()[patchI];
+
+    const indirectPrimitivePatch& cpp = mesh.globalData().coupledPatch();
+
+    {
+        OBJstream str(runTime.path()/"edgePatchNormals.obj");
+
+        labelList patchEdges;
+        labelList coupledEdges;
+        PackedBoolList sameEdgeOrientation;
+        PatchTools::matchEdges
+        (
+            pp,
+            cpp,
+            patchEdges,
+            coupledEdges,
+            sameEdgeOrientation
+        );
+
+        const pointField en
+        (
+            PatchTools::edgeNormals
+            (
+                mesh,
+                pp,
+                patchEdges,
+                coupledEdges
+            )
+        );
+
+        forAll(en, patchEdgeI)
+        {
+            const edge& patchE = pp.edges()[patchEdgeI];
+            //str.write(pp.localPoints()[pointI], en[pointI]);
+            const point pt = patchE.centre(pp.localPoints());
+            str.write(linePointRef(pt, pt + 0.1*en[patchEdgeI]));
+        }
+    }
+
+
+    return 0;
+
+
+//    {
+//        OBJstream str(runTime.path()/"unweightedPatchNormals.obj");
+//
+//        const pointField pn
+//        (
+//            PatchTools::pointNormals
+//            (
+//                mesh,
+//                pp,
+//                identity(pp.size())+pp.start()
+//            )
+//        );
+//        forAll(pn, pointI)
+//        {
+//            str.write(linePointRef(pp.localPoints()[pointI], pn[pointI]));
+//        }
+//    }
+//    {
+//        OBJstream str(runTime.path()/"areaWeightedPatchNormals.obj");
+//
+//        const pointField pn
+//        (
+//            areaPointNormals
+//            (
+//                mesh,
+//                pp,
+//                identity(pp.size())+pp.start()
+//            )
+//        );
+//        forAll(pn, pointI)
+//        {
+//            str.write(linePointRef(pp.localPoints()[pointI], pn[pointI]));
+//        }
+//    }
+
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/advanced/collapseEdges/collapseDict

@@ -14,6 +14,11 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// If on, after collapsing check the quality of the mesh. If bad faces are
+// generated then redo the collapsing with stricter filtering.
+controlMeshQuality      on;
+
+
 collapseEdgesCoeffs
 {
     // Edges shorter than this absolute value will be merged
@@ -22,10 +27,6 @@ collapseEdgesCoeffs
     // The maximum angle between two edges that share a point attached to
     // no other edges
     maximumMergeAngle   30;
-
-    // The amount that minimumEdgeLength will be reduced by for each
-    // edge if that edge's collapse generates a poor quality face
-    reductionFactor     0.5;
 }
 
 
@@ -34,10 +35,6 @@ collapseFacesCoeffs
     // The initial face length factor
     initialFaceLengthFactor     0.5;
 
-    // The amount that initialFaceLengthFactor will be reduced by for each
-    // face if its collapse generates a poor quality face
-    reductionFactor                         $initialFaceLengthFactor;
-
     // If the face can't be collapsed to an edge, and it has a span less than
     // the target face length multiplied by this coefficient, collapse it
     // to a point.
@@ -63,12 +60,20 @@ collapseFacesCoeffs
 }
 
 
-meshQualityCoeffs
+controlMeshQualityCoeffs
 {
     // Name of the dictionary that has the mesh quality coefficients used
     // by motionSmoother::checkMesh
     #include                    "meshQualityDict";
 
+    // The amount that minimumEdgeLength will be reduced by for each
+    // edge if that edge's collapse generates a poor quality face
+    edgeReductionFactor         0.5;
+
+    // The amount that initialFaceLengthFactor will be reduced by for each
+    // face if its collapse generates a poor quality face
+    faceReductionFactor         $initialFaceLengthFactor;
+
     // Maximum number of smoothing iterations for the reductionFactors
     maximumSmoothingIterations  2;
 

applications/utilities/miscellaneous/foamHelp/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+wclean libso helpTypes
+wclean
+
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/miscellaneous/foamHelp/helpTypes/doxygenXmlParser/doxygenXmlParser.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,38 +176,32 @@ Foam::doxygenXmlParser::doxygenXmlParser
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::doxygenXmlParser::skipBlock(IFstream& is, const word blockName) const
+void Foam::doxygenXmlParser::skipBlock
+(
+    IFstream& is,
+    const word& blockName
+) const
 {
     // recurse to move forward in 'is' until come across </blockName>
-
     string closeName = "";
 
-    // fast-forward until we reach a '<'
     char c;
+    while (is.good() && (closeName != blockName))
+    {
+        // fast-forward until we reach a '<'
         while (is.get(c) && c  != '<')
         {}
 
         // check to see if this is a closing block
         if (is.get(c) && c  == '/')
         {
+            closeName = "";
+
             while (is.get(c) && c != '>')
             {
                 closeName += c;
             }
-
-        if (closeName == blockName)
-        {
-            // finished reading block
-            return;
         }
-        else
-        {
-            skipBlock(is, blockName);
-        }
-    }
-    else
-    {
-        skipBlock(is, blockName);
     }
 }
 
@@ -215,7 +209,7 @@ void Foam::doxygenXmlParser::skipBlock(IFstream& is, const word blockName) const
 void Foam::doxygenXmlParser::skipForward
 (
     IFstream& is,
-    const word blockName
+    const word& blockName
 ) const
 {
     // recurse to move forward in 'is' until come across <blockName>

applications/utilities/miscellaneous/foamHelp/helpTypes/doxygenXmlParser/doxygenXmlParser.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,10 @@ Class
     Foam::doxygenXmlParser
 
 Description
+    Parser for doxygen XML
 
 SourceFiles
+    doxygenXmlParser.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -65,10 +67,10 @@ public:
     // Member functions
 
         //- Skip past a block
-        void skipBlock(IFstream& is, const word blockName) const;
+        void skipBlock(IFstream& is, const word& blockName) const;
 
         //- Skip forward to block
-        void skipForward(IFstream& is, const word blockName) const;
+        void skipForward(IFstream& is, const word& blockName) const;
 
         //- Return the entry
         template<class Type>

applications/utilities/miscellaneous/foamHelp/helpTypes/helpType/helpType.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,7 +148,12 @@ void Foam::helpType::displayDoc
     {
         FatalErrorIn
         (
-            "void Foam::helpType::displayDoc(const word, const string)"
+            "void Foam::helpType::displayDoc"
+            "("
+                "const word&, "
+                "const string&, "
+                "const bool"
+            ")"
         )
             << "No help for type " << className << " found."
             << "  Valid options include:" << SortableList<word>(parser.toc())

etc/codeTemplates/source/_Template.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_Template.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateApp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/controlDict

@@ -46,6 +46,8 @@ InfoSwitches
 
 OptimisationSwitches
 {
+    // On NFS mounted file system: maximum wait for files to appear/get
+    // updated. Set to 0 on distributed case.
     fileModificationSkew 10;
 
     //- Modification checking:
@@ -1005,6 +1007,22 @@ DimensionSets
         // but not scaled (only supported for dimensionedScalar, etc)
         //writeUnits (kg m s K mol A Cd);
     }
+
+    USCSCoeffs
+    {
+        // Basic units
+        lb  lb  [ 1 0 0 0 0 0 0 ] 1.0;
+        ft  ft  [ 0 1 0 0 0 0 0 ] 1.0;
+        s   s   [ 0 0 1 0 0 0 0 ] 1.0;
+        R   R   [ 0 0 0 1 0 0 0 ] 1.0;
+        mol mol [ 0 0 0 0 1 0 0 ] 1.0;
+        A   A   [ 0 0 0 0 0 1 0 ] 1.0;
+        Cd  Cd  [ 0 0 0 0 0 0 1 ] 1.0;
+
+        // Set of units used for printing. Can be any basic or derived
+        // but not scaled (only supported for dimensionedScalar, etc)
+        //writeUnits (lb ft s R mol A Cd);
+    }
 }
 
 

src/OpenFOAM/db/IOobject/IOobjectReadHeader.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,7 +86,7 @@ bool Foam::IOobject::readHeader(Istream& is)
     }
     else
     {
-        SeriousIOErrorIn("IOobject::readHeader(Istream&)", is)
+        IOWarningIn("IOobject::readHeader(Istream&)", is)
             << "First token could not be read or is not the keyword 'FoamFile'"
             << nl << nl << "Check header is of the form:" << nl << endl;
 

src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,6 +64,7 @@ Foam::label Foam::dictionaryEntry::startLineNumber() const
     }
 }
 
+
 Foam::label Foam::dictionaryEntry::endLineNumber() const
 {
     if (size())

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,16 +27,9 @@ License
 #include "addToMemberFunctionSelectionTable.H"
 #include "IStringStream.H"
 #include "OStringStream.H"
-#include "IOstreams.H"
-#include "stringOps.H"
 #include "dynamicCode.H"
 #include "dynamicCodeContext.H"
-#include "dlLibraryTable.H"
-#include "OSspecific.H"
 #include "Time.H"
-#include "PstreamReduceOps.H"
-
-#include "long.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -149,7 +142,9 @@ Foam::functionEntries::codeStream::getFunction
     // create library if required
     if (!lib)
     {
-        bool create = Pstream::master();
+        bool create =
+            Pstream::master()
+         || (regIOobject::fileModificationSkew <= 0);   // not NFS
 
         if (create)
         {
@@ -196,9 +191,52 @@ Foam::functionEntries::codeStream::getFunction
 
         //- Only block if we're not doing master-only reading. (flag set by
         //  regIOobject::read, IOdictionary constructor)
-        if (!regIOobject::masterOnlyReading)
+        if
+        (
+           !regIOobject::masterOnlyReading
+         && regIOobject::fileModificationSkew > 0
+        )
         {
-            reduce(create, orOp<bool>());
+            //- Since the library has only been compiled on the master the
+            //  other nodes need to pick this library up through NFS
+            //  We do this by just polling a few times using the
+            //  fileModificationSkew.
+
+            off_t mySize = Foam::fileSize(libPath);
+            off_t masterSize = mySize;
+            Pstream::scatter(masterSize);
+
+            if (mySize < masterSize)
+            {
+                Pout<< "Local file " << libPath
+                    << " not of same size (" << mySize
+                    << ") as master ("
+                    << masterSize << "). Waiting for "
+                    << regIOobject::fileModificationSkew
+                    << " seconds." << endl;
+                Foam::sleep(regIOobject::fileModificationSkew);
+
+                // Recheck local size
+                mySize = Foam::fileSize(libPath);
+
+                if (mySize < masterSize)
+                {
+                    FatalIOErrorIn
+                    (
+                        "functionEntries::codeStream::execute(..)",
+                        parentDict
+                    )   << "Cannot read (NFS mounted) library " << nl
+                        << libPath << nl
+                        << "on processor " << Pstream::myProcNo()
+                        << " detected size " << mySize
+                        << " whereas master size is " << masterSize
+                        << " bytes." << nl
+                        << "If your case is not NFS mounted"
+                        << " (so distributed) set fileModificationSkew"
+                        << " to 0"
+                        << exit(FatalIOError);
+                }
+            }
         }
 
         if (isA<IOdictionary>(topDict(parentDict)))

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -221,8 +221,7 @@ public:
     );
 
 
-    //- Return parallel consistent point normals for patches (on boundary faces)
+    //- Return parallel consistent point normals for patches using mesh points.
-    //  using mesh points.
     template
     <
         class Face,
@@ -231,10 +230,27 @@ public:
         class PointType
     >
     static tmp<pointField> pointNormals
+    (
+        const polyMesh&,
+        const PrimitivePatch<Face, FaceList, PointField, PointType>&
+    );
+
+
+    //- Return parallel consistent edge normals for patches using mesh points.
+    //  Supply with patch matching info from matchEdges.
+    template
+    <
+        class Face,
+        template<class> class FaceList,
+        class PointField,
+        class PointType
+    >
+    static tmp<pointField> edgeNormals
     (
         const polyMesh&,
         const PrimitivePatch<Face, FaceList, PointField, PointType>&,
-        const labelList& meshFaces
+        const labelList& patchEdges,
+        const labelList& coupledEdges
     );
 
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,8 +79,7 @@ Foam::tmp<Foam::pointField>
 Foam::PatchTools::pointNormals
 (
     const polyMesh& mesh,
-    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p
-    const labelList& meshFaces
 )
 {
     // Assume patch is smaller than the globalData().coupledPatch() (?) so
@@ -224,4 +223,90 @@ Foam::PatchTools::pointNormals
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+
+Foam::tmp<Foam::pointField>
+Foam::PatchTools::edgeNormals
+(
+    const polyMesh& mesh,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    const labelList& patchEdges,
+    const labelList& coupledEdges
+)
+{
+    // 1. Start off with local normals
+
+    tmp<pointField> tedgeNormals(new pointField(p.nEdges(), vector::zero));
+    pointField& edgeNormals = tedgeNormals();
+    {
+        const labelListList& edgeFaces = p.edgeFaces();
+        const vectorField& faceNormals = p.faceNormals();
+
+        forAll(edgeFaces, edgeI)
+        {
+            const labelList& eFaces = edgeFaces[edgeI];
+            forAll(eFaces, i)
+            {
+                edgeNormals[edgeI] += faceNormals[eFaces[i]];
+            }
+        }
+        edgeNormals /= mag(edgeNormals)+VSMALL;
+    }
+
+
+
+    const globalMeshData& globalData = mesh.globalData();
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+
+
+    // Convert patch-edge data into cpp-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //- Construct with all data in consistent orientation
+    pointField cppEdgeData(map.constructSize(), vector::zero);
+
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+        cppEdgeData[coupledEdgeI] = edgeNormals[patchEdgeI];
+    }
+
+
+    // Synchronise
+    // ~~~~~~~~~~~
+
+    globalData.syncData
+    (
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),
+        map,
+        globalData.globalTransforms(),
+        plusEqOp<point>(),              // add since normalised later on
+        mapDistribute::transform()
+    );
+    cppEdgeData /= mag(cppEdgeData)+VSMALL;
+
+
+    // Back from cpp-edge to patch-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+        edgeNormals[patchEdgeI] = cppEdgeData[coupledEdgeI];
+    }
+
+    return tedgeNormals;
+}
+
+
 // ************************************************************************* //

src/dynamicMesh/polyMeshFilter/polyMeshFilter.C

@@ -427,9 +427,13 @@ Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
             IOobject::NO_WRITE
         )
     ),
+    controlMeshQuality_
+    (
+        dict_.lookupOrDefault<Switch>("controlMeshQuality", false)
+    ),
     collapseEdgesCoeffDict_(dict_.subDict("collapseEdgesCoeffs")),
-    collapseFacesCoeffDict_(dict_.subDict("collapseFacesCoeffs")),
+    collapseFacesCoeffDict_(dict_.subOrEmptyDict("collapseFacesCoeffs")),
-    meshQualityCoeffDict_(dict_.subDict("meshQualityCoeffs")),
+    meshQualityCoeffDict_(dict_.subOrEmptyDict("controlMeshQualityCoeffs")),
     minLen_(readScalar(collapseEdgesCoeffDict_.lookup("minimumEdgeLength"))),
     maxCos_
     (
@@ -443,27 +447,39 @@ Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
     ),
     edgeReductionFactor_
     (
-        readScalar(collapseEdgesCoeffDict_.lookup("reductionFactor"))
+        meshQualityCoeffDict_.lookupOrDefault<scalar>("edgeReductionFactor", -1)
     ),
     maxIterations_
     (
-        readLabel(meshQualityCoeffDict_.lookup("maximumIterations"))
+        meshQualityCoeffDict_.lookupOrAddDefault<label>("maximumIterations", 1)
     ),
     maxSmoothIters_
     (
-        readLabel(meshQualityCoeffDict_.lookup("maximumSmoothingIterations"))
+        meshQualityCoeffDict_.lookupOrAddDefault<label>
+        (
+            "maximumSmoothingIterations",
+            0
+        )
     ),
     initialFaceLengthFactor_
     (
-        readScalar(collapseFacesCoeffDict_.lookup("initialFaceLengthFactor"))
+        collapseFacesCoeffDict_.lookupOrAddDefault<scalar>
+        (
+            "initialFaceLengthFactor",
+            -1
+        )
     ),
     faceReductionFactor_
     (
-        readScalar(collapseFacesCoeffDict_.lookup("reductionFactor"))
+        meshQualityCoeffDict_.lookupOrAddDefault<scalar>
+        (
+            "faceReductionFactor",
+            -1
+        )
     ),
     maxPointErrorCount_
     (
-        readLabel(meshQualityCoeffDict_.lookup("maxPointErrorCount"))
+        meshQualityCoeffDict_.lookupOrAddDefault<label>("maxPointErrorCount", 0)
     ),
     minEdgeLen_(),
     faceFilterFactor_()
@@ -547,23 +563,10 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
                 Map<point> collapsePointToLocation(newMesh.nPoints());
 
                 // Mark points on boundary
-                const labelList boundaryPoint = findBoundaryPoints
+                const labelList boundaryPoint = findBoundaryPoints(newMesh);
-                (
-                    newMesh//,
-    //                    boundaryIOPts
-                );
 
                 edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
 
-                // Per face collapse status:
-                //    -1 : not collapsed
-                //  >= 0 : index of point in face to collapse to
-                List<Map<point> > faceCollapseToPoints
-                (
-                    newMesh.nFaces(),
-                    Map<point>()
-                );
-
                 {
                     // Collapse faces
                     labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
@@ -832,6 +835,8 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
         // Do not allow collapses in regions of error.
         // Updates minEdgeLen, nRelaxedEdges
 
+        if (controlMeshQuality_)
+        {
             PackedBoolList isErrorPoint(newMesh.nPoints());
             nBadFaces = edgeCollapser::checkMeshQuality
             (
@@ -868,6 +873,11 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
                 pointErrorCount
             );
         }
+        else
+        {
+            return -1;
+        }
+    }
 
     return nBadFaces;
 }
@@ -931,11 +941,7 @@ Foam::label Foam::polyMeshFilter::filterEdges
             Map<point> collapsePointToLocation(newMesh.nPoints());
 
             // Mark points on boundary
-            const labelList boundaryPoint = findBoundaryPoints
+            const labelList boundaryPoint = findBoundaryPoints(newMesh);
-            (
-                newMesh//,
-//                    boundaryIOPts
-            );
 
             edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
 
@@ -1059,6 +1065,8 @@ Foam::label Foam::polyMeshFilter::filterEdges
         // Do not allow collapses in regions of error.
         // Updates minEdgeLen, nRelaxedEdges
 
+        if (controlMeshQuality_)
+        {
             PackedBoolList isErrorPoint(newMesh.nPoints());
             nBadFaces = edgeCollapser::checkMeshQuality
             (
@@ -1087,6 +1095,11 @@ Foam::label Foam::polyMeshFilter::filterEdges
                 pointErrorCount
             );
         }
+        else
+        {
+            return -1;
+        }
+    }
 
     return nBadFaces;
 }
@@ -1174,6 +1187,8 @@ Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
         // Do not allow collapses in regions of error.
         // Updates minEdgeLen, nRelaxedEdges
 
+        if (controlMeshQuality_)
+        {
             PackedBoolList isErrorPoint(newMesh.nPoints());
             nBadFaces = edgeCollapser::checkMeshQuality
             (
@@ -1187,6 +1202,7 @@ Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
                 << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
                 << endl;
         }
+    }
 
     return nBadFaces;
 }

src/dynamicMesh/polyMeshFilter/polyMeshFilter.H

@@ -70,14 +70,18 @@ class polyMeshFilter
         //- Dictionary containing the coefficient sub-dictionaries
         const IOdictionary dict_;
 
+        //- After collapsing, check the mesh quality and redo the collapsing
+        //  iteration if there are too many bad faces in the mesh
+        Switch controlMeshQuality_;
+
         //- Coefficients for collapsing edges
         const dictionary& collapseEdgesCoeffDict_;
 
         //- Coefficients for collapsing faces
-        const dictionary& collapseFacesCoeffDict_;
+        dictionary collapseFacesCoeffDict_;
 
         //- Coefficients for controlling the mesh quality
-        const dictionary& meshQualityCoeffDict_;
+        dictionary meshQualityCoeffDict_;
 
         //- Remove edges shorter than this length
         const scalar minLen_;

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -1226,10 +1226,13 @@ Foam::edgeCollapser::edgeCollapser
 )
 :
     mesh_(mesh),
-    guardFraction_(readScalar(dict.lookup("guardFraction"))),
+    guardFraction_
+    (
+        dict.lookupOrDefault<scalar>("guardFraction", 0)
+    ),
     maxCollapseFaceToPointSideLengthCoeff_
     (
-        readScalar(dict.lookup("maxCollapseFaceToPointSideLengthCoeff"))
+        dict.lookupOrDefault<scalar>("maxCollapseFaceToPointSideLengthCoeff", 0)
     ),
     allowEarlyCollapseToPoint_
     (
@@ -1237,7 +1240,7 @@ Foam::edgeCollapser::edgeCollapser
     ),
     allowEarlyCollapseCoeff_
     (
-        readScalar(dict.lookup("allowEarlyCollapseCoeff"))
+        dict.lookupOrDefault<scalar>("allowEarlyCollapseCoeff", 0)
     )
 {}
 

src/engine/enginePiston/enginePiston.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,7 +63,7 @@ Foam::enginePiston::enginePiston
     (
         coordinateSystem::New
         (
-            "coordinateSystem",
+            mesh_,
             dict.subDict("coordinateSystem")
         )
     ),

src/engine/engineValve/engineValve.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,7 +125,7 @@ Foam::engineValve::engineValve
     (
         coordinateSystem::New
         (
-            "coordinateSystem",
+            mesh_,
             dict.subDict("coordinateSystem")
         )
     ),

src/finiteVolume/Make/files

@@ -381,8 +381,8 @@ $(porosity)/porosityModel/porosityModelNew.C
 $(porosity)/porosityModel/porosityModelList.C
 $(porosity)/porosityModel/IOporosityModelList.C
 $(porosity)/DarcyForchheimer/DarcyForchheimer.C
-$(porosity)/powerLaw/powerLaw.C
 $(porosity)/fixedCoeff/fixedCoeff.C
+$(porosity)/powerLaw/powerLaw.C
 
 MRF = $(general)/MRF
 $(MRF)/MRFZone.C

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,66 +52,65 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    coordSys_(coeffs_, mesh),
+    D_(cellZoneIds_.size()),
-    D_("D", dimless/sqr(dimLength), tensor::zero),
+    F_(cellZoneIds_.size()),
-    F_("F", dimless/dimLength, tensor::zero),
     rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
     muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
     nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
 {
-    // local-to-global transformation tensor
-    const tensor& E = coordSys_.R();
 
     dimensionedVector d(coeffs_.lookup("d"));
-    if (D_.dimensions() != d.dimensions())
+    dimensionedVector f(coeffs_.lookup("f"));
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for d: " << d.dimensions()
-            << " should be " << D_.dimensions()
-            << exit(FatalIOError);
-    }
 
     adjustNegativeResistance(d);
-
-    D_.value().xx() = d.value().x();
-    D_.value().yy() = d.value().y();
-    D_.value().zz() = d.value().z();
-    D_.value() = (E & D_ & E.T()).value();
-
-    dimensionedVector f(coeffs_.lookup("f"));
-    if (F_.dimensions() != f.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for f: " << f.dimensions()
-            << " should be " << F_.dimensions()
-            << exit(FatalIOError);
-    }
-
     adjustNegativeResistance(f);
 
+    if (coordSys_.R().uniform())
+    {
+        forAll (cellZoneIds_, zoneI)
+        {
+            D_[zoneI].setSize(1, tensor::zero);
+            F_[zoneI].setSize(1, tensor::zero);
+
+            D_[zoneI][0].xx() = d.value().x();
+            D_[zoneI][0].yy() = d.value().y();
+            D_[zoneI][0].zz() = d.value().z();
+
+            D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]);
+
             // leading 0.5 is from 1/2*rho
-    F_.value().xx() = 0.5*f.value().x();
+            F_[zoneI][0].xx() = 0.5*f.value().x();
-    F_.value().yy() = 0.5*f.value().y();
+            F_[zoneI][0].yy() = 0.5*f.value().y();
-    F_.value().zz() = 0.5*f.value().z();
+            F_[zoneI][0].zz() = 0.5*f.value().z();
-    F_.value() = (E & F_ & E.T()).value();
+
+            F_[zoneI][0] = coordSys_.R().transformTensor(F_[zoneI][0]);
+        }
+
+    }
+    else
+    {
+        forAll(cellZoneIds_, zoneI)
+        {
+            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+            D_[zoneI].setSize(cells.size(), tensor::zero);
+            F_[zoneI].setSize(cells.size(), tensor::zero);
+
+            forAll(cells, i)
+            {
+                D_[zoneI][i].xx() = d.value().x();
+                D_[zoneI][i].yy() = d.value().y();
+                D_[zoneI][i].zz() = d.value().z();
+
+                F_[zoneI][i].xx() = f.value().x();
+                F_[zoneI][i].yy() = f.value().y();
+                F_[zoneI][i].zz() = f.value().z();
+            }
+
+            D_[zoneI] = coordSys_.R().transformTensor(D_[zoneI], cells);
+            F_[zoneI] = coordSys_.R().transformTensor(F_[zoneI], cells);
+        }
+    }
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,7 +53,6 @@ SourceFiles
 #define DarcyForchheimer_H
 
 #include "porosityModel.H"
-#include "coordinateSystem.H"
 #include "dimensionedTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -75,14 +74,12 @@ private:
 
     // Private data
 
-        //- Local co-ordinate system
-        coordinateSystem coordSys_;
 
         //- Darcy coefficient [1/m2]
-        dimensionedTensor D_;
+        List<tensorField> D_;
 
         //- Forchheimer coefficient [1/m]
-        dimensionedTensor F_;
+        List<tensorField> F_;
 
         //- Name of density field
         word rhoName_;

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,18 +36,19 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& dZones = D_[zoneI];
+        const tensorField& fZones = F_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
-
+            const label j = this->fieldIndex(i);
-            const tensor Cd = mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+            const tensor Cd =
+                mu[cellI]*dZones[j] + (rho[cellI]*mag(U[cellI]))*fZones[j];
 
             const scalar isoCd = tr(Cd);
 
@@ -67,16 +68,20 @@ void Foam::porosityModels::DarcyForchheimer::apply
     const vectorField& U
 ) const
 {
-    const tensor& D = D_.value();
-    const tensor& F = F_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& dZones = D_[zoneI];
+        const tensorField& fZones = F_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
+            const label j = this->fieldIndex(i);
+            const tensor D = dZones[j];
+            const tensor F = fZones[j];
+
             AU[cellI] += mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
         }
     }

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,19 +50,18 @@ void Foam::porosityModels::fixedCoeff::apply
     const scalar rho
 ) const
 {
-    const tensor& alpha = alpha_.value();
-    const tensor& beta = beta_.value();
-
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& alphaZones = alpha_[zoneI];
+        const tensorField& betaZones = beta_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
-
+            const label j = fieldIndex(i);
-            const tensor Cd = rho*(alpha + beta*mag(U[cellI]));
+            const tensor Cd = rho*(alphaZones[j] + betaZones[j]*mag(U[cellI]));
-
             const scalar isoCd = tr(Cd);
 
             Udiag[cellI] += V[cellI]*isoCd;
@@ -79,16 +78,21 @@ void Foam::porosityModels::fixedCoeff::apply
     const scalar rho
 ) const
 {
-    const tensor& alpha = alpha_.value();
-    const tensor& beta = beta_.value();
 
     forAll(cellZoneIds_, zoneI)
     {
+        const tensorField& alphaZones = alpha_[zoneI];
+        const tensorField& betaZones = beta_[zoneI];
+
         const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
 
         forAll(cells, i)
         {
             const label cellI = cells[i];
+            const label j = fieldIndex(i);
+            const tensor alpha = alphaZones[j];
+            const tensor beta = betaZones[j];
+
             AU[cellI] += rho*(alpha + beta*mag(U[cellI]));
         }
     }
@@ -107,62 +111,59 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
 )
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
-    coordSys_(coeffs_, mesh),
+    alpha_(cellZoneIds_.size()),
-    alpha_("alpha", dimless/dimTime, tensor::zero),
+    beta_(cellZoneIds_.size())
-    beta_("beta", dimless/dimLength, tensor::zero)
 {
-    // local-to-global transformation tensor
-    const tensor& E = coordSys_.R();
-
     dimensionedVector alpha(coeffs_.lookup("alpha"));
-    if (alpha_.dimensions() != alpha.dimensions())
+    dimensionedVector beta(coeffs_.lookup("beta"));
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::fixedCoeff::fixedCoeff"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for alpha: " << alpha.dimensions()
-            << " should be " << alpha_.dimensions()
-            << exit(FatalIOError);
-    }
 
     adjustNegativeResistance(alpha);
-
-    alpha_.value().xx() = alpha.value().x();
-    alpha_.value().yy() = alpha.value().y();
-    alpha_.value().zz() = alpha.value().z();
-    alpha_.value() = (E & alpha_ & E.T()).value();
-
-    dimensionedVector beta(coeffs_.lookup("beta"));
-    if (beta_.dimensions() != beta.dimensions())
-    {
-        FatalIOErrorIn
-        (
-            "Foam::porosityModels::fixedCoeff::fixedCoeff"
-            "("
-                "const word&, "
-                "const word&, "
-                "const fvMesh&, "
-                "const dictionary&"
-            ")",
-            coeffs_
-        )   << "incorrect dimensions for beta: " << beta.dimensions()
-            << " should be " << beta_.dimensions()
-            << exit(FatalIOError);
-    }
-
     adjustNegativeResistance(beta);
 
-    beta_.value().xx() = beta.value().x();
+    if (coordSys_.R().uniform())
-    beta_.value().yy() = beta.value().y();
+    {
-    beta_.value().zz() = beta.value().z();
+        forAll (cellZoneIds_, zoneI)
-    beta_.value() = (E & beta_ & E.T()).value();
+        {
+            alpha_[zoneI].setSize(1, tensor::zero);
+            beta_[zoneI].setSize(1, tensor::zero);
+
+            alpha_[zoneI][0].xx() = alpha.value().x();
+            alpha_[zoneI][0].yy() = alpha.value().y();
+            alpha_[zoneI][0].zz() = alpha.value().z();
+            alpha_[zoneI][0] = coordSys_.R().transformTensor(alpha_[zoneI][0]);
+
+            beta_[zoneI][0].xx() = beta.value().x();
+            beta_[zoneI][0].yy() = beta.value().y();
+            beta_[zoneI][0].zz() = beta.value().z();
+            beta_[zoneI][0] = coordSys_.R().transformTensor(beta_[zoneI][0]);
+        }
+    }
+    else
+    {
+        forAll(cellZoneIds_, zoneI)
+        {
+            const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+            alpha_[zoneI].setSize(cells.size(), tensor::zero);
+            beta_[zoneI].setSize(cells.size(), tensor::zero);
+
+            forAll(cells, i)
+            {
+                alpha_[zoneI][i].xx() = alpha.value().x();
+                alpha_[zoneI][i].yy() = alpha.value().y();
+                alpha_[zoneI][i].zz() = alpha.value().z();
+
+                beta_[zoneI][i].xx() = beta.value().x();
+                beta_[zoneI][i].yy() = beta.value().y();
+                beta_[zoneI][i].zz() = beta.value().z();
+            }
+
+            alpha_[zoneI] =
+                coordSys_.R().transformTensor(alpha_[zoneI], cells);
+
+            beta_[zoneI] = coordSys_.R().transformTensor(beta_[zoneI], cells);
+        }
+    }
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,6 @@ SourceFiles
 #define fixedCoeff_H
 
 #include "porosityModel.H"
-#include "coordinateSystem.H"
 #include "dimensionedTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -65,14 +64,11 @@ private:
 
     // Private data
 
-        //- Local co-ordinate system
-        coordinateSystem coordSys_;
-
         //- Model alpha coefficient [1/s]
-        dimensionedTensor alpha_;
+        List<tensorField> alpha_;
 
         //- Model beta coefficient [1/m]
-        dimensionedTensor beta_;
+        List<tensorField> beta_;
 
 
     // Private Member Functions

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,6 +66,17 @@ void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
 }
 
 
+Foam::label Foam::porosityModel::fieldIndex(const label i) const
+{
+    label index = 0;
+    if (!coordSys_.R().uniform())
+    {
+        index = i;
+    }
+    return index;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::porosityModel::porosityModel
@@ -83,7 +94,8 @@ Foam::porosityModel::porosityModel
     coeffs_(dict.subDict(modelType + "Coeffs")),
     active_(true),
     zoneName_(cellZoneName),
-    cellZoneIds_()
+    cellZoneIds_(),
+    coordSys_(coordinateSystem::New(mesh, coeffs_))
 {
     if (zoneName_ == word::null)
     {

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,7 @@ SourceFiles
 #include "dictionary.H"
 #include "fvMatricesFwd.H"
 #include "runTimeSelectionTables.H"
+#include "coordinateSystem.H"
 #include "dimensionedVector.H"
 #include "keyType.H"
 
@@ -90,6 +91,9 @@ protected:
         //- Cell zone Ids
         labelList cellZoneIds_;
 
+        //- Local co-ordinate system
+        coordinateSystem coordSys_;
+
 
     // Protected Member Functions
 
@@ -111,6 +115,9 @@ protected:
             volTensorField& AU
         ) const = 0;
 
+        //- Return label index
+        label fieldIndex(const label index) const;
+
 
 public:
 

src/fvOptions/include/createFvOptions.H

@@ -1,2 +1 @@
-Info<< "Creating finite volume options" << endl;
 fv::IOoptionList fvOptions(mesh);

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C

@@ -92,7 +92,7 @@ void Foam::fv::rotorDiskSource::checkData()
                     (
                         readScalar(coeffs_.lookup("inletNormalVelocity"))
                     );
-                    inletVelocity_ = -coordSys_.e3()*UIn;
+                    inletVelocity_ = -coordSys_.R().e3()*UIn;
                     break;
                 }
                 case ifLocal:
@@ -345,9 +345,9 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
         << "    - disk diameter = " << diameter << nl
         << "    - disk area     = " << sumArea << nl
         << "    - origin        = " << coordSys_.origin() << nl
-        << "    - r-axis        = " << coordSys_.e1() << nl
+        << "    - r-axis        = " << coordSys_.R().e1() << nl
-        << "    - psi-axis      = " << coordSys_.e2() << nl
+        << "    - psi-axis      = " << coordSys_.R().e2() << nl
-        << "    - z-axis        = " << coordSys_.e3() << endl;
+        << "    - z-axis        = " << coordSys_.R().e3() << endl;
 }
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C

@@ -59,9 +59,9 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs
     const List<point>& x = rotor_.x();
 
     const vector& origin = rotor_.coordSys().origin();
-    const vector& rollAxis = rotor_.coordSys().e1();
+    const vector& rollAxis = rotor_.coordSys().R().e1();
-    const vector& pitchAxis = rotor_.coordSys().e2();
+    const vector& pitchAxis = rotor_.coordSys().R().e2();
-    const vector& yawAxis = rotor_.coordSys().e3();
+    const vector& yawAxis = rotor_.coordSys().R().e3();
 
     scalar coeff1 = alpha_*sqr(rotor_.omega())*mathematical::pi;
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -471,7 +471,7 @@ public:
             inline scalar rotationalKineticEnergyOfSystem() const;
 
             //- Penetration for fraction [0-1] of the current total mass
-            inline scalar penetration(const scalar& fraction) const;
+            inline scalar penetration(const scalar fraction) const;
 
             //- Mean diameter Dij
             inline scalar Dij(const label i, const label j) const;

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -344,7 +344,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 template<class CloudType>
 inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
 (
-    const scalar& fraction
+    const scalar fraction
 ) const
 {
     if ((fraction < 0) || (fraction > 1))
@@ -353,7 +353,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
         (
             "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
             "("
-                "const scalar&"
+                "const scalar"
             ") const"
         )
             << "fraction should be in the range 0 < fraction < 1"

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,10 +127,10 @@ bool Foam::PilchErdman<CloudType>::update
 
         scalar rho12 = sqrt(rhoc/rho);
 
-        scalar Vd = Urmag*rho12*(B1_*taubBar * B2_*taubBar*taubBar);
+        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
         scalar Vd1 = sqr(1.0 - Vd/Urmag);
         Vd1 = max(Vd1, SMALL);
-        scalar Ds = 2.0*Wec*sigma*Vd1/(Vd1*rhoc*sqr(Urmag));
+        scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
         scalar A = Urmag*rho12/d;
 
         scalar taub = taubBar/A;

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -1464,15 +1464,7 @@ void Foam::autoLayerDriver::getPatchDisplacement
     // Determine pointNormal
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    pointField pointNormals
+    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-    (
-        PatchTools::pointNormals
-        (
-            mesh,
-            pp,
-            pp.addressing()
-        )
-    );
 
 
     // Determine local length scale on patch

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -826,15 +826,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     // Determine pointNormal
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    pointField pointNormals
+    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-    (
-        PatchTools::pointNormals
-        (
-            mesh,
-            pp,
-            pp.addressing()
-        )
-    );
 
     // pointNormals
     if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
@@ -1074,15 +1066,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                         << featureAngle << " degrees." << endl;
 
                     scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
-                    pointField pointNormals
+                    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-                    (
-                        PatchTools::pointNormals
-                        (
-                            mesh,
-                            pp,
-                            identity(pp.size())+pp.start()
-                        )
-                    );
 
                     forAll(meshPoints, i)
                     {

src/meshTools/Make/files

@@ -14,12 +14,13 @@ $(csys)/coordinateSystem.C
 $(csys)/coordinateSystemNew.C
 $(csys)/coordinateSystems.C
 $(csys)/cylindricalCS.C
-$(csys)/sphericalCS.C
+$(csys)/cartesianCS.C
-$(csys)/parabolicCylindricalCS.C
+$(csys)/coordinateRotation/axesRotation.C
-$(csys)/toroidalCS.C
 $(csys)/coordinateRotation/coordinateRotation.C
+$(csys)/coordinateRotation/coordinateRotationNew.C
 $(csys)/coordinateRotation/EulerCoordinateRotation.C
 $(csys)/coordinateRotation/STARCDCoordinateRotation.C
+$(csys)/coordinateRotation/localAxesRotation.C
 
 edgeFaceCirculator/edgeFaceCirculator.C
 
@@ -187,6 +188,9 @@ mappedPatches/mappedPolyPatch/mappedVariableThicknessWallPolyPatch.C
 mappedPatches/mappedPointPatch/mappedPointPatch.C
 mappedPatches/mappedPointPatch/mappedWallPointPatch.C
 
+meshStructure/meshStructure.C
+meshStructure/topoDistanceData.C
+meshStructure/pointTopoDistanceData.C
 
 regionCoupled/patches/regionCoupledPolyPatch/regionCoupledBase.C
 regionCoupled/patches/regionCoupledPolyPatch/regionCoupledPolyPatch.C

src/meshTools/coordinateSystems/cartesianCS.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cartesianCS.H"
+
+#include "one.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cartesianCS, 0);
+    addToRunTimeSelectionTable(coordinateSystem, cartesianCS, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cartesianCS::cartesianCS()
+:
+    coordinateSystem()
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const coordinateSystem& cs
+)
+:
+    coordinateSystem(cs)
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const word& name,
+    const coordinateSystem& cs
+)
+:
+    coordinateSystem(name, cs)
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const word& name,
+    const point& origin,
+    const coordinateRotation& cr
+)
+:
+    coordinateSystem(name, origin, cr)
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const word& name,
+    const point& origin,
+    const vector& axis,
+    const vector& dirn
+)
+:
+    coordinateSystem(name, origin, axis, dirn)
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const word& name,
+    const dictionary& dict
+)
+:
+    coordinateSystem(name, dict)
+{}
+
+
+Foam::cartesianCS::cartesianCS
+(
+    const objectRegistry& obr,
+    const dictionary& dict
+)
+:
+    coordinateSystem(obr, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cartesianCS::~cartesianCS()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::vector Foam::cartesianCS::localToGlobal
+(
+    const vector& local,
+    bool translate
+) const
+{
+    return coordinateSystem::localToGlobal(local, translate);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::cartesianCS::localToGlobal
+(
+    const vectorField& local,
+    bool translate
+) const
+{
+    return coordinateSystem::localToGlobal(local, translate);
+}
+
+
+Foam::vector Foam::cartesianCS::globalToLocal
+(
+    const vector& global,
+    bool translate
+) const
+{
+    return coordinateSystem::globalToLocal(global, translate);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::cartesianCS::globalToLocal
+(
+    const vectorField& global,
+    bool translate
+) const
+{
+    return coordinateSystem::globalToLocal(global, translate);
+}
+
+
+
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/cartesianCS.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,20 +22,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::sphericalCS
+    Foam::cartesianCS
 
 Description
-    Spherical coordinate system
+    Cylindrical coordinate system
 
 SourceFiles
-    sphericalCS.C
+    cartesianCS.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef sphericalCS_H
+#ifndef cartesianCS_H
-#define sphericalCS_H
+#define cartesianCS_H
 
 #include "coordinateSystem.H"
+#include "typeInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,23 +44,18 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class sphericalCS Declaration
+                        Class cartesianCS Declaration
 \*---------------------------------------------------------------------------*/
 
-class sphericalCS
+class cartesianCS
 :
     public coordinateSystem
 {
-    // Private data members
-
-        //- Are angles in degrees? (default = true)
-        bool inDegrees_;
-
-
 protected:
 
     // Protected Member Functions
 
+
         //- Convert from local coordinate system to the global Cartesian system
         //  with optional translation for the origin
         virtual vector localToGlobal(const vector&, bool translate) const;
@@ -88,61 +84,58 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("spherical");
+    TypeName("cartesian");
 
 
     // Constructors
 
         //- Construct null
-        sphericalCS(const bool inDegrees=true);
+        cartesianCS();
 
         //- Construct copy
-        sphericalCS
+        cartesianCS
         (
-            const coordinateSystem&,
+            const coordinateSystem&
-            const bool inDegrees=true
         );
 
         //- Construct copy with a different name
-        sphericalCS
+        cartesianCS
         (
             const word& name,
-            const coordinateSystem&,
+            const coordinateSystem&
-            const bool inDegrees=true
         );
 
         //- Construct from origin and rotation
-        sphericalCS
+        cartesianCS
         (
             const word& name,
             const point& origin,
-            const coordinateRotation&,
+            const coordinateRotation&
-            const bool inDegrees=true
         );
 
         //- Construct from origin and 2 axes
-        sphericalCS
+        cartesianCS
         (
             const word& name,
             const point& origin,
             const vector& axis,
-            const vector& dirn,
+            const vector& dirn
-            const bool inDegrees=true
         );
 
         //- Construct from dictionary
-        sphericalCS(const word& name, const dictionary&);
+        cartesianCS(const word&, const dictionary&);
 
 
-    // Member Functions
+        //- Construct from dictionary and objectRegistry
+        cartesianCS(const objectRegistry&, const dictionary&);
 
-        //- Are angles in degrees?
-        bool  inDegrees() const;
 
-        //- Non-const access to inDegrees
+    //- Destructor
-        bool& inDegrees();
+    virtual ~cartesianCS();
+
 };
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,8 +39,134 @@ namespace Foam
         EulerCoordinateRotation,
         dictionary
     );
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        EulerCoordinateRotation,
+        objectRegistry
+    );
 }
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector Foam::EulerCoordinateRotation::transform(const vector& st) const
+{
+    return (R_ & st);
+}
+
+
+Foam::vector Foam::EulerCoordinateRotation::invTransform
+(
+    const vector& st
+) const
+{
+    return (Rtr_ & st);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::EulerCoordinateRotation::transform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField> Foam::EulerCoordinateRotation:: "
+        "transform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::EulerCoordinateRotation::invTransform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField>  Foam::EulerCoordinateRotation::"
+        "invTransform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+const Foam::tensorField& Foam::EulerCoordinateRotation::Tr() const
+{
+    notImplemented
+    (
+        "const tensorField& EulerCoordinateRotation::Tr() const"
+    );
+    return *reinterpret_cast<const tensorField*>(0);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::EulerCoordinateRotation::transformTensor
+(
+    const tensorField& st
+) const
+{
+     notImplemented
+    (
+        "const tensorField& EulerCoordinateRotation::transformTensor() const"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+
+Foam::tensor Foam::EulerCoordinateRotation::transformTensor
+(
+    const tensor& st
+) const
+{
+    return (R_ & st & Rtr_);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::EulerCoordinateRotation::transformTensor
+(
+    const tensorField& st,
+    const labelList& cellMap
+) const
+{
+    notImplemented
+    (
+        "tmp<Foam::tensorField> EulerCoordinateRotation::transformTensor "
+        " const tensorField& st,"
+        " const labelList& cellMap "
+        ") const"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::EulerCoordinateRotation::
+transformVector
+(
+    const vectorField& st
+) const
+{
+    tmp<symmTensorField> tfld(new symmTensorField(st.size()));
+    symmTensorField& fld = tfld();
+
+    forAll(fld, i)
+    {
+        fld[i] = transformPrincipal(R_, st[i]);
+    }
+    return tfld;
+}
+
+
+Foam::symmTensor Foam::EulerCoordinateRotation::transformVector
+(
+    const vector& st
+) const
+{
+    return transformPrincipal(R_, st);
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::EulerCoordinateRotation::calcTransform
@@ -62,7 +188,7 @@ void Foam::EulerCoordinateRotation::calcTransform
         psi   *= constant::mathematical::pi/180.0;
     }
 
-    tensor::operator=
+    R_ =
     (
         tensor
         (
@@ -79,6 +205,8 @@ void Foam::EulerCoordinateRotation::calcTransform
             cos(theta)
         )
     );
+
+    Rtr_ = R_.T();
 }
 
 
@@ -86,7 +214,8 @@ void Foam::EulerCoordinateRotation::calcTransform
 
 Foam::EulerCoordinateRotation::EulerCoordinateRotation()
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {}
 
 
@@ -96,7 +225,8 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
     const bool inDegrees
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     calcTransform
     (
@@ -116,7 +246,8 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
     const bool inDegrees
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     calcTransform(phiAngle, thetaAngle, psiAngle, inDegrees);
 }
@@ -127,7 +258,8 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
     const dictionary& dict
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     vector rotation(dict.lookup("rotation"));
 
@@ -141,4 +273,32 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
 }
 
 
+Foam::EulerCoordinateRotation::EulerCoordinateRotation
+(
+    const dictionary& dict,
+    const objectRegistry&
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(R_)
+{
+    vector rotation(dict.lookup("rotation"));
+
+    calcTransform
+    (
+        rotation.component(vector::X),
+        rotation.component(vector::Y),
+        rotation.component(vector::Z),
+        dict.lookupOrDefault("degrees", true)
+    );
+}
+
+
+void Foam::EulerCoordinateRotation::write(Ostream& os) const
+{
+     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
+     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
+     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,16 @@ class EulerCoordinateRotation
 :
     public coordinateRotation
 {
+
+    // Private Member Data
+
+        //- Local-to-global transformation tensor
+        tensor R_;
+
+        //- Global-to-Local transformation tensor
+        tensor Rtr_;
+
+
     // Private Member Functions
 
         //- Calculate transformation tensor
@@ -107,6 +117,93 @@ public:
         //- Construct from dictionary
         EulerCoordinateRotation(const dictionary&);
 
+        //- Construct from dictionary and mesh
+        EulerCoordinateRotation(const dictionary&, const objectRegistry&);
+
+
+    // Member Functions
+
+        //- Reset rotation to an identity rotation
+        virtual void clear()
+        {
+            R_ = sphericalTensor::I;
+            Rtr_ = sphericalTensor::I;
+        }
+
+        //- Return local-to-global transformation tensor
+        virtual const tensor& R() const
+        {
+            return R_;
+        }
+
+        //- Return global-to-local transformation tensor
+        virtual const tensor& Rtr() const
+        {
+            return Rtr_;
+        };
+
+        //- Return local Cartesian x-axis
+        virtual const vector e1() const
+        {
+            return R_.x();
+        }
+
+        //- Return local Cartesian y-axis
+        virtual const vector e2() const
+        {
+            return R_.y();
+        }
+
+        //- Return local Cartesian z-axis
+        virtual const vector e3() const
+        {
+            return R_.z();
+        }
+
+        //- Return transformation tensor field
+        virtual const tensorField& Tr() const;
+
+        //- Transform vectorField using transformation tensor field
+        virtual tmp<vectorField> transform(const vectorField& st) const;
+
+        //- Transform vector using transformation tensor
+        virtual vector transform(const vector& st) const;
+
+        //- Inverse transform vectorField using transformation tensor field
+        virtual tmp<vectorField> invTransform(const vectorField& st) const;
+
+        //- Inverse transform vector using transformation tensor
+        virtual vector invTransform(const vector& st) const;
+
+        //- Transform tensor field using transformation tensorField
+        virtual tmp<tensorField> transformTensor(const tensorField& st) const;
+
+        //- Transform tensor using transformation tensorField
+        virtual tensor transformTensor(const tensor& st) const;
+
+        //- Transform tensor sub-field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st,
+            const labelList& cellMap
+        ) const;
+
+        //- Transform vectorField using transformation tensorField and return
+        // symmetrical tensorField
+        virtual tmp<symmTensorField> transformVector
+        (
+            const vectorField& st
+        ) const;
+
+        //- Transform vector using transformation tensor and return
+        // symmetrical tensor
+        virtual symmTensor transformVector(const vector& st) const;
+
+
+    // Write
+
+        //- Write
+        virtual void write(Ostream&) const;
 };
 
 

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,9 +39,134 @@ namespace Foam
         STARCDCoordinateRotation,
         dictionary
     );
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        STARCDCoordinateRotation,
+        objectRegistry
+    );
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector Foam::STARCDCoordinateRotation::transform(const vector& st) const
+{
+    return (R_ & st);
+}
+
+
+Foam::vector Foam::STARCDCoordinateRotation::invTransform
+(
+    const vector& st
+) const
+{
+    return (Rtr_ & st);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::STARCDCoordinateRotation::transform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField> Foam::STARCDCoordinateRotation:: "
+        "transform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::STARCDCoordinateRotation::invTransform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField>  Foam::STARCDCoordinateRotation::"
+        "invTransform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+const Foam::tensorField& Foam::STARCDCoordinateRotation::Tr() const
+{
+    notImplemented
+    (
+        "const tensorField& STARCDCoordinateRotatio::Tr() const"
+    );
+     return *reinterpret_cast<const tensorField*>(0);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::STARCDCoordinateRotation::transformTensor
+(
+    const tensorField& st
+) const
+{
+     notImplemented
+    (
+        "tmp<Foam::tensorField> STARCDCoordinateRotation::transformTensor()"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+
+Foam::tensor Foam::STARCDCoordinateRotation::transformTensor
+(
+    const tensor& st
+) const
+{
+    return (R_ & st & Rtr_);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::STARCDCoordinateRotation::transformTensor
+(
+    const tensorField& st,
+    const labelList& cellMap
+) const
+{
+    notImplemented
+    (
+        "tmp<Foam::tensorField> STARCDCoordinateRotation::transformTensor "
+        " const tensorField& st,"
+        " const labelList& cellMap "
+        ") const"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::STARCDCoordinateRotation::
+transformVector
+(
+    const vectorField& st
+) const
+{
+    tmp<symmTensorField> tfld(new symmTensorField(st.size()));
+    symmTensorField& fld = tfld();
+
+    forAll(fld, i)
+    {
+        fld[i] = transformPrincipal(R_, st[i]);
+    }
+    return tfld;
+}
+
+
+Foam::symmTensor Foam::STARCDCoordinateRotation::transformVector
+(
+    const vector& st
+) const
+{
+    return transformPrincipal(R_, st);
+}
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::STARCDCoordinateRotation::calcTransform
@@ -63,7 +188,7 @@ void Foam::STARCDCoordinateRotation::calcTransform
         z *= constant::mathematical::pi/180.0;
     }
 
-    tensor::operator=
+    R_ =
     (
         tensor
         (
@@ -80,6 +205,8 @@ void Foam::STARCDCoordinateRotation::calcTransform
             cos(x)*cos(y)
         )
     );
+
+    Rtr_ = R_.T();
 }
 
 
@@ -87,7 +214,8 @@ void Foam::STARCDCoordinateRotation::calcTransform
 
 Foam::STARCDCoordinateRotation::STARCDCoordinateRotation()
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {}
 
 
@@ -97,7 +225,8 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
     const bool inDegrees
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     calcTransform
     (
@@ -117,7 +246,8 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
     const bool inDegrees
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     calcTransform(rotZ, rotX, rotY, inDegrees);
 }
@@ -128,7 +258,8 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
     const dictionary& dict
 )
 :
-    coordinateRotation()
+    R_(sphericalTensor::I),
+    Rtr_(R_)
 {
     vector rotation(dict.lookup("rotation"));
 
@@ -141,4 +272,30 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
     );
 }
 
+
+Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
+(
+    const dictionary& dict,
+    const objectRegistry&
+)
+{
+    vector rotation(dict.lookup("rotation"));
+
+    calcTransform
+    (
+        rotation.component(vector::X),
+        rotation.component(vector::Y),
+        rotation.component(vector::Z),
+        dict.lookupOrDefault("degrees", true)
+    );
+}
+
+
+void Foam::STARCDCoordinateRotation::write(Ostream& os) const
+{
+     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
+     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
+     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,6 +62,16 @@ class STARCDCoordinateRotation
 :
     public coordinateRotation
 {
+
+    // Private Member Data
+
+        //- Local-to-Global transformation tensor
+        tensor R_;
+
+        //- Global-to-Local transformation tensor
+        tensor Rtr_;
+
+
     // Private Member Functions
 
         //- Calculate transformation tensor
@@ -104,6 +114,92 @@ public:
         //- Construct from dictionary
         STARCDCoordinateRotation(const dictionary&);
 
+        //- Construct from dictionary and mesh
+        STARCDCoordinateRotation(const dictionary&, const objectRegistry&);
+
+         // Member Functions
+
+        //- Reset rotation to an identity rotation
+        virtual void clear()
+        {
+            R_ = sphericalTensor::I;
+            Rtr_ = sphericalTensor::I;
+        }
+
+        //- Return local-to-global transformation tensor
+        virtual const tensor& R() const
+        {
+            return R_;
+        }
+
+        //- Return global-to-local transformation tensor
+        virtual const tensor& Rtr() const
+        {
+            return Rtr_;
+        };
+
+        //- Return local Cartesian x-axis
+        virtual const vector e1() const
+        {
+            return R_.x();
+        }
+
+        //- Return local Cartesian y-axis
+        virtual const vector e2() const
+        {
+            return R_.y();
+        }
+
+        //- Return local Cartesian z-axis
+        virtual const vector e3() const
+        {
+            return R_.z();
+        }
+
+        //- Return transformation tensor field
+        virtual const tensorField& Tr() const;
+
+        //- Transform vectorField using transformation tensor field
+        virtual tmp<vectorField> transform(const vectorField& st) const;
+
+        //- Transform vector using transformation tensor
+        virtual vector transform(const vector& st) const;
+
+        //- Inverse transform vectorField using transformation tensor field
+        virtual tmp<vectorField> invTransform(const vectorField& st) const;
+
+        //- Inverse transform vector using transformation tensor
+        virtual vector invTransform(const vector& st) const;
+
+        //- Transform tensor field using transformation tensorField
+        virtual tmp<tensorField> transformTensor(const tensorField& st) const;
+
+        //- Transform tensor using transformation tensorField
+        virtual tensor transformTensor(const tensor& st) const;
+
+        //- Transform tensor sub-field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st,
+            const labelList& cellMap
+        ) const;
+
+        //- Transform vectorField using transformation tensorField and return
+        // symmetrical tensorField
+        virtual tmp<symmTensorField> transformVector
+        (
+            const vectorField& st
+        ) const;
+
+        //- Transform vector using transformation tensor and return
+        // symmetrical tensor
+        virtual symmTensor transformVector(const vector& st) const;
+
+
+    // Write
+
+        //- Write
+        virtual void write(Ostream&) const;
 };
 
 

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C

@@ -0,0 +1,324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "axesRotation.H"
+#include "dictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(axesRotation, 0);
+    addToRunTimeSelectionTable(coordinateRotation, axesRotation, dictionary);
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        axesRotation,
+        objectRegistry
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::axesRotation::calcTransform
+(
+    const vector& axis1,
+    const vector& axis2,
+    const axisOrder& order
+)
+{
+    vector a = axis1 / mag(axis1);
+    vector b = axis2;
+
+    // Absorb minor nonorthogonality into axis2
+    b = b - (b & a)*a;
+
+    if (mag(b) < SMALL)
+    {
+        FatalErrorIn("axesRotation::calcTransform()")
+            << "axis1, axis2 appear co-linear: "
+            << axis1 << ", " << axis2 << endl
+            << abort(FatalError);
+    }
+
+    b = b / mag(b);
+    vector c = a ^ b;
+
+    tensor Rtr;
+    switch (order)
+    {
+        case e1e2:
+            Rtr = tensor(a, b, c);
+            break;
+
+        case e2e3:
+            Rtr = tensor(c, a, b);
+            break;
+
+        case e3e1:
+            Rtr = tensor(b, c, a);
+            break;
+
+        default:
+            FatalErrorIn("axesRotation::calcTransform()")
+                << "programmer error" << endl
+                << abort(FatalError);
+
+            Rtr = tensor::zero;
+            break;
+    }
+
+    // the global -> local transformation
+    Rtr_ = Rtr;
+    // the local -> global transformation
+    R_ = Rtr.T();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::axesRotation::axesRotation()
+:
+    R_(sphericalTensor::I),
+    Rtr_(R_)
+{}
+
+
+Foam::axesRotation::axesRotation
+(
+    const vector& axis,
+    const vector& dir
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(R_)
+{
+    calcTransform(axis, dir, e3e1);
+}
+
+
+Foam::axesRotation::axesRotation
+(
+    const dictionary& dict
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(R_)
+{
+    operator=(dict);
+}
+
+
+Foam::axesRotation::axesRotation
+(
+    const dictionary& dict,
+    const objectRegistry& obr
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(R_)
+{
+    operator=(dict);
+}
+
+
+Foam::axesRotation::axesRotation(const tensor& R)
+:
+    R_(R),
+    Rtr_(R_.T())
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::vector Foam::axesRotation::transform(const vector& st) const
+{
+    return (R_ & st);
+}
+
+
+Foam::vector Foam::axesRotation::invTransform(const vector& st) const
+{
+    return (Rtr_ & st);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::axesRotation::transform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField> Foam::axesRotation:: "
+        "transform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::axesRotation::invTransform
+(
+    const vectorField& st
+) const
+{
+    notImplemented
+    (
+        "tmp<vectorField>  Foam::axesRotation::"
+        "invTransform(const vectorField& st) const"
+    );
+    return tmp<vectorField>(NULL);
+}
+
+
+const Foam::tensorField& Foam::axesRotation::Tr() const
+{
+    notImplemented
+    (
+        "const Foam::tensorField& axesRotation::Tr() const"
+    );
+    return *reinterpret_cast<const tensorField*>(0);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::axesRotation::transformTensor
+(
+    const tensorField& st
+) const
+{
+    notImplemented
+    (
+        "const tensorField& axesRotation::transformTensor() const"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+
+Foam::tensor Foam::axesRotation::transformTensor
+(
+    const tensor& st
+) const
+{
+    return (R_ & st & Rtr_);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::axesRotation::transformTensor
+(
+    const tensorField& st,
+    const labelList& cellMap
+) const
+{
+    notImplemented
+    (
+        "tmp<Foam::tensorField> axesRotation::transformTensor "
+        " const tensorField& st,"
+        " const labelList& cellMap "
+        ") const"
+    );
+    return tmp<tensorField>(NULL);
+}
+
+Foam::tmp<Foam::symmTensorField> Foam::axesRotation::transformVector
+(
+    const vectorField& st
+) const
+{
+    tmp<symmTensorField> tfld(new symmTensorField(st.size()));
+    symmTensorField& fld = tfld();
+
+    forAll(fld, i)
+    {
+        fld[i] = transformPrincipal(R_, st[i]);
+    }
+    return tfld;
+}
+
+
+Foam::symmTensor Foam::axesRotation::transformVector
+(
+    const vector& st
+) const
+{
+    return transformPrincipal(R_, st);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::axesRotation::operator=(const dictionary& dict)
+{
+    if (debug)
+    {
+        Pout<< "axesRotation::operator=(const dictionary&) : "
+            << "assign from " << dict << endl;
+    }
+
+    vector axis1, axis2;
+    axisOrder order(e3e1);
+
+    if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2))
+    {
+        order = e1e2;
+    }
+    else if (dict.readIfPresent("e2", axis1)&& dict.readIfPresent("e3", axis2))
+    {
+        order = e2e3;
+    }
+    else if (dict.readIfPresent("e3", axis1)&& dict.readIfPresent("e1", axis2))
+    {
+        order = e3e1;
+    }
+    else if (dict.found("axis") || dict.found("direction"))
+    {
+        // let it bomb if only one of axis/direction is defined
+        order = e3e1;
+        dict.lookup("axis") >> axis1;
+        dict.lookup("direction") >> axis2;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "axesRotation::operator=(const dictionary&) "
+        )   << "not entry of the type (e1, e2) or (e2, e3) or (e3, e1) "
+            << "found "
+            << exit(FatalError);
+    }
+
+    calcTransform(axis1, axis2, order);
+}
+
+
+void Foam::axesRotation::write(Ostream& os) const
+{
+     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
+     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
+     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H

@@ -0,0 +1,229 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    Foam::axesRotation
+
+Description
+    A coordinate rotation specified using global axis
+
+    The rotation is defined by a combination of vectors (e1/e2), (e2/e3)
+    or (e3/e1). Any nonorthogonality will be absorbed into the second vector.
+
+    \verbatim
+        axesRotation
+        {
+            type        axesRotation;
+            e1          (1 0 0);
+            e2          (0 1 0);
+        }
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef axesRotation_H
+#define axesRotation_H
+
+#include "vector.H"
+#include "coordinateRotation.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class axesRotation Declaration
+\*---------------------------------------------------------------------------*/
+
+class axesRotation
+:
+    public coordinateRotation
+{
+    // Private data
+
+        //- Local-to-Global transformation tensor
+        tensor R_;
+
+        //- Global-to-Local transformation tensor
+        tensor Rtr_;
+
+        //- the combination of local axes to be used
+        enum axisOrder
+        {
+            e1e2,
+            e2e3,
+            e3e1
+        };
+
+    // Private Member Functions
+
+        //- Calculate transformation tensor
+        void calcTransform
+        (
+            const vector& axis1,
+            const vector& axis2,
+            const axisOrder& order = e3e1
+        );
+
+public:
+
+    //- Runtime type information
+    TypeName("axesRotation");
+
+    // Constructors
+
+        //- Construct null
+        axesRotation();
+
+        //- Construct from 2 axes
+        axesRotation
+        (
+            const vector& axis,
+            const vector& dir
+        );
+
+        //- Construct from dictionary
+        axesRotation(const dictionary&);
+
+        //- Construct from components
+        axesRotation(const tensor& R);
+
+        //- Construct from dictionary and mesh
+        axesRotation(const dictionary&, const objectRegistry&);
+
+        //- Return clone
+        autoPtr<axesRotation> clone() const
+        {
+            return autoPtr<axesRotation>(new axesRotation(*this));
+        }
+
+
+    //- Destructor
+    virtual ~axesRotation()
+    {}
+
+
+    // Member Functions
+
+        //- Reset rotation to an identity rotation
+        virtual void clear()
+        {
+            R_ = sphericalTensor::I;
+            Rtr_ = sphericalTensor::I;
+        }
+
+        //- Return local-to-global transformation tensor
+        virtual const tensor& R() const
+        {
+            return R_;
+        }
+
+        //- Return global-to-local transformation tensor
+        virtual const tensor& Rtr() const
+        {
+            return Rtr_;
+        }
+
+        //- Return local Cartesian x-axis
+        virtual const vector e1() const
+        {
+            return R_.x();
+        }
+
+        //- Return local Cartesian y-axis
+        virtual const vector e2() const
+        {
+            return R_.y();
+        }
+
+        //- Return local Cartesian z-axis
+        virtual const vector e3() const
+        {
+            return R_.z();
+        }
+
+        //- Return transformation tensor field
+        virtual const tensorField& Tr() const;
+
+        //- Transform vectorField using transformation tensor field
+        virtual tmp<vectorField> transform(const vectorField& st) const;
+
+        //- Transform vector using transformation tensor
+        virtual vector transform(const vector& st) const;
+
+        //- Inverse transform vectorField using transformation tensor field
+        virtual tmp<vectorField> invTransform(const vectorField& st) const;
+
+        //- Inverse transform vector using transformation tensor
+        virtual vector invTransform(const vector& st) const;
+
+        //- Transform tensor field using transformation tensorField
+        virtual tmp<tensorField> transformTensor(const tensorField& st) const;
+
+        //- Transform tensor using transformation tensorField
+        virtual tensor transformTensor(const tensor& st) const;
+
+        //- Transform tensor sub-field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st,
+            const labelList& cellMap
+        ) const;
+
+        //- Transform vectorField using transformation tensorField and return
+        // symmetrical tensorField
+        virtual tmp<symmTensorField> transformVector
+        (
+            const vectorField& st
+        ) const;
+
+        //- Transform vector using transformation tensor and return
+        // symmetrical tensor
+        virtual symmTensor transformVector(const vector& st) const;
+
+
+    // Member Operators
+
+        //- assign from dictionary
+        void operator=(const dictionary&);
+
+
+    // Write
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,197 +33,45 @@ namespace Foam
 {
     defineTypeNameAndDebug(coordinateRotation, 0);
     defineRunTimeSelectionTable(coordinateRotation, dictionary);
-}
+    defineRunTimeSelectionTable(coordinateRotation, objectRegistry);
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::coordinateRotation::calcTransform
-(
-    const vector& axis1,
-    const vector& axis2,
-    const axisOrder& order
-)
-{
-    vector a = axis1 / mag(axis1);
-    vector b = axis2;
-
-    // Absorb minor nonorthogonality into axis2
-    b = b - (b & a)*a;
-
-    if (mag(b) < SMALL)
-    {
-        FatalErrorIn("coordinateRotation::calcTransform()")
-            << "axis1, axis2 appear co-linear: "
-            << axis1 << ", " << axis2 << endl
-            << abort(FatalError);
-    }
-
-    b = b / mag(b);
-    vector c = a ^ b;
-
-    // the global -> local transformation
-    tensor Rtr;
-    switch (order)
-    {
-        case e1e2:
-            Rtr = tensor(a, b, c);
-            break;
-
-        case e2e3:
-            Rtr = tensor(c, a, b);
-            break;
-
-        case e3e1:
-            Rtr = tensor(b, c, a);
-            break;
-
-        default:
-            FatalErrorIn("coordinateRotation::calcTransform()")
-                << "programmer error" << endl
-                << abort(FatalError);
-            // To satisfy compiler warnings
-            Rtr = tensor::zero;
-            break;
-    }
-
-    // the local -> global transformation
-    tensor::operator=( Rtr.T() );
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::coordinateRotation::coordinateRotation()
-:
-    tensor(sphericalTensor::I)
-{}
-
-
-Foam::coordinateRotation::coordinateRotation
-(
-    const vector& axis,
-    const vector& dir
-)
-:
-    tensor(sphericalTensor::I)
-{
-    calcTransform(axis, dir, e3e1);
-}
-
-
-Foam::coordinateRotation::coordinateRotation
-(
-    const dictionary& dict
-)
-:
-    tensor(sphericalTensor::I)
-{
-    operator=(dict);
-}
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
-(
-    const dictionary& dict
-)
-{
-    if (debug)
-    {
-        Pout<< "coordinateRotation::New(const dictionary&) : "
-            << "constructing coordinateRotation"
-            << endl;
-    }
-
-    // default type is self (alias: "axes")
-    word rotType(typeName_());
-    dict.readIfPresent("type", rotType);
-
-    // can (must) construct base class directly
-    if (rotType == typeName_() || rotType == "axes")
-    {
-        return autoPtr<coordinateRotation>(new coordinateRotation(dict));
-    }
-
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(rotType);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalIOErrorIn
-        (
-            "coordinateRotation::New(const dictionary&)",
-            dict
-        )   << "Unknown coordinateRotation type "
-            << rotType << nl << nl
-            << "Valid coordinateRotation types are :" <<  nl
-            << "[default: axes " << typeName_() << "]"
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalIOError);
-    }
-
-    return autoPtr<coordinateRotation>(cstrIter()(dict));
 }
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::coordinateRotation::clear()
+Foam::symmTensor Foam::coordinateRotation::transformPrincipal
-{
-    this->tensor::operator=(sphericalTensor::I);
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-void Foam::coordinateRotation::operator=(const dictionary& rhs)
-{
-    if (debug)
-    {
-        Pout<< "coordinateRotation::operator=(const dictionary&) : "
-            << "assign from " << rhs << endl;
-    }
-
-    // allow as embedded sub-dictionary "coordinateRotation"
-    const dictionary& dict =
 (
-        rhs.found(typeName_())
+    const tensor& tt,
-      ? rhs.subDict(typeName_())
+    const vector& st
-      : rhs
+) const
+{
+    return symmTensor
+    (
+        tt.xx()*st.x()*tt.xx()
+      + tt.xy()*st.y()*tt.xy()
+      + tt.xz()*st.z()*tt.xz(),
+
+        tt.xx()*st.x()*tt.yx()
+      + tt.xy()*st.y()*tt.yy()
+      + tt.xz()*st.z()*tt.yz(),
+
+        tt.xx()*st.x()*tt.zx()
+      + tt.xy()*st.y()*tt.zy()
+      + tt.xz()*st.z()*tt.zz(),
+
+        tt.yx()*st.x()*tt.yx()
+      + tt.yy()*st.y()*tt.yy()
+      + tt.yz()*st.z()*tt.yz(),
+
+        tt.yx()*st.x()*tt.zx()
+      + tt.yy()*st.y()*tt.zy()
+      + tt.yz()*st.z()*tt.zz(),
+
+        tt.zx()*st.x()*tt.zx()
+      + tt.zy()*st.y()*tt.zy()
+      + tt.zz()*st.z()*tt.zz()
     );
 
-    vector axis1, axis2;
-    axisOrder order(e3e1);
-
-    if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2))
-    {
-        order = e1e2;
 }
-    else if (dict.readIfPresent("e2", axis1) && dict.readIfPresent("e3", axis2))
-    {
-        order = e2e3;
-    }
-    else if (dict.readIfPresent("e3", axis1) && dict.readIfPresent("e1", axis2))
-    {
-        order = e3e1;
-    }
-    else if (dict.found("axis") || dict.found("direction"))
-    {
-        // let it bomb if only one of axis/direction is defined
-        order = e3e1;
-        dict.lookup("axis") >> axis1;
-        dict.lookup("direction") >> axis2;
-    }
-    else
-    {
-        // unspecified axes revert to the global system
-        tensor::operator=(sphericalTensor::I);
-        return;
-    }
-
-    calcTransform(axis1, axis2, order);
-}
-
 
 // ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,46 +25,24 @@ Class
     Foam::coordinateRotation
 
 Description
-    A coordinate rotation specified per local axes and the base class for
+    Abstract base class for coordinate rotation
-    other rotation specifications
-
-    The rotation is defined by a combination of local vectors (e1/e2), (e2/e3)
-    or (e3/e1). Any nonorthogonality will be absorbed into the second vector.
-
-    For convenience, the dictionary constructor forms allow a few shortcuts:
-    - if the \c type is not otherwise specified, the type \c axes
-      is implicit
-    - if an axes specification (eg, e3/e1) is used, the coordinateRotation
-      sub-dictionary can be dropped.
-
-    Specifying the rotation by an EulerCoordinateRotation
-    (type "EulerRotation") or by a STARCDCoordinateRotation
-    (type "STARCDRotation") requires the coordinateRotation sub-dictionary.
 
     \verbatim
         coordinateRotation
         {
-            type        STARCDRotation
+            type        axesRotation
-            rotation    (0 0 90);
+            e1          (1 0 0);
+            e2          (0 1 0);
         }
     \endverbatim
 
-    - the rotation angles are in degrees, unless otherwise explictly specified:
+    Types of coordinateRotation:
+    1) axesRotation
+    2) STARCDRotation
+    3) localAxesRotation
+    4) EulerCoordinateRotation
 
-    \verbatim
-        coordinateRotation
-        {
-            type        STARCDRotation
-            degrees     false;
-            rotation    (0 0 3.141592654);
-        }
-    \endverbatim
 
-Deprecated
-    Specifying the local vectors as an \c axis (corresponding to e3) and a
-    \c direction (corresponding to e1), is allowed for backwards
-    compatibility, but this terminology is generally a bit confusing.
-    (deprecated Apr 2008)
 
 \*---------------------------------------------------------------------------*/
 
@@ -73,8 +51,10 @@ Deprecated
 
 #include "vector.H"
 #include "tensor.H"
+#include "tensorField.H"
 #include "dictionary.H"
 #include "runTimeSelectionTables.H"
+#include "objectRegistry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,57 +66,38 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class coordinateRotation
-:
-    public tensor
 {
-    // Private data
+protected:
 
-        //- the combination of local axes to be used
+    // Protected member functions
-        enum axisOrder
-        {
-            e1e2,
-            e2e3,
-            e3e1
-        };
 
-    // Private Member Functions
+        //- Transform principal
+        symmTensor transformPrincipal(const tensor&, const vector&) const;
 
-        //- Calculate transformation tensor
-        void calcTransform
-        (
-            const vector& axis1,
-            const vector& axis2,
-            const axisOrder& order = e3e1
-        );
 
 public:
 
+
     //- Runtime type information
     TypeName("coordinateRotation");
 
-    // Constructors
-
-        //- Construct null
-        coordinateRotation();
-
-        //- Construct from 2 axes
-        coordinateRotation
-        (
-            const vector& axis,
-            const vector& dir
-        );
-
-        //- Construct from dictionary
-        coordinateRotation(const dictionary&);
-
-        //- Return clone
-        autoPtr<coordinateRotation> clone() const
-        {
-            return autoPtr<coordinateRotation>(new coordinateRotation(*this));
-        }
 
     // Declare run-time constructor selection table
+    // for constructors with dictionary and objectRegistry
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        coordinateRotation,
+        objectRegistry,
+        (
+            const dictionary& dict, const objectRegistry& obr
+        ),
+        (dict, obr)
+    );
 
+
+    // Declare run-time constructor selection table
+    // for constructors with dictionary
     declareRunTimeSelectionTable
     (
         autoPtr,
@@ -151,7 +112,13 @@ public:
 
     // Selectors
 
-        //- Select constructed from Istream
+        //- Select constructed from dictionary and objectRegistry
+        static autoPtr<coordinateRotation> New
+        (
+            const dictionary& dict, const objectRegistry& obr
+        );
+
+        //- Select constructed from dictionary
         static autoPtr<coordinateRotation> New
         (
             const dictionary& dict
@@ -166,37 +133,76 @@ public:
     // Member Functions
 
         //- Reset rotation to an identity rotation
-        virtual void clear();
+        virtual void clear() = 0;
 
         //- Return local-to-global transformation tensor
-        const tensor& R() const
+        virtual const tensor& R() const = 0;
-        {
+
-            return (*this);
+        //- Return global-to-local transformation tensor
-        }
+        virtual const tensor& Rtr() const = 0;
 
         //- Return local Cartesian x-axis
-        const vector e1() const
+        virtual const vector e1() const = 0;
-        {
-            return tensor::T().x();
-        }
 
         //- Return local Cartesian y-axis
-        const vector e2() const
+        virtual const vector e2() const = 0;
-        {
-            return tensor::T().y();
-        }
 
         //- Return local Cartesian z-axis
-        const vector e3() const
+        virtual const vector e3() const = 0;
+
+        //- Return local-to-global transformation tensor
+        virtual const tensorField& Tr() const = 0;
+
+        //- Return true if the rotation tensor is uniform
+        virtual bool uniform() const
         {
-            return tensor::T().z();
+            return true;
         }
 
+        //- Transform vectorField using transformation tensor field
+        virtual tmp<vectorField> transform(const vectorField& st) const = 0;
 
-    // Member Operators
+        //- Transform vector using transformation tensor
+        virtual vector transform(const vector& st) const = 0;
 
-        //- assign from dictionary
+        //- Inverse transform vectorField using transformation tensor field
-        void operator=(const dictionary&);
+        virtual tmp<vectorField> invTransform(const vectorField& st) const = 0;
+
+        //- Inverse transform vector using transformation tensor
+        virtual vector invTransform(const vector& st) const = 0;
+
+        //- Transform tensor field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st
+        ) const = 0;
+
+        //- Transform tensor sub-field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st,
+            const labelList& cellMap
+        ) const = 0;
+
+        //- Transform tensor using transformation tensorField
+        virtual tensor transformTensor(const tensor& st) const = 0;
+
+        //- Transform vectorField using transformation tensorField and return
+        // symmetrical tensorField
+        virtual tmp<symmTensorField> transformVector
+        (
+            const vectorField& st
+        ) const = 0;
+
+        //- Transform vector using transformation tensor and return
+        // symmetrical tensor
+        virtual symmTensor transformVector(const vector& st) const = 0;
+
+
+    // Write
+
+            //- Write
+            virtual void write(Ostream&) const = 0;
 
 };
 

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotationNew.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coordinateRotation.H"
+#include "objectRegistry.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
+(
+    const dictionary& dict, const objectRegistry& obr
+)
+{
+    if (debug)
+    {
+        Pout<< "coordinateRotation::New(const dictionary&) : "
+            << "constructing coordinateRotation"
+            << endl;
+    }
+
+    word rotType = dict.lookup("type");
+
+    objectRegistryConstructorTable::iterator cstrIter =
+        objectRegistryConstructorTablePtr_->find(rotType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorIn
+        (
+            "coordinateRotation::New"
+            "("
+            "   const dictionary&, "
+            "   const objectRegistry& "
+            ")",
+            dict
+        )   << "Unknown coordinateRotation type "
+            << rotType << nl << nl
+            << "Valid coordinateRotation types are :" <<  nl
+            << "[default: axes ]"
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalIOError);
+    }
+
+    return autoPtr<coordinateRotation>(cstrIter()(dict, obr));
+}
+
+
+Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
+(
+    const dictionary& dict
+)
+{
+    if (debug)
+    {
+        Pout<< "coordinateRotation::New(const dictionary&) : "
+            << "constructing coordinateRotation"
+            << endl;
+    }
+
+    word rotType = dict.lookup("type");
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(rotType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorIn
+        (
+            "coordinateRotation::New"
+            "("
+            "   const dictionary&, "
+            ")",
+            dict
+        )   << "Unknown coordinateRotation type "
+            << rotType << nl << nl
+            << "Valid coordinateRotation types are :" <<  nl
+            << "[default: axes ]"
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalIOError);
+    }
+
+    return autoPtr<coordinateRotation>(cstrIter()(dict));
+}
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.C

@@ -0,0 +1,267 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "localAxesRotation.H"
+#include "axesRotation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "polyMesh.H"
+#include "tensorIOField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(localAxesRotation, 0);
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        localAxesRotation,
+        dictionary
+    );
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        localAxesRotation,
+        objectRegistry
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::localAxesRotation::localAxesRotation
+(
+    const dictionary& dict, const objectRegistry& orb
+)
+:
+    Rptr_(),
+    origin_(point::zero),
+    e3_(vector::zero)
+{
+    // If origin is specified in the coordinateSystem
+    if (dict.parent().found("origin"))
+    {
+        dict.parent().lookup("origin") >> origin_;
+    }
+
+    // rotation axis
+    dict.lookup("e3") >> e3_;
+
+    const polyMesh& mesh = refCast<const polyMesh>(orb);
+
+    Rptr_.reset
+    (
+        new tensorField(mesh.nCells())
+    );
+    init(dict, orb);
+}
+
+
+Foam::localAxesRotation::localAxesRotation
+(
+    const dictionary& dict
+)
+:
+    Rptr_(),
+    origin_(),
+    e3_()
+{
+    FatalErrorIn("localAxesRotation(const dictionary&)")
+        << " localAxesRotation can not be contructed from  dictionary "
+        << " use the construtctor : "
+           "("
+           "    const dictionary& dict, const objectRegistry& orb"
+           ")"
+        << exit(FatalIOError);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::localAxesRotation::clear()
+{
+    if (!Rptr_.empty())
+    {
+        Rptr_.clear();
+    }
+}
+
+
+Foam::vector Foam::localAxesRotation::transform(const vector& st) const
+{
+    notImplemented
+    (
+        "vector Foam::localAxesRotation::transform(const vector&) const"
+    );
+    return vector::zero;
+}
+
+
+Foam::vector Foam::localAxesRotation::invTransform(const vector& st) const
+{
+    notImplemented
+    (
+        "vector Foam::localAxesRotation::invTransform(const vector&) const"
+    );
+    return vector::zero;
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::localAxesRotation::transform
+(
+    const vectorField& st
+) const
+{
+    if (Rptr_->size() != st.size())
+    {
+        FatalErrorIn("localAxesRotation::transform(const vectorField&)")
+            << "vectorField st has different size to tensorField "
+            << abort(FatalError);
+    }
+
+    return (Rptr_() & st);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::localAxesRotation::invTransform
+(
+    const vectorField& st
+) const
+{
+    return (Rptr_().T() & st);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
+(
+    const tensorField& st
+) const
+{
+    if (Rptr_->size() != st.size())
+    {
+        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
+            << "tensorField st has different size to tensorField Tr"
+            << abort(FatalError);
+    }
+    return (Rptr_() & st & Rptr_().T());
+}
+
+
+Foam::tensor Foam::localAxesRotation::transformTensor
+(
+    const tensor& st
+) const
+{
+    notImplemented("tensor localAxesRotation::transformTensor() const");
+    return tensor::zero;
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
+(
+    const tensorField& st,
+    const labelList& cellMap
+) const
+{
+    if (cellMap.size() != st.size())
+    {
+        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
+            << "tensorField st has different size to tensorField Tr"
+            << abort(FatalError);
+    }
+
+    const tensorField Rtr(Rptr_().T());
+    tmp<tensorField> tt(new tensorField(cellMap.size()));
+    tensorField& t = tt();
+    forAll(cellMap, i)
+    {
+        const label cellI = cellMap[i];
+        t[i] = Rptr_()[cellI] & st[i] & Rtr[cellI];
+    }
+    return tt;
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::localAxesRotation::transformVector
+(
+    const vectorField& st
+) const
+{
+    if (Rptr_->size() != st.size())
+    {
+        FatalErrorIn("localAxesRotation::transformVector(const vectorField&)")
+            << "tensorField st has different size to tensorField Tr"
+            << abort(FatalError);
+    }
+
+    tmp<symmTensorField> tfld(new symmTensorField(Rptr_->size()));
+    symmTensorField& fld = tfld();
+
+    forAll(fld, i)
+    {
+        fld[i] = transformPrincipal(Rptr_()[i], st[i]);
+    }
+    return tfld;
+}
+
+
+Foam::symmTensor Foam::localAxesRotation::transformVector
+(
+    const vector& st
+) const
+{
+    notImplemented
+    (
+        "tensor localAxesRotation::transformVector(const vector&) const"
+    );
+    return symmTensor::zero;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::localAxesRotation::init
+(
+    const dictionary& dict,
+    const objectRegistry& obr
+)
+{
+    const polyMesh& mesh = refCast<const polyMesh>(obr);
+    forAll(mesh.cellCentres(), cellI)
+    {
+        vector dir = mesh.cellCentres()[cellI] - origin_;
+        dir /= mag(dir) + VSMALL;
+
+        Rptr_()[cellI] = axesRotation(e3_, dir).R();
+    }
+}
+
+
+void Foam::localAxesRotation::write(Ostream& os) const
+{
+     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.H

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    Foam::localAxesRotation
+
+Description
+    A local coordinate rotation.
+    Each rotational tensor is defined with two vectors (dir and e3)
+    where dir =  cellC - origin and e3 is the rotation axis.
+    Per each cell an axesRotation type of rotation is created
+
+    \verbatim
+        localAxesRotation
+        {
+            type        localAxes;
+            e3          (0 0 1);
+        }
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef localAxesRotation_H
+#define localAxesRotation_H
+
+#include "point.H"
+#include "vector.H"
+#include "coordinateRotation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class localAxesRotation Declaration
+\*---------------------------------------------------------------------------*/
+
+class localAxesRotation
+:
+    public coordinateRotation
+{
+    // Private data
+
+        //- AutoPtr to transformation tensor
+        autoPtr<tensorField> Rptr_;
+
+        //- Origin of the coordinate system
+        point origin_;
+
+        //- Rotation axis
+        vector e3_;
+
+
+    // Private members
+
+        //- Init transformation tensor
+        void init(const dictionary& dict, const objectRegistry& obr);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("localAxesRotation");
+
+    // Constructors
+
+        //- Construct from dictionary and objectRegistry
+        localAxesRotation(const dictionary&, const objectRegistry&);
+
+        //- Construct from dictionary
+        localAxesRotation(const dictionary&);
+
+        //- Return clone
+        autoPtr<localAxesRotation> clone() const
+        {
+            return autoPtr<localAxesRotation>(new localAxesRotation(*this));
+        }
+
+
+    //- Destructor
+    virtual ~localAxesRotation()
+    {}
+
+
+    // Member Functions
+
+        //- Reset rotation to an identity rotation
+        virtual void clear();
+
+        //- Return local-to-global transformation tensor
+        virtual const tensor& R() const
+        {
+            notImplemented("const tensor& localAxesRotation::R() const");
+            return tensor::zero;
+        }
+
+        //- Return global-to-local transformation tensor
+        virtual const tensor& Rtr() const
+        {
+             notImplemented("const tensor& localAxesRotation::Rtr() const");
+             return tensor::zero;
+        }
+
+        //- Return local Cartesian x-axis
+        virtual const vector e1() const
+        {
+            notImplemented("const tensor& localAxesRotation::e1() const");
+            return vector::zero;
+        }
+
+        //- Return local Cartesian y-axis
+        virtual const vector e2() const
+        {
+            notImplemented("const tensor& localAxesRotation::e2() const");
+            return vector::zero;
+        }
+
+        //- Return local Cartesian z-axis
+        virtual const vector e3() const
+        {
+            return e3_;
+        }
+
+        virtual const tensorField& Tr() const
+        {
+            return Rptr_();
+        }
+
+        //- Transform vectorField using transformation tensor field
+        virtual tmp<vectorField> transform(const vectorField& st) const;
+
+        //- Transform vector using transformation tensor
+        virtual vector transform(const vector& st) const;
+
+        //- Inverse transform vectorField using transformation tensor field
+        virtual tmp<vectorField> invTransform(const vectorField& st) const;
+
+        //- Inverse transform vector using transformation tensor
+        virtual vector invTransform(const vector& st) const;
+
+        //- Return if the rotation is uniform
+        virtual bool uniform() const
+        {
+            return false;
+        }
+
+        //- Transform tensor field using transformation tensorField
+        virtual tmp<tensorField> transformTensor(const tensorField& st) const;
+
+        //- Transform tensor using transformation tensorField
+        virtual tensor transformTensor(const tensor& st) const;
+
+        //- Transform tensor sub-field using transformation tensorField
+        virtual tmp<tensorField> transformTensor
+        (
+            const tensorField& st,
+            const labelList& cellMap
+        ) const;
+
+        //- Transform vectorField using transformation tensorField and return
+        // symmetrical tensorField
+        virtual tmp<symmTensorField> transformVector
+        (
+            const vectorField& st
+        ) const;
+
+        //- Transform vector using transformation tensor and return
+        // symmetrical tensor (R & st & R.T())
+        virtual symmTensor transformVector(const vector& st) const;
+
+
+    // Write
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/coordinateSystems/coordinateSystem.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "IOstream.H"
+#include "axesRotation.H"
 #include "coordinateSystem.H"
 #include "coordinateSystems.H"
 #include "addToRunTimeSelectionTable.H"
@@ -34,18 +35,16 @@ namespace Foam
 {
     defineTypeNameAndDebug(coordinateSystem, 0);
     defineRunTimeSelectionTable(coordinateSystem, dictionary);
-    defineRunTimeSelectionTable(coordinateSystem, origRotation);
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::coordinateSystem::coordinateSystem()
 :
-    name_(type()),
+    name_(),
     note_(),
     origin_(point::zero),
-    R_(),
+    R_(new axesRotation(sphericalTensor::I))
-    Rtr_(sphericalTensor::I)
 {}
 
 
@@ -58,8 +57,7 @@ Foam::coordinateSystem::coordinateSystem
     name_(name),
     note_(),
     origin_(cs.origin_),
-    R_(cs.R_),
+    R_(const_cast<coordinateRotation*>(&cs.R()))
-    Rtr_(R_.T())
 {}
 
 
@@ -73,8 +71,7 @@ Foam::coordinateSystem::coordinateSystem
     name_(name),
     note_(),
     origin_(origin),
-    R_(cr),
+    R_(const_cast<coordinateRotation*>(&cr))
-    Rtr_(R_.T())
 {}
 
 
@@ -89,8 +86,7 @@ Foam::coordinateSystem::coordinateSystem
     name_(name),
     note_(),
     origin_(origin),
-    R_(axis, dirn),
+    R_(new axesRotation(axis, dirn))
-    Rtr_(R_.T())
 {}
 
 
@@ -103,37 +99,35 @@ Foam::coordinateSystem::coordinateSystem
     name_(name),
     note_(),
     origin_(point::zero),
-    R_(),
+    R_()
-    Rtr_(sphericalTensor::I)
 {
-    operator=(dict);
+    init(dict);
 }
 
 
 Foam::coordinateSystem::coordinateSystem(const dictionary& dict)
 :
-    name_(type()),
+    name_(),
     note_(),
     origin_(point::zero),
-    R_(),
+    R_()
-    Rtr_(sphericalTensor::I)
 {
-    operator=(dict);
+    init(dict);
 }
 
 
 Foam::coordinateSystem::coordinateSystem
 (
-    const dictionary& dict,
+    const objectRegistry& obr,
-    const objectRegistry& obr
+    const dictionary& dict
 )
 :
-    name_(type()),
+    name_(),
     note_(),
     origin_(point::zero),
-    R_(),
+    R_()
-    Rtr_(sphericalTensor::I)
 {
+
     const entry* entryPtr = dict.lookupEntryPtr(typeName_(), false, false);
 
     // non-dictionary entry is a lookup into global coordinateSystems
@@ -170,7 +164,7 @@ Foam::coordinateSystem::coordinateSystem
     }
     else
     {
-        operator=(dict);
+        init(dict, obr);
     }
 }
 
@@ -180,11 +174,10 @@ Foam::coordinateSystem::coordinateSystem(Istream& is)
     name_(is),
     note_(),
     origin_(point::zero),
-    R_(),
+    R_()
-    Rtr_(sphericalTensor::I)
 {
     dictionary dict(is);
-    operator=(dict);
+    init(dict);
 }
 
 
@@ -215,8 +208,8 @@ Foam::dictionary Foam::coordinateSystem::dict(bool ignoreType) const
     }
 
     dict.add("origin", origin_);
-    dict.add("e1", e1());
+    dict.add("e1", R_->e1());
-    dict.add("e3", e3());
+    dict.add("e3", R_->e3());
 
     return dict;
 }
@@ -230,11 +223,11 @@ Foam::vector Foam::coordinateSystem::localToGlobal
 {
     if (translate)
     {
-        return (R_ & local) + origin_;
+        return (R_->transform(local)) + origin_;
     }
     else
     {
-        return (R_ & local);
+        return R_->transform(local);
     }
 }
 
@@ -247,11 +240,11 @@ Foam::tmp<Foam::vectorField> Foam::coordinateSystem::localToGlobal
 {
     if (translate)
     {
-        return (R_ & local) + origin_;
+        return (R_->transform(local)) + origin_;
     }
     else
     {
-        return (R_ & local);
+        return R_->transform(local);
     }
 }
 
@@ -264,11 +257,11 @@ Foam::vector Foam::coordinateSystem::globalToLocal
 {
     if (translate)
     {
-        return (Rtr_ & (global - origin_));
+        return R_->invTransform(global - origin_);
     }
     else
     {
-        return (Rtr_ & global);
+        return R_->invTransform(global);
     }
 }
 
@@ -281,11 +274,11 @@ Foam::tmp<Foam::vectorField> Foam::coordinateSystem::globalToLocal
 {
     if (translate)
     {
-        return (Rtr_ & (global - origin_));
+        return R_->invTransform(global - origin_);
     }
     else
     {
-        return (Rtr_ & global);
+        return R_->invTransform(global);
     }
 }
 
@@ -294,15 +287,14 @@ void Foam::coordinateSystem::clear()
 {
     note_.clear();
     origin_ = point::zero;
-    R_.clear();
+    R_->clear();
-    Rtr_ = sphericalTensor::I;
 }
 
 
 void Foam::coordinateSystem::write(Ostream& os) const
 {
-    os  << type()
+    os  << type() << " origin: " << origin() << nl;
-        << " origin: " << origin() << " e1: " << e1() << " e3: " << e3();
+    R_->write(os);
 }
 
 
@@ -314,11 +306,8 @@ void Foam::coordinateSystem::writeDict(Ostream& os, bool subDict) const
             << indent << token::BEGIN_BLOCK << incrIndent << nl;
     }
 
-    // only write type for derived types
-    if (type() != typeName_())
-    {
     os.writeKeyword("type") << type() << token::END_STATEMENT << nl;
-    }
+
 
     // The note entry is optional
     if (note_.size())
@@ -327,8 +316,7 @@ void Foam::coordinateSystem::writeDict(Ostream& os, bool subDict) const
     }
 
     os.writeKeyword("origin") << origin_ << token::END_STATEMENT << nl;
-    os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
+    R_->write(os);
-    os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
 
     if (subDict)
     {
@@ -339,7 +327,20 @@ void Foam::coordinateSystem::writeDict(Ostream& os, bool subDict) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void Foam::coordinateSystem::operator=(const dictionary& rhs)
+void Foam::coordinateSystem::init(const dictionary& rhs)
+{
+    rhs.lookup("origin") >> origin_;
+    note_.clear();
+    rhs.readIfPresent("note", note_);
+    R_.reset(coordinateRotation::New(rhs.subDict("coordinateRotation")).ptr());
+}
+
+
+void Foam::coordinateSystem::init
+(
+    const dictionary& rhs,
+    const objectRegistry& obr
+)
 {
     if (debug)
     {
@@ -347,34 +348,16 @@ void Foam::coordinateSystem::operator=(const dictionary& rhs)
             << "assign from " << rhs << endl;
     }
 
-    // allow as embedded sub-dictionary "coordinateSystem"
+    rhs.lookup("origin") >> origin_;
-    const dictionary& dict =
-    (
-        rhs.found(typeName_())
-      ? rhs.subDict(typeName_())
-      : rhs
-    );
-
-    // unspecified origin is (0 0 0)
-    origin_ = point::zero;
-    dict.readIfPresent("origin", origin_);
 
     // The note entry is optional
     note_.clear();
     rhs.readIfPresent("note", note_);
 
-    // specify via coordinateRotation sub-dictionary
+    R_.reset
-    if (dict.found("coordinateRotation"))
+    (
-    {
+        coordinateRotation::New(rhs.subDict("coordinateRotation"), obr).ptr()
-        R_  = coordinateRotation::New(dict.subDict("coordinateRotation"))();
+    );
-    }
-    else
-    {
-        // let coordinateRotation constructor extract the axes specification
-        R_ = coordinateRotation(dict);
-    }
-
-    Rtr_ = R_.T();
 }
 
 
@@ -382,7 +365,12 @@ void Foam::coordinateSystem::operator=(const dictionary& rhs)
 
 bool Foam::operator!=(const coordinateSystem& a, const coordinateSystem& b)
 {
-    return (a.origin() != b.origin() || a.R() != b.R() || a.type() != b.type());
+    return
+    (
+        a.origin() != b.origin()
+     || a.R().R() != b.R().R()
+     || a.type() != b.type()
+    );
 }
 
 

src/meshTools/coordinateSystems/coordinateSystem.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,60 +25,15 @@ Class
     Foam::coordinateSystem
 
 Description
-    A cartesian coordinate system and the base class for other coordinate
+    Base class for other coordinate
     system specifications.
 
     All systems are defined by an origin point and a coordinateRotation.
-    For convenience, the dictionary constructor forms allow a few shortcuts:
-    - the default origin corresponds to <em>(0 0 0)</em>
-    - if the \c type is not otherwise specified, a Cartesian coordinateSystem
-      is implicit
 
     \verbatim
-        flipped
-        {
-            origin  (0 0 0);
-            coordinateRotation
-            {
-                type        STARCDRotation;
-                rotation    (0 0 90);
-            }
-        }
-    \endverbatim
-
-    - if an axes specification (eg, e3/e1) is used, the coordinateRotation
-      sub-dictionary can be dropped.
-
-    \verbatim
-        flipped     // the same, specified as axes
-        {
-            origin  (0 0 0);
-            coordinateRotation
-            {
-                type    axes;
-                e3      (1 0 0);
-                e1      (0 0 -1);
-            }
-        }
-        flipped     // the same, using all the shortcuts
-        {
-            e3      (1 0 0);
-            e1      (0 0 -1);
-        }
-    \endverbatim
-
-    - if a sub-dictionary coordinateSystem is found within the dictionary, it
-      will be used. This provides a convenient means of embedding
-      coordinateSystem information in another dictionary.
-      This is used, for example, in the porousZones:
-
-    \verbatim
-        1
-        (
-        cat1
-        {
         coordinateSystem
         {
+            type    cartesian;
             origin  (0 0 0);
             coordinateRotation
             {
@@ -86,39 +41,20 @@ Description
                 rotation    (0 0 90);
             }
         }
-            porosity        0.781;
-            Darcy
-            {
-                d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
-                f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
-            }
-        }
-        )
     \endverbatim
 
-    - additionally, if the coordinateSystem points to a plain entry,
+    Types of coordinateRotation:
-      it can be used to reference one of the global coordinateSystems
+    1) axesRotation
+    2) STARCDRotation
+    3) localAxesRotation
+    4) EulerCoordinateRotation
 
-    \verbatim
+    Type of coordinates:
-        1
+    1) cartesian
-        (
+    2) cylindricalCS
-        cat1
-        {
-            coordinateSystem  system_10;
-            porosity        0.781;
-            Darcy
-            {
-                d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
-                f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
-            }
-        }
-        )
-    \endverbatim
-    For this to work correctly, the coordinateSystem constructor must be
-    supplied with both a dictionary and an objectRegistry.
 
 See Also
-    coordinateSystems and coordinateSystems::New
+    coordinateSystem and coordinateSystem::New
 
 SourceFiles
     coordinateSystem.C
@@ -136,6 +72,7 @@ SourceFiles
 #include "tmp.H"
 #include "coordinateRotation.H"
 #include "objectRegistry.H"
+#include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -151,19 +88,16 @@ class coordinateSystem
     // Private data
 
         //- Name of coordinate system
-        mutable word name_;
+        word name_;
 
         //- Optional note
-        mutable string note_;
+        string note_;
 
         //- Origin
-        mutable point origin_;
+        point origin_;
 
         //- Local-to-Global transformation tensor
-        coordinateRotation R_;
+        autoPtr<coordinateRotation> R_;
-
-        //- Global-to-Local transformation tensor
-        tensor Rtr_;
 
 
 protected:
@@ -194,6 +128,12 @@ protected:
             bool translate
         ) const;
 
+        //- Init from dict and obr
+        void init(const dictionary&);
+
+        //- Init from dictionary
+        void init(const dictionary&, const objectRegistry&);
+
 
 public:
 
@@ -238,12 +178,13 @@ public:
 
         //- Construct from dictionary (default name)
         //  With the ability to reference global coordinateSystems
-        coordinateSystem(const dictionary&, const objectRegistry&);
+        coordinateSystem(const objectRegistry&, const dictionary&);
 
         //- Construct from Istream
         //  The Istream contains a word followed by a dictionary
         coordinateSystem(Istream&);
 
+
     //- Return clone
     autoPtr<coordinateSystem> clone() const
     {
@@ -252,49 +193,32 @@ public:
 
 
     // Declare run-time constructor selection table
-
     declareRunTimeSelectionTable
     (
         autoPtr,
         coordinateSystem,
         dictionary,
         (
-                const word& name,
+            const objectRegistry& obr,
             const dictionary& dict
         ),
-            (name, dict)
+        (obr, dict)
-        );
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            coordinateSystem,
-            origRotation,
-            (
-                const word& name,
-                const point& origin,
-                const coordinateRotation& cr
-            ),
-            (name, origin, cr)
     );
 
 
     // Selectors
 
+        //- Select constructed from dictionary and objectRegistry
+        static autoPtr<coordinateSystem> New
+        (
+            const objectRegistry& obr,
+            const dictionary& dict
+        );
+
         //- Select constructed from dictionary
         static autoPtr<coordinateSystem> New
         (
-            const word& name,
+            const dictionary& dict
-            const dictionary&
-        );
-
-        //- Select constructed from origin and rotation
-        static autoPtr<coordinateSystem> New
-        (
-            const word& coordType,
-            const word& name,
-            const point& origin,
-            const coordinateRotation&
         );
 
         //- Select constructed from Istream
@@ -307,6 +231,7 @@ public:
 
     // Member Functions
 
+
         // Access
 
             //- Return name
@@ -333,49 +258,18 @@ public:
                 return origin_;
             }
 
-            //- Return coordinate rotation
+            //- Return const reference to coordinate rotation
-            const coordinateRotation& rotation() const
+            const coordinateRotation& R() const
             {
-                return R_;
+                return R_();
             }
 
-            //- Return local-to-global transformation tensor
+            //- Return non const reference to coordinate rotation
-            const tensor& R() const
+            coordinateRotation& R()
             {
-                return R_;
+                return R_();
             }
 
-            //- Return local Cartesian x-axis
-            const vector e1() const
-            {
-                return Rtr_.x();
-            }
-
-            //- Return local Cartesian y-axis
-            const vector e2() const
-            {
-                return Rtr_.y();
-            }
-
-            //- Return local Cartesian z-axis
-            const vector e3() const
-            {
-                return Rtr_.z();
-            }
-
-            //- Return axis (e3: local Cartesian z-axis)
-            // \deprecated method e3 is preferred (deprecated Apr 2008)
-            const vector axis() const
-            {
-                return Rtr_.z();
-            }
-
-            //- Return direction (e1: local Cartesian x-axis)
-            // \deprecated method e1 is preferred (deprecated Apr 2008)
-            const vector direction() const
-            {
-                return Rtr_.x();
-            }
 
             //- Return as dictionary of entries
             //  \param [in] ignoreType drop type (cartesian, cylindrical, etc)
@@ -386,7 +280,7 @@ public:
         // Edit
 
             //- Rename
-            virtual void rename(const word& newName)
+            void rename(const word& newName)
             {
                 name_ = newName;
             }
@@ -408,7 +302,7 @@ public:
             virtual void write(Ostream&) const;
 
             //- Write dictionary
-            virtual void writeDict(Ostream&, bool subDict=true) const;
+            void writeDict(Ostream&, bool subDict=true) const;
 
 
         // Transformations
@@ -472,10 +366,6 @@ public:
 
     // Member Operators
 
-        //- assign from dictionary
-        void operator=(const dictionary&);
-
-
         // friend Operators
 
             friend bool operator!=

src/meshTools/coordinateSystems/coordinateSystemNew.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,7 @@ License
 
 Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
 (
-    const word& name,
+    const objectRegistry& obr,
     const dictionary& dict
 )
 {
@@ -41,17 +41,8 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
             << endl;
     }
 
-    // construct base class directly, also allow 'cartesian' as an alias
+    const dictionary& coordDict = dict.subDict(typeName_());
-    word coordType(typeName_());
+    word coordType = coordDict.lookup("type");
-    if
-    (
-        !dict.readIfPresent("type", coordType)
-     || coordType == typeName_()
-     || coordType == "cartesian"
-    )
-    {
-        return autoPtr<coordinateSystem>(new coordinateSystem(name, dict));
-    }
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(coordType);
@@ -60,20 +51,54 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
     {
         FatalIOErrorIn
         (
-            "coordinateSystem::New(const word&, const dictionary&)",
+            "coordinateSystem::New(const objectRegistry&, const dictionary&)",
             dict
         )   << "Unknown coordinateSystem type "
             << coordType << nl << nl
             << "Valid coordinateSystem types are :" << nl
-            << "[default: " << typeName_() << "]"
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalIOError);
     }
 
-    return autoPtr<coordinateSystem>(cstrIter()(name, dict));
+    return autoPtr<coordinateSystem>(cstrIter()(obr, coordDict));
 }
 
 
+Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
+(
+    const dictionary& dict
+)
+{
+    if (debug)
+    {
+        Pout<< "coordinateSystem::New(cconst dictionary&) : "
+            << "constructing coordinateSystem"
+            << endl;
+    }
+
+    const dictionary& coordDict = dict.subDict(typeName_());
+    word coordType = coordDict.lookup("type");
+/*
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(coordType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorIn
+        (
+            "coordinateSystem::New(const dictionary&)",
+            dict
+        )   << "Unknown coordinateSystem type "
+            << coordType << nl << nl
+            << "Valid coordinateSystem types are :" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalIOError);
+    }
+*/
+    return autoPtr<coordinateSystem>(new coordinateSystem(coordDict));
+}
+
+/*
 Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
 (
     const word& coordType,
@@ -109,7 +134,7 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
 
     return autoPtr<coordinateSystem>(cstrIter()(name, origin, cr));
 }
-
+*/
 
 Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
 (
@@ -119,7 +144,26 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
     const word name(is);
     const dictionary dict(is);
 
-    return New(name, dict);
+    word coordType = dict.lookup("type");
+/*
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(coordType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorIn
+        (
+            "coordinateSystem::New(Istream& is)",
+            dict
+        )   << "Unknown coordinateSystem type "
+            << coordType << nl << nl
+            << "Valid coordinateSystem types are :" << nl
+            << "[default: " << typeName_() << "]"
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalIOError);
+    }
+*/
+    return autoPtr<coordinateSystem>(new coordinateSystem(name, dict));
 }
 
 

src/meshTools/coordinateSystems/coordinateSystems.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,25 @@ Note
     Mixing normal constructors and the coordinateSystems::New constructor
     may yield unexpected results.
 
+    \verbatim
+        1
+        (
+        cat1
+        {
+            coordinateSystem  system_10;
+            porosity        0.781;
+            Darcy
+            {
+                d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
+                f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
+            }
+        }
+        )
+    \endverbatim
+
+    For this to work correctly, the coordinateSystem constructor must be
+    supplied with both a dictionary and an objectRegistry.
+
 SourceFiles
     coordinateSystems.C
 

src/meshTools/coordinateSystems/cylindricalCS.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,6 @@ namespace Foam
 {
     defineTypeNameAndDebug(cylindricalCS, 0);
     addToRunTimeSelectionTable(coordinateSystem, cylindricalCS, dictionary);
-    addToRunTimeSelectionTable(coordinateSystem, cylindricalCS, origRotation);
 }
 
 
@@ -109,6 +108,23 @@ Foam::cylindricalCS::cylindricalCS
 {}
 
 
+Foam::cylindricalCS::cylindricalCS
+(
+    const objectRegistry& obr,
+    const dictionary& dict
+)
+:
+    coordinateSystem(obr, dict),
+    inDegrees_(dict.lookupOrDefault("degrees", true))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cylindricalCS::~cylindricalCS()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::cylindricalCS::inDegrees() const

src/meshTools/coordinateSystems/cylindricalCS.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,7 @@ protected:
 
     // Protected Member Functions
 
+
         //- Convert from local coordinate system to the global Cartesian system
         //  with optional translation for the origin
         virtual vector localToGlobal(const vector&, bool translate) const;
@@ -131,8 +132,15 @@ public:
             const bool inDegrees=true
         );
 
-        //- Construct from dictionary
+        //- Construct from dictionary and name
-        cylindricalCS(const word& name, const dictionary&);
+        cylindricalCS(const word&, const dictionary&);
+
+        //- Construct from dictionary and objectRegistry
+        cylindricalCS(const objectRegistry&, const dictionary&);
+
+
+    //- Destructor
+    virtual ~cylindricalCS();
 
 
     // Member Functions

src/meshTools/coordinateSystems/parabolicCylindricalCS.C

@@ -1,177 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "parabolicCylindricalCS.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(parabolicCylindricalCS, 0);
-    addToRunTimeSelectionTable
-    (
-        coordinateSystem,
-        parabolicCylindricalCS,
-        dictionary
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::parabolicCylindricalCS::parabolicCylindricalCS()
-:
-    coordinateSystem()
-{}
-
-
-Foam::parabolicCylindricalCS::parabolicCylindricalCS
-(
-    const word& name,
-    const point& origin,
-    const coordinateRotation& cr
-)
-:
-    coordinateSystem(name, origin, cr)
-{}
-
-
-Foam::parabolicCylindricalCS::parabolicCylindricalCS
-(
-    const word& name,
-    const dictionary& dict
-)
-:
-    coordinateSystem(name, dict)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::vector Foam::parabolicCylindricalCS::localToGlobal
-(
-    const vector& local,
-    bool translate
-) const
-{
-    // Notation: u = local.x() v = local.y() z = local.z();
-    if (local.y() < 0.0)
-    {
-        FatalErrorIn
-        (
-            "parabolicCylindricalCS::localToGlobal(const vector&, bool) const"
-        )
-            << "parabolic cylindrical coordinates v < 0"
-            << abort(FatalError);
-    }
-
-    return coordinateSystem::localToGlobal
-    (
-        vector
-        (
-            0.5*(sqr(local.x()) - sqr(local.y())),
-            local.x()*local.y(),
-            local.z()
-        ),
-        translate
-    );
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::parabolicCylindricalCS::localToGlobal
-(
-    const vectorField& local,
-    bool translate
-) const
-{
-    if (min(local.component(vector::Y)) < 0.0)
-    {
-        FatalErrorIn
-        (
-            "parabolicCylindricalCS::localToGlobal"
-            "(const vectorField&, bool) const"
-        )   << "parabolic cylindrical coordinates v < 0"
-            << abort(FatalError);
-    }
-
-    vectorField lc(local.size());
-    lc.replace
-    (
-        vector::X,
-        0.5*
-        (
-            sqr(local.component(vector::X))
-          - sqr(local.component(vector::Y))
-        )
-    );
-
-    lc.replace
-    (
-        vector::Y,
-        local.component(vector::X) * local.component(vector::Y)
-    );
-
-    lc.replace
-    (
-        vector::Z,
-        local.component(vector::Z)
-    );
-
-    return coordinateSystem::localToGlobal(lc, translate);
-}
-
-
-Foam::vector Foam::parabolicCylindricalCS::globalToLocal
-(
-    const vector& global,
-    bool translate
-) const
-{
-    notImplemented
-    (
-        "parabolicCylindricalCS::globalToLocal(const vector&, bool) const"
-    );
-
-    return vector::zero;
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::parabolicCylindricalCS::globalToLocal
-(
-    const vectorField& global,
-    bool translate
-) const
-{
-    notImplemented
-    (
-        "parabolicCylindricalCS::globalToLocal(const vectorField&, bool) const"
-    );
-
-    return tmp<vectorField>(vectorField::null());
-}
-
-
-// ************************************************************************* //

src/meshTools/coordinateSystems/parabolicCylindricalCS.H

@@ -1,120 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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/>.
-
-Class
-    Foam::parabolicCylindricalCS
-
-Description
-    Parabolic cylindrical coordinate system.
-
-    Notation: u = a.x() v = a.y() z = a.z();
-
-Note
-    The maintenance of this class may lag that of the main types.
-
-SourceFiles
-    parabolicCylindricalCS.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef parabolicCylindricalCS_H
-#define parabolicCylindricalCS_H
-
-#include "coordinateSystem.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class parabolicCylindricalCS Declaration
-\*---------------------------------------------------------------------------*/
-
-class parabolicCylindricalCS
-:
-    public coordinateSystem
-{
-
-protected:
-
-    // Protected Member Functions
-
-        //- Convert from local coordinate system to the global Cartesian system
-        //  with optional translation for the origin
-        virtual vector localToGlobal(const vector&, bool translate) const;
-
-        //- Convert from local coordinate system to the global Cartesian system
-        //  with optional translation for the origin
-        virtual tmp<vectorField> localToGlobal
-        (
-            const vectorField&,
-            bool translate
-        ) const;
-
-        //- Convert from global Cartesian system to the local coordinate system
-        //  with optional translation for the origin
-        virtual vector globalToLocal(const vector&, bool translate) const;
-
-        //- Convert from global Cartesian system to the local coordinate system
-        //  with optional translation for the origin
-        virtual tmp<vectorField> globalToLocal
-        (
-            const vectorField&,
-            bool translate
-        ) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("parabolicCylindrical");
-
-
-    // Constructors
-
-        //- Construct null
-        parabolicCylindricalCS();
-
-        //- Construct from origin and rotation
-        parabolicCylindricalCS
-        (
-            const word& name,
-            const point& origin,
-            const coordinateRotation&
-        );
-
-        //- Construct from dictionary
-        parabolicCylindricalCS(const word&, const dictionary&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/meshTools/coordinateSystems/sphericalCS.C

@@ -1,244 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sphericalCS.H"
-
-#include "one.H"
-#include "mathematicalConstants.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(sphericalCS, 0);
-    addToRunTimeSelectionTable(coordinateSystem, sphericalCS, dictionary);
-    addToRunTimeSelectionTable(coordinateSystem, sphericalCS, origRotation);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::sphericalCS::sphericalCS(const bool inDegrees)
-:
-    coordinateSystem(),
-    inDegrees_(inDegrees)
-{}
-
-
-Foam::sphericalCS::sphericalCS
-(
-    const coordinateSystem& cs,
-    const bool inDegrees
-)
-:
-    coordinateSystem(cs),
-    inDegrees_(inDegrees)
-{}
-
-
-Foam::sphericalCS::sphericalCS
-(
-    const word& name,
-    const coordinateSystem& cs,
-    const bool inDegrees
-)
-:
-    coordinateSystem(name, cs),
-    inDegrees_(inDegrees)
-{}
-
-
-Foam::sphericalCS::sphericalCS
-(
-    const word& name,
-    const point& origin,
-    const coordinateRotation& cr,
-    const bool inDegrees
-)
-:
-    coordinateSystem(name, origin, cr),
-    inDegrees_(inDegrees)
-{}
-
-
-Foam::sphericalCS::sphericalCS
-(
-    const word& name,
-    const point& origin,
-    const vector& axis,
-    const vector& dirn,
-    const bool inDegrees
-)
-:
-    coordinateSystem(name, origin, axis, dirn),
-    inDegrees_(inDegrees)
-{}
-
-
-Foam::sphericalCS::sphericalCS
-(
-    const word& name,
-    const dictionary& dict
-)
-:
-    coordinateSystem(name, dict),
-    inDegrees_(dict.lookupOrDefault("degrees", true))
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sphericalCS::inDegrees() const
-{
-    return inDegrees_;
-}
-
-
-bool& Foam::sphericalCS::inDegrees()
-{
-    return inDegrees_;
-}
-
-
-Foam::vector Foam::sphericalCS::localToGlobal
-(
-    const vector& local,
-    bool translate
-) const
-{
-    scalar r = local.x();
-    const scalar theta
-    (
-        local.y()
-       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
-    );
-    const scalar phi
-    (
-        local.z()
-       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
-    );
-
-    return coordinateSystem::localToGlobal
-    (
-        vector(r*cos(theta)*sin(phi), r*sin(theta)*sin(phi), r*cos(phi)),
-        translate
-    );
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::sphericalCS::localToGlobal
-(
-    const vectorField& local,
-    bool translate
-) const
-{
-    const scalarField r(local.component(vector::X));
-    const scalarField theta
-    (
-        local.component(vector::Y)
-       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
-    );
-    const scalarField phi
-    (
-        local.component(vector::Z)
-       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
-    );
-
-    vectorField lc(local.size());
-    lc.replace(vector::X, r*cos(theta)*sin(phi));
-    lc.replace(vector::Y, r*sin(theta)*sin(phi));
-    lc.replace(vector::Z, r*cos(phi));
-
-    return coordinateSystem::localToGlobal(lc, translate);
-}
-
-
-Foam::vector Foam::sphericalCS::globalToLocal
-(
-    const vector& global,
-    bool translate
-) const
-{
-    const vector lc = coordinateSystem::globalToLocal(global, translate);
-    const scalar r = mag(lc);
-
-    return vector
-    (
-        r,
-        atan2
-        (
-            lc.y(), lc.x()
-        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0),
-        acos
-        (
-            lc.z()/(r + SMALL)
-        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0)
-    );
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::sphericalCS::globalToLocal
-(
-    const vectorField& global,
-    bool translate
-) const
-{
-    const vectorField lc(coordinateSystem::globalToLocal(global, translate));
-    const scalarField r(mag(lc));
-
-    tmp<vectorField> tresult(new vectorField(lc.size()));
-    vectorField& result = tresult();
-
-    result.replace
-    (
-        vector::X, r
-
-    );
-
-    result.replace
-    (
-        vector::Y,
-        atan2
-        (
-            lc.component(vector::Y),
-            lc.component(vector::X)
-        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0)
-    );
-
-    result.replace
-    (
-        vector::Z,
-        acos
-        (
-            lc.component(vector::Z)/(r + SMALL)
-        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0)
-    );
-
-    return tresult;
-}
-
-
-// ************************************************************************* //

src/meshTools/coordinateSystems/toroidalCS.C

@@ -1,184 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "toroidalCS.H"
-#include "addToRunTimeSelectionTable.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(toroidalCS, 0);
-    addToRunTimeSelectionTable(coordinateSystem, toroidalCS, dictionary);
-}
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-Foam::toroidalCS::toroidalCS
-(
-    const word& name,
-    const point& origin,
-    const coordinateRotation& cr,
-    const scalar radius
-)
-:
-    coordinateSystem(name, origin, cr),
-    radius_(radius)
-{}
-
-
-Foam::toroidalCS::toroidalCS
-(
-    const word& name,
-    const dictionary& dict
-)
-:
-    coordinateSystem(name, dict),
-    radius_(readScalar(dict.lookup("radius")))
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::vector Foam::toroidalCS::localToGlobal
-(
-    const vector& local,
-    bool translate
-) const
-{
-    // Notation: r = local.x()
-    scalar theta = degToRad(local.y());
-    scalar phi = degToRad(local.z());
-
-    scalar rprime = radius_ + local.x()*sin(phi);
-
-    if ((local.x()*sin(phi)) > (radius_))
-    {
-        FatalErrorIn("toroidalCS::toGlobal(vector) const")
-            << "Badly defined toroidal coordinates"
-            << abort(FatalError);
-    }
-
-    return coordinateSystem::localToGlobal
-    (
-        vector(rprime*cos(theta), rprime*sin(theta), local.x()*cos(phi)),
-        translate
-    );
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::toroidalCS::localToGlobal
-(
-    const vectorField& local,
-    bool translate
-) const
-{
-    const scalarField r
-    (
-        local.component(vector::X)
-    );
-
-    const scalarField theta
-    (
-        local.component(vector::Y)*constant::mathematical::pi/180.0
-    );
-
-    const scalarField phi
-    (
-        local.component(vector::Z)*constant::mathematical::pi/180.0
-    );
-
-    const scalarField rprime
-    (
-        radius_ + r*sin(phi)
-    );
-
-    vectorField lc(local.size());
-    lc.replace(vector::X, rprime*cos(theta));
-    lc.replace(vector::Y, rprime*sin(theta));
-    lc.replace(vector::Z, r*cos(phi));
-
-    return coordinateSystem::localToGlobal(lc, translate);
-}
-
-
-Foam::vector Foam::toroidalCS::globalToLocal
-(
-    const vector& global,
-    bool translate
-) const
-{
-    notImplemented
-    (
-        "toroidalCS::globalToLocal(const vector&, bool) const"
-    );
-
-    return vector::zero;
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::toroidalCS::globalToLocal
-(
-    const vectorField& global,
-    bool translate
-) const
-{
-    notImplemented
-    (
-        "toroidalCS::globalToLocal(const vectorField&, bool) const"
-    );
-
-    return tmp<vectorField>(vectorField::null());
-}
-
-
-void Foam::toroidalCS::write(Ostream& os) const
-{
-    coordinateSystem::write(os);
-    os << "radius: " << radius() << endl;
-}
-
-
-void Foam::toroidalCS::writeDict(Ostream& os, bool subDict) const
-{
-    if (subDict)
-    {
-        os  << indent << name() << nl
-            << indent << token::BEGIN_BLOCK << incrIndent << nl;
-    }
-
-    coordinateSystem::writeDict(os, false);
-    os.writeKeyword("radius") << radius() << token::END_STATEMENT << nl;
-
-    if (subDict)
-    {
-        os << decrIndent << indent << token::END_BLOCK << endl;
-    }
-}
-
-
-// ************************************************************************* //

src/meshTools/coordinateSystems/toroidalCS.H

@@ -1,133 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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/>.
-
-Class
-    Foam::toroidalCS
-
-Description
-    Toroidal coordinate system, always in degrees
-
-Note
-    The maintenance of this class may lag that of the main types.
-
-SourceFiles
-    toroidalCS.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef toroidalCS_H
-#define toroidalCS_H
-
-#include "coordinateSystem.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class toroidalCS Declaration
-\*---------------------------------------------------------------------------*/
-
-class toroidalCS
-:
-    public coordinateSystem
-{
-    // Private data
-
-        //- Radius of the torus
-        scalar radius_;
-
-    // Private Member Functions
-
-        //- Convert from local coordinate system to the global Cartesian system
-        //  with optional translation for the origin
-        virtual vector localToGlobal(const vector&, bool translate) const;
-
-        //- Convert from local coordinate system to the global Cartesian system
-        //  with optional translation for the origin
-        virtual tmp<vectorField> localToGlobal
-        (
-            const vectorField&,
-            bool translate
-        ) const;
-
-        //- Convert from global Cartesian system to the local coordinate system
-        //  with optional translation for the origin
-        virtual vector globalToLocal(const vector&, bool translate) const;
-
-        //- Convert from global Cartesian system to the local coordinate system
-        //  with optional translation for the origin
-        virtual tmp<vectorField> globalToLocal
-        (
-            const vectorField&,
-            bool translate
-        ) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("toroidal");
-
-
-    // Constructors
-
-        //- Construct from origin, rotation and radius
-        toroidalCS
-        (
-            const word& name,
-            const point& origin,
-            const coordinateRotation&,
-            const scalar radius
-        );
-
-        //- Construct from dictionary
-        toroidalCS(const word& name, const dictionary&);
-
-
-    // Member Functions
-
-        //- Return radius
-        scalar radius() const
-        {
-            return radius_;
-        }
-
-        //- Write
-        virtual void write(Ostream&) const;
-
-        //- Write dictionary
-        virtual void writeDict(Ostream&, bool subDict=true) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/meshTools/meshStructure/meshStructure.C

@@ -0,0 +1,387 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshStructure.H"
+#include "FaceCellWave.H"
+#include "topoDistanceData.H"
+#include "pointTopoDistanceData.H"
+#include "PointEdgeWave.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(meshStructure, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::meshStructure::isStructuredCell
+(
+    const polyMesh& mesh,
+    const label layerI,
+    const label cellI
+) const
+{
+    const cell& cFaces = mesh.cells()[cellI];
+
+    // Count number of side faces
+    label nSide = 0;
+    forAll(cFaces, i)
+    {
+        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
+        {
+            nSide++;
+        }
+    }
+
+    if (nSide != cFaces.size()-2)
+    {
+        return false;
+    }
+
+    // Check that side faces have correct point layers
+    forAll(cFaces, i)
+    {
+        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
+        {
+            const face& f = mesh.faces()[cFaces[i]];
+
+            label nLayer = 0;
+            label nLayerPlus1 = 0;
+            forAll(f, fp)
+            {
+                label pointI = f[fp];
+                if (pointLayer_[pointI] == layerI)
+                {
+                    nLayer++;
+                }
+                else if (pointLayer_[pointI] == layerI+1)
+                {
+                    nLayerPlus1++;
+                }
+            }
+
+            if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
+            {
+                return false;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+void Foam::meshStructure::correct
+(
+    const polyMesh& mesh,
+    const uindirectPrimitivePatch& pp
+)
+{
+    // Field on cells and faces.
+    List<topoDistanceData> cellData(mesh.nCells());
+    List<topoDistanceData> faceData(mesh.nFaces());
+
+    {
+        if (debug)
+        {
+            Info<< typeName << " : seeding "
+                << returnReduce(pp.size(), sumOp<label>()) << " patch faces"
+                << nl << endl;
+        }
+
+
+        // Start of changes
+        labelList patchFaces(pp.size());
+        List<topoDistanceData> patchData(pp.size());
+        forAll(pp, patchFaceI)
+        {
+            patchFaces[patchFaceI] = pp.addressing()[patchFaceI];
+            patchData[patchFaceI] = topoDistanceData(patchFaceI, 0);
+        }
+
+
+        // Propagate information inwards
+        FaceCellWave<topoDistanceData> distanceCalc
+        (
+            mesh,
+            patchFaces,
+            patchData,
+            faceData,
+            cellData,
+            mesh.globalData().nTotalCells()+1
+        );
+
+
+        // Determine cells from face-cell-walk
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        cellToPatchFaceAddressing_.setSize(mesh.nCells());
+        cellLayer_.setSize(mesh.nCells());
+        forAll(cellToPatchFaceAddressing_, cellI)
+        {
+            cellToPatchFaceAddressing_[cellI] = cellData[cellI].data();
+            cellLayer_[cellI] = cellData[cellI].distance();
+        }
+
+
+
+        // Determine faces from face-cell-walk
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        faceToPatchFaceAddressing_.setSize(mesh.nFaces());
+        faceToPatchEdgeAddressing_.setSize(mesh.nFaces());
+        faceToPatchEdgeAddressing_ = labelMin;
+        faceLayer_.setSize(mesh.nFaces());
+
+        forAll(faceToPatchFaceAddressing_, faceI)
+        {
+            label own = mesh.faceOwner()[faceI];
+            label patchFaceI = faceData[faceI].data();
+            label patchDist = faceData[faceI].distance();
+
+            if (mesh.isInternalFace(faceI))
+            {
+                label nei = mesh.faceNeighbour()[faceI];
+
+                if (cellData[own].distance() == cellData[nei].distance())
+                {
+                    // side face
+                    faceToPatchFaceAddressing_[faceI] = 0;
+                    faceLayer_[faceI] = cellData[own].distance();
+                }
+                else if (cellData[own].distance() < cellData[nei].distance())
+                {
+                    // unturned face
+                    faceToPatchFaceAddressing_[faceI] = patchFaceI+1;
+                    faceToPatchEdgeAddressing_[faceI] = -1;
+                    faceLayer_[faceI] = patchDist;
+                }
+                else
+                {
+                    // turned face
+                    faceToPatchFaceAddressing_[faceI] = -(patchFaceI+1);
+                    faceToPatchEdgeAddressing_[faceI] = -1;
+                    faceLayer_[faceI] = patchDist;
+                }
+            }
+            else if (patchDist == cellData[own].distance())
+            {
+                // starting face
+                faceToPatchFaceAddressing_[faceI] = -(patchFaceI+1);
+                faceToPatchEdgeAddressing_[faceI] = -1;
+                faceLayer_[faceI] = patchDist;
+            }
+            else
+            {
+                // unturned face or side face. Cannot be determined until
+                // we determine the point layers. Problem is that both are
+                // the same number of steps away from the initial seed face.
+            }
+        }
+    }
+
+
+    // Determine points from separate walk on point-edge
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        pointToPatchPointAddressing_.setSize(mesh.nPoints());
+        pointLayer_.setSize(mesh.nPoints());
+
+        if (debug)
+        {
+            Info<< typeName << " : seeding "
+                << returnReduce(pp.nPoints(), sumOp<label>()) << " patch points"
+                << nl << endl;
+        }
+
+        // Field on edges and points.
+        List<pointTopoDistanceData> edgeData(mesh.nEdges());
+        List<pointTopoDistanceData> pointData(mesh.nPoints());
+
+        // Start of changes
+        labelList patchPoints(pp.nPoints());
+        List<pointTopoDistanceData> patchData(pp.nPoints());
+        forAll(pp.meshPoints(), patchPointI)
+        {
+            patchPoints[patchPointI] = pp.meshPoints()[patchPointI];
+            patchData[patchPointI] = pointTopoDistanceData(patchPointI, 0);
+        }
+
+
+        // Walk
+        PointEdgeWave<pointTopoDistanceData> distanceCalc
+        (
+            mesh,
+            patchPoints,
+            patchData,
+
+            pointData,
+            edgeData,
+            mesh.globalData().nTotalPoints()  // max iterations
+        );
+
+        forAll(pointData, pointI)
+        {
+            pointToPatchPointAddressing_[pointI] = pointData[pointI].data();
+            pointLayer_[pointI] = pointData[pointI].distance();
+        }
+
+
+        // Derive from originating patch points what the patch edges were.
+        EdgeMap<label> pointsToEdge(pp.nEdges());
+        forAll(pp.edges(), edgeI)
+        {
+            pointsToEdge.insert(pp.edges()[edgeI], edgeI);
+        }
+
+        // Look up on faces
+        forAll(faceToPatchEdgeAddressing_, faceI)
+        {
+            if (faceToPatchEdgeAddressing_[faceI] == labelMin)
+            {
+                // Face not yet done. Check if all points on same level
+                // or if not see what edge it originates from
+
+                const face& f = mesh.faces()[faceI];
+
+                label levelI = pointLayer_[f[0]];
+                for (label fp = 1; fp < f.size(); fp++)
+                {
+                    if (pointLayer_[f[fp]] != levelI)
+                    {
+                        levelI = -1;
+                        break;
+                    }
+                }
+
+                if (levelI != -1)
+                {
+                    // All same level
+                    //Pout<< "Horizontal boundary face " << faceI
+                    //    << " at:" << mesh.faceCentres()[faceI]
+                    //    << " data:" << faceData[faceI]
+                    //    << " pointDatas:"
+                    //    << UIndirectList<pointTopoDistanceData>(pointData, f)
+                    //    << endl;
+
+                    label patchFaceI = faceData[faceI].data();
+                    label patchDist = faceData[faceI].distance();
+
+                    faceToPatchEdgeAddressing_[faceI] = -1;
+                    faceToPatchFaceAddressing_[faceI] = patchFaceI+1;
+                    faceLayer_[faceI] = patchDist;
+                }
+                else
+                {
+                    // Points of face on different levels
+
+                    // See if there is any edge
+                    forAll(f, fp)
+                    {
+                        label pointI = f[fp];
+                        label nextPointI = f.nextLabel(fp);
+
+                        EdgeMap<label>::const_iterator fnd = pointsToEdge.find
+                        (
+                            edge
+                            (
+                                pointData[pointI].data(),
+                                pointData[nextPointI].data()
+                            )
+                        );
+                        if (fnd != pointsToEdge.end())
+                        {
+                            faceToPatchEdgeAddressing_[faceI] = fnd();
+                            faceToPatchFaceAddressing_[faceI] = 0;
+                            label own = mesh.faceOwner()[faceI];
+                            faceLayer_[faceI] = cellData[own].distance();
+
+                            // Note: could test whether the other edges on the
+                            // face are consistent
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // Use maps to find out mesh structure.
+    {
+        label nLayers = gMax(cellLayer_)+1;
+        labelListList layerToCells(invertOneToMany(nLayers, cellLayer_));
+
+        structured_ = true;
+        forAll(layerToCells, layerI)
+        {
+            const labelList& lCells = layerToCells[layerI];
+
+            forAll(lCells, lCellI)
+            {
+                label cellI = lCells[lCellI];
+
+                structured_ = isStructuredCell
+                (
+                    mesh,
+                    layerI,
+                    cellI
+                );
+
+                if (!structured_)
+                {
+                    break;
+                }
+            }
+
+            if (!structured_)
+            {
+                break;
+            }
+        }
+
+        reduce(structured_, andOp<bool>());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshStructure::meshStructure
+(
+    const polyMesh& mesh,
+    const uindirectPrimitivePatch& pp
+)
+{
+    correct(mesh, pp);
+}
+
+
+// ************************************************************************* //

src/meshTools/meshStructure/meshStructure.H

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    Foam::meshStructure
+
+Description
+    Detect extruded mesh structure given a set of patch faces.
+
+SourceFiles
+    meshStructure.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshStructure_H
+#define meshStructure_H
+
+#include "labelList.H"
+#include "uindirectPrimitivePatch.H"
+#include "className.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+
+
+/*---------------------------------------------------------------------------*\
+                         Class meshStructure Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshStructure
+{
+    // Private data
+
+        //- Cell to patch face
+        labelList cellToPatchFaceAddressing_;
+
+        //- Cell to layer
+        labelList cellLayer_;
+
+        //- Face to patch face
+        labelList faceToPatchFaceAddressing_;
+
+        //- Face to patch edge
+        labelList faceToPatchEdgeAddressing_;
+
+        //- Face to layer
+        labelList faceLayer_;
+
+        //- Point to patch point
+        labelList pointToPatchPointAddressing_;
+
+        //- Point to layer
+        labelList pointLayer_;
+
+        //- Is mesh structured?
+        bool structured_;
+
+
+   // Private Member Functions
+
+        //- Is cell structured
+        bool isStructuredCell
+        (
+            const polyMesh& mesh,
+            const label layerI,
+            const label cellI
+        ) const;
+
+        //- Calculate all maps.
+        void correct
+        (
+            const polyMesh& mesh,
+            const uindirectPrimitivePatch& pp
+        );
+
+
+public:
+
+    // Declare name of the class and its debug switch
+    ClassName("meshStructure");
+
+    // Constructors
+
+        //- Construct null
+        meshStructure(const polyMesh& mesh, const uindirectPrimitivePatch&);
+
+
+    // Member Functions
+
+        //- Is mesh structured?
+        inline bool structured() const;
+
+        //- Cell to patch face
+        inline const labelList& cellToPatchFaceAddressing() const;
+
+        //- Cell to patch face
+        inline labelList& cellToPatchFaceAddressing();
+
+        //- Cell to layer
+        inline const labelList& cellLayer() const;
+
+        //- Cell to layer
+        inline labelList& cellLayer();
+
+        //- Face to patch face
+        inline const labelList& faceToPatchFaceAddressing() const;
+
+        //- Face to patch face
+        inline labelList& faceToPatchFaceAddressing();
+
+        //- Face to patch edge
+        inline const labelList& faceToPatchEdgeAddressing() const;
+
+        //- Face to patch edge
+        inline labelList& faceToPatchEdgeAddressing();
+
+        //- Face to layer
+        inline const labelList& faceLayer() const;
+
+        //- Face to layer
+        inline labelList& faceLayer();
+
+        //- Point to patch point
+        inline const labelList& pointToPatchPointAddressing() const;
+
+        //- Point to patch point
+        inline labelList& pointToPatchPointAddressing();
+
+        //- Point to layer
+        inline const labelList& pointLayer() const;
+
+        //- Point to layer
+        inline labelList& pointLayer();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "meshStructureI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/meshStructure/meshStructureI.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshStructure.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::meshStructure::structured() const
+{
+    return structured_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::cellToPatchFaceAddressing() const
+{
+    return cellToPatchFaceAddressing_;
+}
+
+
+Foam::labelList& Foam::meshStructure::cellToPatchFaceAddressing()
+{
+    return cellToPatchFaceAddressing_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::cellLayer() const
+{
+    return cellLayer_;
+}
+
+
+Foam::labelList& Foam::meshStructure::cellLayer()
+{
+    return cellLayer_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::faceToPatchFaceAddressing() const
+{
+    return faceToPatchFaceAddressing_;
+}
+
+
+Foam::labelList& Foam::meshStructure::faceToPatchFaceAddressing()
+{
+    return faceToPatchFaceAddressing_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::faceToPatchEdgeAddressing() const
+{
+    return faceToPatchEdgeAddressing_;
+}
+
+
+Foam::labelList& Foam::meshStructure::faceToPatchEdgeAddressing()
+{
+    return faceToPatchEdgeAddressing_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::faceLayer() const
+{
+    return faceLayer_;
+}
+
+
+Foam::labelList& Foam::meshStructure::faceLayer()
+{
+    return faceLayer_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::pointToPatchPointAddressing() const
+{
+    return pointToPatchPointAddressing_;
+}
+
+
+Foam::labelList& Foam::meshStructure::pointToPatchPointAddressing()
+{
+    return pointToPatchPointAddressing_;
+}
+
+
+const Foam::labelList& Foam::meshStructure::pointLayer() const
+{
+    return pointLayer_;
+}
+
+
+Foam::labelList& Foam::meshStructure::pointLayer()
+{
+    return pointLayer_;
+}
+
+
+// ************************************************************************* //

src/meshTools/meshStructure/pointTopoDistanceData.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointTopoDistanceData.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::pointTopoDistanceData& wDist
+)
+{
+    return os << wDist.data_ << token::SPACE << wDist.distance_;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::pointTopoDistanceData& wDist
+)
+{
+    return is >> wDist.data_ >> wDist.distance_;
+}
+
+
+// ************************************************************************* //

src/meshTools/meshStructure/pointTopoDistanceData.H

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    Foam::pointTopoDistanceData
+
+Description
+    For use with PointEdgeWave. Determines topological distance to
+    starting points
+
+SourceFiles
+    pointTopoDistanceDataI.H
+    pointTopoDistanceData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointTopoDistanceData_H
+#define pointTopoDistanceData_H
+
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class pointTopoDistanceData Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointTopoDistanceData
+{
+    // Private data
+
+        //- Starting data
+        label data_;
+
+        //- Distance
+        label distance_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointTopoDistanceData();
+
+        //- Construct from count
+        inline pointTopoDistanceData
+        (
+            const label data,
+            const label distance
+        );
+
+
+    // Member Functions
+
+        // Access
+
+
+            inline label data() const
+            {
+                return data_;
+            }
+            inline label distance() const
+            {
+                return distance_;
+            }
+
+
+        // Needed by PointEdgeWave
+
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            template<class TrackingData>
+            inline bool sameGeometry
+            (
+                const pointTopoDistanceData&,
+                const scalar tol,
+                TrackingData& td
+            ) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            template<class TrackingData>
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            template<class TrackingData>
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Apply rotation matrix to origin
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
+
+            //- Influence of edge on point
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointTopoDistanceData& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointTopoDistanceData& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const pointTopoDistanceData& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of point on edge.
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointTopoDistanceData& pointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const pointTopoDistanceData&, TrackingData&)
+            const;
+
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const pointTopoDistanceData&) const;
+        inline bool operator!=(const pointTopoDistanceData&) const;
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointTopoDistanceData&);
+        friend Istream& operator>>(Istream&, pointTopoDistanceData&);
+};
+
+
+//- Data associated with pointTopoDistanceData type are contiguous
+template<>
+inline bool contiguous<pointTopoDistanceData>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointTopoDistanceDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/meshStructure/pointTopoDistanceDataI.H

@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointTopoDistanceData::pointTopoDistanceData()
+:
+    data_(-1),
+    distance_(-1)
+{}
+
+
+// Construct from components
+inline Foam::pointTopoDistanceData::pointTopoDistanceData
+(
+    const label data,
+    const label distance
+)
+:
+    data_(data),
+    distance_(distance)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class TrackingData>
+inline bool Foam::pointTopoDistanceData::valid(TrackingData& td) const
+{
+    return distance_ != -1;
+}
+
+
+// No geometric data so never any problem on cyclics
+template <class TrackingData>
+inline bool Foam::pointTopoDistanceData::sameGeometry
+(
+    const pointTopoDistanceData&,
+    const scalar tol,
+    TrackingData& td
+) const
+{
+    return true;
+}
+
+
+// No geometric data.
+template <class TrackingData>
+inline void Foam::pointTopoDistanceData::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord,
+    TrackingData& td
+)
+{}
+
+
+// No geometric data.
+template <class TrackingData>
+inline void Foam::pointTopoDistanceData::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
+{}
+
+
+// No geometric data.
+template <class TrackingData>
+inline void Foam::pointTopoDistanceData::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord,
+    TrackingData& td
+)
+{}
+
+
+// Update this with information from connected edge
+template <class TrackingData>
+inline bool Foam::pointTopoDistanceData::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointTopoDistanceData& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (distance_ == -1)
+    {
+        data_ = edgeInfo.data_;
+        distance_ = edgeInfo.distance_ + 1;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with new information on same point
+template<class TrackingData>
+inline bool Foam::pointTopoDistanceData::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointTopoDistanceData& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (distance_ == -1)
+    {
+        operator=(newPointInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with new information on same point. No extra information.
+template<class TrackingData>
+inline bool Foam::pointTopoDistanceData::updatePoint
+(
+    const pointTopoDistanceData& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (distance_ == -1)
+    {
+        operator=(newPointInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with information from connected point
+template<class TrackingData>
+inline bool Foam::pointTopoDistanceData::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointTopoDistanceData& pointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (distance_ == -1)
+    {
+        operator=(pointInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template <class TrackingData>
+inline bool Foam::pointTopoDistanceData::equal
+(
+    const pointTopoDistanceData& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointTopoDistanceData::operator==
+(
+    const Foam::pointTopoDistanceData& rhs
+) const
+{
+    return data() == rhs.data() && distance() == rhs.distance();
+}
+
+
+inline bool Foam::pointTopoDistanceData::operator!=
+(
+    const Foam::pointTopoDistanceData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/meshTools/meshStructure/topoDistanceData.C

@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/meshTools/meshStructure/topoDistanceData.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License