From 976ad36776844a8b87ab215c52f3704c6096dca7 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Mon, 24 Apr 2017 10:34:05 +0100
Subject: [PATCH] ENH: Initial attempt to track oriented surface fields

---
 .../rhoCentralDyMFoam/rhoCentralDyMFoam.C     |   3 +
 .../rhoCentralFoam/rhoCentralFoam.C           |   2 +
 .../simpleReactingParcelFoam.C                |   1 -
 .../multiphase/interFoam/alphaEqnSubCycle.H   |   2 +
 src/OpenFOAM/Make/files                       |   2 +
 .../DimensionedField/DimensionedField.C       |  42 +-
 .../DimensionedField/DimensionedField.H       |  14 +-
 .../DimensionedFieldFunctionsM.C              |  44 +-
 .../DimensionedField/DimensionedFieldI.H      |  22 +-
 .../DimensionedField/DimensionedFieldIO.C     |  19 +-
 .../GeometricField/GeometricField.C           |   5 +-
 .../GeometricField/GeometricFieldFunctions.C  |  11 +
 .../GeometricField/GeometricFieldFunctionsM.C |   9 +
 src/OpenFOAM/orientedType/orientedType.C      | 498 ++++++++++++++++++
 src/OpenFOAM/orientedType/orientedType.H      | 168 ++++++
 .../gaussLaplacianScheme.C                    |   1 +
 .../snGradSchemes/snGradScheme/snGradScheme.C |   1 +
 .../fvMatrices/fvMatrix/fvMatrix.C            |   2 +
 src/finiteVolume/fvMesh/fvMeshGeometry.C      |   5 +
 .../surfaceInterpolation.C                    |  23 +-
 .../surfaceInterpolationScheme.C              |   5 +
 21 files changed, 836 insertions(+), 43 deletions(-)
 create mode 100644 src/OpenFOAM/orientedType/orientedType.C
 create mode 100644 src/OpenFOAM/orientedType/orientedType.H

diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
index 96413e02017..4a58e5d3fce 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
@@ -120,11 +120,13 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
+cSf_pos.oriented().oriented() = true;
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
+cSf_neg.oriented().oriented() = true;
 
         surfaceScalarField ap
         (
@@ -269,4 +271,5 @@ int main(int argc, char *argv[])
     return 0;
 }
 
+
 // ************************************************************************* //
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 65ab9d2db3e..b212fd30ea8 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@@ -101,11 +101,13 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
+cSf_pos.oriented().oriented() = true;
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
+cSf_neg.oriented().oriented() = true;
 
         surfaceScalarField ap
         (
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
index b3504da91a5..bbec4a41ab9 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
@@ -38,7 +38,6 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "IOporosityModelList.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "simpleControl.H"
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 772de0b97af..dd1e69b2cc6 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -13,6 +13,8 @@ if (nAlphaSubCycles > 1)
         dimensionedScalar("0", rhoPhi.dimensions(), 0)
     );
 
+    rhoPhiSum.oriented().oriented() = true;
+
     tmp<volScalarField> trSubDeltaT;
 
     if (LTS)
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index b5c58d17291..8c674c627b8 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -265,6 +265,8 @@ dimensionedTypes/dimensionedSphericalTensor/dimensionedSphericalTensor.C
 dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.C
 dimensionedTypes/dimensionedTensor/dimensionedTensor.C
 
+orientedType/orientedType.C
+
 matrices/solution/solution.C
 
 scalarMatrices = matrices/scalarMatrices
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
index 166f7941faa..9901226ba23 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
@@ -59,7 +59,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(field),
     mesh_(mesh),
-    dimensions_(dims)
+    dimensions_(dims),
+    oriented_(false)
 {
     if (field.size() && field.size() != GeoMesh::size(mesh))
     {
@@ -84,7 +85,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(GeoMesh::size(mesh)),
     mesh_(mesh),
-    dimensions_(dims)
+    dimensions_(dims),
+    oriented_(false)
 {
     if (checkIOFlags)
     {
@@ -105,7 +107,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(GeoMesh::size(mesh), dt.value()),
     mesh_(mesh),
-    dimensions_(dt.dimensions())
+    dimensions_(dt.dimensions()),
+    oriented_(false)
 {
     if (checkIOFlags)
     {
@@ -123,7 +126,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(df),
     Field<Type>(df),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -137,7 +141,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(df, reuse),
     Field<Type>(df, reuse),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -150,7 +155,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(df(), true),
     Field<Type>(df),
     mesh_(df->mesh_),
-    dimensions_(df->dimensions_)
+    dimensions_(df->dimensions_),
+    oriented_(df->oriented_)
 {}
 
 
@@ -168,7 +174,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
         tdf.isTmp()
     ),
     mesh_(tdf().mesh_),
-    dimensions_(tdf().dimensions_)
+    dimensions_(tdf().dimensions_),
+    oriented_(tdf().oriented_)
 {
     tdf.clear();
 }
@@ -185,7 +192,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(df),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -200,7 +208,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io, df),
     Field<Type>(df, reuse),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -214,7 +223,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(newName, df, newName == df.name()),
     Field<Type>(df),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -229,7 +239,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(newName, df, true),
     Field<Type>(df, reuse),
     mesh_(df.mesh_),
-    dimensions_(df.dimensions_)
+    dimensions_(df.dimensions_),
+    oriented_(df.oriented_)
 {}
 
 
@@ -243,7 +254,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(newName, df, true),
     Field<Type>(df),
     mesh_(df->mesh_),
-    dimensions_(df->dimensions_)
+    dimensions_(df->dimensions_),
+    oriented_(df->oriented_)
 {}
 
 
@@ -262,7 +274,8 @@ DimensionedField<Type, GeoMesh>::DimensionedField
         tdf.isTmp()
     ),
     mesh_(tdf().mesh_),
-    dimensions_(tdf().dimensions_)
+    dimensions_(tdf().dimensions_),
+    oriented_(tdf().oriented_)
 {
     tdf.clear();
 }
@@ -437,6 +450,7 @@ void DimensionedField<Type, GeoMesh>::operator=
     checkField(*this, df, "=");
 
     dimensions_ = df.dimensions();
+    oriented_ = df.oriented();
     Field<Type>::operator=(df);
 }
 
@@ -460,6 +474,7 @@ void DimensionedField<Type, GeoMesh>::operator=
     checkField(*this, df, "=");
 
     dimensions_ = df.dimensions();
+    oriented_ = df.oriented();
     this->transfer(const_cast<DimensionedField<Type, GeoMesh>&>(df));
     tdf.clear();
 }
@@ -487,6 +502,7 @@ void DimensionedField<Type, GeoMesh>::operator op                              \
     checkField(*this, df, #op);                                                \
                                                                                \
     dimensions_ op df.dimensions();                                            \
+    oriented_ op df.oriented();                                                \
     Field<Type>::operator op(df);                                              \
 }                                                                              \
                                                                                \
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
index a302e541710..ff894271c2a 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,7 @@ SourceFiles
 #include "regIOobject.H"
 #include "Field.H"
 #include "dimensionedType.H"
+#include "orientedType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +66,7 @@ template<class Type, class GeoMesh> Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                           Class DimensionedField Declaration
+                      Class DimensionedField Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type, class GeoMesh>
@@ -99,6 +100,9 @@ private:
         //- Dimension set for this field
         dimensionSet dimensions_;
 
+        //- Oriented flag
+        orientedType oriented_;
+
 
     // Private Member Functions
 
@@ -262,6 +266,12 @@ public:
         //- Return non-const access to dimensions
         inline dimensionSet& dimensions();
 
+        //- Return oriented flag
+        inline const orientedType& oriented() const;
+
+        //- Return non-const access to the oriented flag
+        inline orientedType& oriented();
+
         inline const Field<Type>& field() const;
 
         inline Field<Type>& field();
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctionsM.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctionsM.C
index 6da09d945c4..8e69d617454 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctionsM.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctionsM.C
@@ -52,6 +52,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field());                                     \
                                                                                \
+    tRes.ref().oriented() = Dfunc(df1.oriented());                             \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -75,6 +77,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field());                                     \
                                                                                \
+    tRes.ref().oriented() = Dfunc(df1.oriented());                             \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -108,6 +112,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field());                             \
                                                                                \
+    tRes.ref().oriented() = Dfunc(df1.oriented());                             \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -131,6 +137,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field());                             \
                                                                                \
+    tRes.ref().oriented() = Dfunc(df1.oriented());                             \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -165,6 +173,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = Func(df1.oriented(), df2.oriented());              \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -189,6 +199,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = Func(df1.oriented(), df2.oriented());              \
+                                                                               \
     tdf2.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -215,6 +227,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = Func(df1.oriented(), df2.oriented());              \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -244,6 +258,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = Func(df1.oriented(), df2.oriented());              \
+                                                                               \
     tdf1.clear();                                                              \
     tdf2.clear();                                                              \
                                                                                \
@@ -279,6 +295,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), dt1.value(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = df2.oriented();                                    \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -314,6 +332,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), dt1.value(), df2.field());                        \
                                                                                \
+    tRes.ref().oriented() = df2.oriented();                                    \
+                                                                               \
     tdf2.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -356,6 +376,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), dt2.value());                        \
                                                                                \
+    tRes.ref().oriented() = df1.oriented();                                    \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -391,6 +413,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> Func                                \
                                                                                \
     Func(tRes.ref().field(), df1.field(), dt2.value());                        \
                                                                                \
+    tRes.ref().oriented() = df1.oriented();                                    \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -440,6 +464,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field(), df2.field());                \
                                                                                \
+    tRes.ref().oriented() = df1.oriented() Op df2.oriented();                  \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -464,6 +490,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field(), df2.field());                \
                                                                                \
+    tRes.ref().oriented() = df1.oriented() Op df2.oriented();                  \
+                                                                               \
     tdf2.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -490,6 +518,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field(), df2.field());                \
                                                                                \
+    tRes.ref().oriented() = df1.oriented() Op df2.oriented();                  \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -519,6 +549,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field(), df2.field());                \
                                                                                \
+    tRes.ref().oriented() = df1.oriented() Op df2.oriented();                  \
+                                                                               \
     tdf1.clear();                                                              \
     tdf2.clear();                                                              \
                                                                                \
@@ -528,7 +560,7 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define BINARY_TYPE_OPERATOR_SF(ReturnType, Type1, Type2, Op, OpName, OpFunc) \
+#define BINARY_TYPE_OPERATOR_SF(ReturnType, Type1, Type2, Op, OpName, OpFunc)  \
                                                                                \
 TEMPLATE                                                                       \
 tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
@@ -552,6 +584,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
         )                                                                      \
     );                                                                         \
                                                                                \
+    tRes.ref().oriented() = df2.oriented();                                    \
+                                                                               \
     Foam::OpFunc(tRes.ref().field(), dt1.value(), df2.field());                \
                                                                                \
     return tRes;                                                               \
@@ -589,6 +623,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), dt1.value(), tdf2().field());             \
                                                                                \
+    tRes.ref().oriented() = df2.oriented();                                    \
+                                                                               \
     tdf2.clear();                                                              \
                                                                                \
     return tRes;                                                               \
@@ -605,7 +641,7 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
 }
 
 
-#define BINARY_TYPE_OPERATOR_FS(ReturnType, Type1, Type2, Op, OpName, OpFunc) \
+#define BINARY_TYPE_OPERATOR_FS(ReturnType, Type1, Type2, Op, OpName, OpFunc)  \
                                                                                \
 TEMPLATE                                                                       \
 tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
@@ -631,6 +667,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), df1.field(), dt2.value());                \
                                                                                \
+    tRes.ref().oriented() = df1.oriented();                                    \
+                                                                               \
     return tRes;                                                               \
 }                                                                              \
                                                                                \
@@ -666,6 +704,8 @@ tmp<DimensionedField<ReturnType, GeoMesh>> operator Op                         \
                                                                                \
     Foam::OpFunc(tRes.ref().field(), tdf1().field(), dt2.value());             \
                                                                                \
+    tRes.ref().oriented() = df1.oriented();                                    \
+                                                                               \
     tdf1.clear();                                                              \
                                                                                \
     return tRes;                                                               \
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldI.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldI.H
index 64cfd315895..0de4d18faa0 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldI.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,14 +48,29 @@ Foam::DimensionedField<Type, GeoMesh>::dimensions() const
     return dimensions_;
 }
 
+
 template<class Type, class GeoMesh>
-inline Foam::dimensionSet&
-Foam::DimensionedField<Type, GeoMesh>::dimensions()
+inline Foam::dimensionSet& Foam::DimensionedField<Type, GeoMesh>::dimensions()
 {
     return dimensions_;
 }
 
 
+template<class Type, class GeoMesh>
+inline const Foam::orientedType&
+Foam::DimensionedField<Type, GeoMesh>::oriented() const
+{
+    return oriented_;
+}
+
+
+template<class Type, class GeoMesh>
+inline Foam::orientedType& Foam::DimensionedField<Type, GeoMesh>::oriented()
+{
+    return oriented_;
+}
+
+
 template<class Type, class GeoMesh>
 inline const Foam::Field<Type>&
 Foam::DimensionedField<Type, GeoMesh>::field() const
@@ -63,6 +78,7 @@ Foam::DimensionedField<Type, GeoMesh>::field() const
     return *this;
 }
 
+
 template<class Type, class GeoMesh>
 inline Foam::Field<Type>&
 Foam::DimensionedField<Type, GeoMesh>::field()
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
index 74da71025c3..5fe2ed249f9 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C
@@ -37,6 +37,7 @@ void Foam::DimensionedField<Type, GeoMesh>::readField
 )
 {
     dimensions_.reset(dimensionSet(fieldDict.lookup("dimensions")));
+    fieldDict.template readIfPresent<bool>("oriented", oriented_.oriented());
 
     Field<Type> f(fieldDictEntry, fieldDict, GeoMesh::size(mesh_));
     this->transfer(f);
@@ -74,7 +75,8 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(0),
     mesh_(mesh),
-    dimensions_(dimless)
+    dimensions_(dimless),
+    oriented_(false)
 {
     readField(dictionary(readStream(typeName)), fieldDictEntry);
 }
@@ -92,7 +94,8 @@ Foam::DimensionedField<Type, GeoMesh>::DimensionedField
     regIOobject(io),
     Field<Type>(0),
     mesh_(mesh),
-    dimensions_(dimless)
+    dimensions_(dimless),
+    oriented_(false)
 {
     readField(fieldDict, fieldDictEntry);
 }
@@ -107,17 +110,15 @@ bool Foam::DimensionedField<Type, GeoMesh>::writeData
     const word& fieldDictEntry
 ) const
 {
-    os.writeKeyword("dimensions") << dimensions() << token::END_STATEMENT
-        << nl << nl;
+    os.writeEntry("dimensions", dimensions());
+    if (oriented_.oriented()) os.writeEntry("oriented", oriented());
+
+    os<< nl << nl;
 
     Field<Type>::writeEntry(fieldDictEntry, os);
 
     // Check state of Ostream
-    os.check
-    (
-        "bool DimensionedField<Type, GeoMesh>::writeData"
-        "(Ostream& os, const word& fieldDictEntry) const"
-    );
+    os.check(FUNCTION_NAME);
 
     return (os.good());
 }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index e7d6ea91402..6cc6dc61293 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -1175,6 +1175,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
     // Only assign field contents not ID
 
     this->dimensions() = gf.dimensions();
+    this->oriented() = gf.oriented();
 
     // Transfer the storage from the tmp
     primitiveFieldRef().transfer
@@ -1239,7 +1240,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op              \
 {                                                                              \
     checkField(*this, gf, #op);                                                \
                                                                                \
-    ref() op gf();            \
+    ref() op gf();                                                             \
     boundaryFieldRef() op gf.boundaryField();                                  \
 }                                                                              \
                                                                                \
@@ -1259,7 +1260,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op              \
     const dimensioned<TYPE>& dt                                                \
 )                                                                              \
 {                                                                              \
-    ref() op dt;                                       \
+    ref() op dt;                                                               \
     boundaryFieldRef() op dt.value();                                          \
 }
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index 49a7ce2d622..4113b1d752b 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
@@ -51,6 +51,7 @@ void component
 {
     component(gcf.primitiveFieldRef(), gf.primitiveField(), d);
     component(gcf.boundaryFieldRef(), gf.boundaryField(), d);
+    gcf.oriented() = gf.oriented();
 }
 
 
@@ -63,6 +64,7 @@ void T
 {
     T(gf.primitiveFieldRef(), gf1.primitiveField());
     T(gf.boundaryFieldRef(), gf1.boundaryField());
+    gf.oriented() = gf1.oriented();
 }
 
 
@@ -81,6 +83,7 @@ void pow
 {
     pow(gf.primitiveFieldRef(), gf1.primitiveField(), r);
     pow(gf.boundaryFieldRef(), gf1.boundaryField(), r);
+    gf.oriented() = pow(gf1.oriented(), r);
 }
 
 template
@@ -175,6 +178,7 @@ void sqr
 {
     sqr(gf.primitiveFieldRef(), gf1.primitiveField());
     sqr(gf.boundaryFieldRef(), gf1.boundaryField());
+    gf.oriented() = sqr(gf1.oriented());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -263,6 +267,7 @@ void magSqr
 {
     magSqr(gsf.primitiveFieldRef(), gf.primitiveField());
     magSqr(gsf.boundaryFieldRef(), gf.boundaryField());
+    gsf.oriented() = magSqr(gf.oriented());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -335,6 +340,7 @@ void mag
 {
     mag(gsf.primitiveFieldRef(), gf.primitiveField());
     mag(gsf.boundaryFieldRef(), gf.boundaryField());
+    gsf.oriented() = mag(gf.oriented());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -412,6 +418,7 @@ void cmptAv
 {
     cmptAv(gcf.primitiveFieldRef(), gf.primitiveField());
     cmptAv(gcf.boundaryFieldRef(), gf.boundaryField());
+    gcf.oriented() = cmptAv(gf.oriented());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -611,6 +618,8 @@ void opFunc                                                                    \
         gf1.boundaryField(),                                                   \
         gf2.boundaryField()                                                    \
     );                                                                         \
+                                                                               \
+    gf.oriented() = gf1.oriented() op gf2.oriented();                          \
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -757,6 +766,7 @@ void opFunc                                                                    \
 {                                                                              \
     Foam::opFunc(gf.primitiveFieldRef(), gf1.primitiveField(), dvs.value());   \
     Foam::opFunc(gf.boundaryFieldRef(), gf1.boundaryField(), dvs.value());     \
+    gf.oriented() = gf1.oriented();                                            \
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -870,6 +880,7 @@ void opFunc                                                                    \
 {                                                                              \
     Foam::opFunc(gf.primitiveFieldRef(), dvs.value(), gf1.primitiveField());   \
     Foam::opFunc(gf.boundaryFieldRef(), dvs.value(), gf1.boundaryField());     \
+    gf.oriented() = gf1.oriented();                                            \
 }                                                                              \
                                                                                \
 template                                                                       \
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
index 377339dae97..9ff2ed8639c 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
@@ -43,6 +43,7 @@ void Func                                                                      \
 {                                                                              \
     Foam::Func(res.primitiveFieldRef(), gf1.primitiveField());                 \
     Foam::Func(res.boundaryFieldRef(), gf1.boundaryField());                   \
+    res.oriented() = gf1.oriented();                                           \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -112,6 +113,7 @@ void OpFunc                                                                    \
 {                                                                              \
     Foam::OpFunc(res.primitiveFieldRef(), gf1.primitiveField());               \
     Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField());                 \
+    res.oriented() = gf1.oriented();                                           \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -192,6 +194,7 @@ void Func                                                                      \
         gf1.boundaryField(),                                                   \
         gf2.boundaryField()                                                    \
     );                                                                         \
+    res.oriented() = Func(gf1.oriented(), gf2.oriented());                     \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -320,6 +323,7 @@ void Func                                                                      \
 {                                                                              \
     Foam::Func(res.primitiveFieldRef(), dt1.value(), gf2.primitiveField());    \
     Foam::Func(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());      \
+    res.oriented() = gf2.oriented();                                           \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -411,6 +415,7 @@ void Func                                                                      \
 {                                                                              \
     Foam::Func(res.primitiveFieldRef(), gf1.primitiveField(), dt2.value());    \
     Foam::Func(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());      \
+    res.oriented() = gf1.oriented();                                           \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -511,6 +516,7 @@ void OpFunc                                                                    \
     (res.primitiveFieldRef(), gf1.primitiveField(), gf2.primitiveField());     \
     Foam::OpFunc                                                               \
     (res.boundaryFieldRef(), gf1.boundaryField(), gf2.boundaryField());        \
+    res.oriented() = gf1.oriented() Op gf2.oriented();                         \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -639,6 +645,8 @@ void OpFunc                                                                    \
 {                                                                              \
     Foam::OpFunc(res.primitiveFieldRef(), dt1.value(), gf2.primitiveField());  \
     Foam::OpFunc(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());    \
+    res.oriented() = gf2.oriented();                                           \
+                                                                               \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -730,6 +738,7 @@ void OpFunc                                                                    \
 {                                                                              \
     Foam::OpFunc(res.primitiveFieldRef(), gf1.primitiveField(), dt2.value());  \
     Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());    \
+    res.oriented() = gf1.oriented();                                           \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
diff --git a/src/OpenFOAM/orientedType/orientedType.C b/src/OpenFOAM/orientedType/orientedType.C
new file mode 100644
index 00000000000..08198c0f049
--- /dev/null
+++ b/src/OpenFOAM/orientedType/orientedType.C
@@ -0,0 +1,498 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "orientedType.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::orientedType::orientedType()
+:
+    oriented_(false)
+{}
+
+
+Foam::orientedType::orientedType(const orientedType& of)
+:
+    oriented_(of.oriented_)
+{}
+
+
+Foam::orientedType::orientedType(const bool oriented)
+:
+    oriented_(oriented)
+{}
+
+
+Foam::orientedType::orientedType(Istream& is)
+:
+    oriented_(readBool(is))
+{
+    is.check(FUNCTION_NAME);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool& Foam::orientedType::oriented()
+{
+    return oriented_;
+}
+
+
+bool Foam::orientedType::oriented() const
+{
+    return oriented_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::orientedType::operator=(const orientedType& of)
+{
+    oriented_ = of.oriented();
+}
+
+
+void Foam::orientedType::operator+=(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    if (oriented_ != of.oriented())
+    {
+        FatalErrorInFunction
+            << "Operator += is undefined for oriented and unoriented types. "
+            << "oriented:" << oriented_ << ", of:" << of.oriented()
+            << abort(FatalError);
+    }
+
+    // No change to oriented_ flag
+}
+
+
+void Foam::orientedType::operator-=(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    if (oriented_ != of.oriented())
+    {
+        FatalErrorInFunction
+            << "Operator -= is undefined for oriented and unoriented types. "
+            << "oriented:" << oriented_ << ", of:" << of.oriented()
+            << abort(FatalError);
+    }
+
+    // No change to oriented_ flag
+}
+
+
+void Foam::orientedType::operator*=(const orientedType& of)
+{
+    oriented_ = oriented_ ^ of.oriented();
+}
+
+
+void Foam::orientedType::operator/=(const orientedType& of)
+{
+    oriented_ = oriented_ ^ of.oriented();
+}
+
+
+void Foam::orientedType::operator*=(const scalar s)
+{
+//InfoInFunction << "oriented_: " << oriented_ << endl;
+    // No change to oriented_ flag
+}
+
+
+void Foam::orientedType::operator/=(const scalar s)
+{
+//InfoInFunction << "oriented_: " << oriented_ << endl;
+    // No change to oriented_ flag
+}
+
+
+bool Foam::orientedType::operator()() const
+{
+//InfoInFunction << "oriented_: " << oriented_ << endl;
+    return oriented_;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+Foam::orientedType Foam::max(const orientedType& of1, const orientedType& of2)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    if (of1.oriented() != of2.oriented())
+    {
+        FatalErrorInFunction
+            << "max is undefined for oriented and unoriented types. "
+            << "of1:" << of1.oriented() << ", of2:" << of2.oriented()
+            << abort(FatalError);
+    }
+
+    return of1;
+}
+
+
+Foam::orientedType Foam::min(const orientedType& of1, const orientedType& of2)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    if (of1.oriented() != of2.oriented())
+    {
+        FatalErrorInFunction
+            << "min is undefined for oriented and unoriented types. "
+            << "of1:" << of1.oriented() << ", of2:" << of2.oriented()
+            << abort(FatalError);
+    }
+
+    return of1;
+}
+
+
+Foam::orientedType Foam::cmptMultiply
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::cmptDivide
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::cmptAv(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType Foam::pow(const orientedType& of, const scalar r)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    // Undefined???
+    // - only defined for integers where:
+    //   - odd powers = oriented_ = yes (if of is oriented)
+    //   - even powers = oriented_ = no
+    return of;
+}
+
+
+Foam::orientedType Foam::sqr(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(false);
+}
+
+
+Foam::orientedType Foam::pow3(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of.oriented());
+}
+
+
+Foam::orientedType Foam::pow4(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(false);
+}
+
+
+Foam::orientedType Foam::pow5(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of.oriented());
+}
+
+
+Foam::orientedType Foam::pow6(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(false);
+}
+
+
+Foam::orientedType Foam::pow025(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of.oriented());
+}
+
+
+Foam::orientedType Foam::sqrt(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType Foam::cbrt(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType Foam::magSqr(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(false);
+}
+
+
+Foam::orientedType  Foam::mag(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(false);
+}
+
+
+Foam::orientedType  Foam::sign(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType  Foam::pos(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType  Foam::neg(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType  Foam::posPart(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType  Foam::negPart(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType  Foam::inv(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType Foam::trans(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+Foam::orientedType Foam::atan2
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    if (of1.oriented() != of2.oriented())
+    {
+        FatalErrorInFunction
+            << "atan2 is undefined for oriented and unoriented types. "
+            << "of1:" << of1.oriented() << ", of2:" << of2.oriented()
+            << abort(FatalError);
+    }
+
+    return of1;
+}
+
+
+Foam::orientedType Foam::transform(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return of;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    is >> of.oriented_;
+
+    is.check(FUNCTION_NAME);
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    os << of.oriented();
+
+    os.check(FUNCTION_NAME);
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+Foam::orientedType Foam::operator+
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    if (of1.oriented() != of2.oriented())
+    {
+        FatalErrorInFunction
+            << "Operator + is undefined for oriented and unoriented types. "
+            << "of1:" << of1.oriented() << ", of2:" << of2.oriented()
+            << abort(FatalError);
+    }
+
+    return orientedType(of1.oriented() || of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator-(const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of);
+}
+
+
+Foam::orientedType Foam::operator-
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    if (of1.oriented() != of2.oriented())
+    {
+        FatalErrorInFunction
+            << "Operator - is undefined for oriented and unoriented types. "
+            << "of1:" << of1.oriented() << ", of2:" << of2.oriented()
+            << abort(FatalError);
+    }
+
+    return orientedType(of1.oriented() || of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator*(const scalar s, const orientedType& of)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of);
+}
+
+
+Foam::orientedType Foam::operator/(const orientedType& of, const scalar s)
+{
+//InfoInFunction << "of:" << of.oriented() << endl;
+    return orientedType(of);
+}
+
+
+Foam::orientedType Foam::operator/
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator*
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator^
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator&
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(of1.oriented() ^ of2.oriented());
+}
+
+
+Foam::orientedType Foam::operator&&
+(
+    const orientedType& of1,
+    const orientedType& of2
+)
+{
+//InfoInFunction << "of1:" << of1.oriented() << ", of2:" << of2.oriented() << endl;
+    return orientedType(false);
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/orientedType/orientedType.H b/src/OpenFOAM/orientedType/orientedType.H
new file mode 100644
index 00000000000..4a26164dee0
--- /dev/null
+++ b/src/OpenFOAM/orientedType/orientedType.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::orientedType
+
+Description
+    Class to determine the 'oriented' status of surface fields
+
+SourceFiles
+    orientedType.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef orientedType_H
+#define orientedType_H
+
+#include "Istream.H"
+#include "Ostream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+class orientedType;
+
+Istream& operator>>(Istream&, orientedType&);
+
+Ostream& operator<<(Ostream&, const orientedType&);
+
+/*---------------------------------------------------------------------------*\
+                           Class dimensioned Declaration
+\*---------------------------------------------------------------------------*/
+
+class orientedType
+{
+    // Private data
+
+        //- Oriented flag
+        bool oriented_;
+
+
+public:
+
+
+    // Constructors
+
+        //- Null constructor - flag initialised to false
+        orientedType();
+
+        //- Copy constructor
+        orientedType(const orientedType& of);
+
+        //- Construct from bool
+        orientedType(const bool oriented);
+
+        //- Construct from Istream
+        orientedType(Istream& is);
+
+
+    // Member functions
+
+        //- Return non-const reference to the oriented flag
+        bool& oriented();
+
+        //- Return const reference to the oriented flag
+        bool oriented() const;
+
+
+    // Member operators
+
+        void operator=(const orientedType& of);
+
+        void operator+=(const orientedType& of);
+        void operator-=(const orientedType& of);
+        void operator*=(const orientedType& of);
+        void operator/=(const orientedType& of);
+        void operator*=(const scalar s);
+        void operator/=(const scalar s);
+        bool operator()() const;
+
+
+    // IOstream operators
+
+        friend Istream& operator>>(Istream& is, orientedType& of);
+
+        friend Ostream& operator<<(Ostream& os, const orientedType& of);
+};
+
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+orientedType max(const orientedType& of1, const orientedType& of2);
+orientedType min(const orientedType& of1, const orientedType& of2);
+orientedType cmptMultiply(const orientedType& of1, const orientedType& of2);
+orientedType cmptDivide(const orientedType& of1, const orientedType& of);
+orientedType cmptAv(const orientedType& of);
+
+
+orientedType pow(const orientedType& of, const scalar r);
+orientedType sqr(const orientedType& of);
+orientedType pow3(const orientedType& of);
+orientedType pow4(const orientedType& of);
+orientedType pow5(const orientedType& of);
+orientedType pow6(const orientedType& of);
+orientedType pow025(const orientedType& of);
+
+
+orientedType sqrt(const orientedType& of);
+orientedType cbrt(const orientedType& of);
+orientedType magSqr(const orientedType& of);
+orientedType mag(const orientedType& of);
+orientedType sign(const orientedType& of);
+orientedType pos(const orientedType& of);
+orientedType neg(const orientedType& of);
+orientedType posPart(const orientedType& of);
+orientedType negPart(const orientedType& of);
+orientedType inv(const orientedType& of);
+
+
+orientedType trans(const orientedType& of);
+orientedType atan2(const orientedType& of1, const orientedType& of2);
+orientedType transform(const orientedType& of);
+
+orientedType operator-(const orientedType& of);
+orientedType operator*(const scalar s, const orientedType& of);
+orientedType operator/(const orientedType& of, const scalar s);
+
+orientedType operator+(const orientedType& of1, const orientedType& of2);
+orientedType operator-(const orientedType& of1, const orientedType& of2);
+orientedType operator/(const orientedType& of1, const orientedType& of2);
+orientedType operator*(const orientedType& of1, const orientedType& of2);
+orientedType operator^(const orientedType& of1, const orientedType& of2);
+orientedType operator&(const orientedType& of1, const orientedType& of2);
+orientedType operator&&(const orientedType& of1, const orientedType& of2);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
index 795fbff9759..5f372b1e841 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
@@ -115,6 +115,7 @@ gaussLaplacianScheme<Type, GType>::gammaSnGradCorr
            *vf.dimensions()*mesh.deltaCoeffs().dimensions()
         )
     );
+    tgammaSnGradCorr.ref().oriented() = SfGammaCorr.oriented();
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
index 2ad75944f44..b6f6fbbef6a 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
@@ -124,6 +124,7 @@ snGradScheme<Type>::snGrad
         )
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tsf.ref();
+    ssf.oriented().oriented() = true;
 
     // set reference to difference factors array
     const scalarField& deltaCoeffs = tdeltaCoeffs();
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 4c2ffb5a1f1..bb87c6d7f75 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -892,6 +892,8 @@ flux() const
     GeometricField<Type, fvsPatchField, surfaceMesh>& fieldFlux =
         tfieldFlux.ref();
 
+    fieldFlux.oriented().oriented() = true;
+
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
         fieldFlux.primitiveFieldRef().replace
diff --git a/src/finiteVolume/fvMesh/fvMeshGeometry.C b/src/finiteVolume/fvMesh/fvMeshGeometry.C
index a05569849f2..6d2976636ad 100644
--- a/src/finiteVolume/fvMesh/fvMeshGeometry.C
+++ b/src/finiteVolume/fvMesh/fvMeshGeometry.C
@@ -67,6 +67,8 @@ void Foam::fvMesh::makeSf() const
         dimArea,
         faceAreas()
     );
+
+    SfPtr_->oriented().oriented() = true;
 }
 
 
@@ -400,6 +402,7 @@ Foam::tmp<Foam::surfaceVectorField> Foam::fvMesh::delta() const
         )
     );
     surfaceVectorField& delta = tdelta.ref();
+    delta.oriented().oriented() = true;
 
     const volVectorField& C = this->C();
     const labelUList& owner = this->owner();
@@ -438,6 +441,8 @@ const Foam::surfaceScalarField& Foam::fvMesh::phi() const
         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
     }
 
+    phiPtr_->oriented().oriented() = true;
+
     return *phiPtr_;
 }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index 6e8f8bc7d45..f09cc74c339 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,8 +73,7 @@ Foam::surfaceInterpolation::~surfaceInterpolation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::surfaceScalarField&
-Foam::surfaceInterpolation::weights() const
+const Foam::surfaceScalarField& Foam::surfaceInterpolation::weights() const
 {
     if (!weights_)
     {
@@ -85,8 +84,7 @@ Foam::surfaceInterpolation::weights() const
 }
 
 
-const Foam::surfaceScalarField&
-Foam::surfaceInterpolation::deltaCoeffs() const
+const Foam::surfaceScalarField& Foam::surfaceInterpolation::deltaCoeffs() const
 {
     if (!deltaCoeffs_)
     {
@@ -156,11 +154,12 @@ void Foam::surfaceInterpolation::makeWeights() const
         dimless
     );
     surfaceScalarField& weights = *weights_;
+    weights.oriented().oriented() = true;
 
     // Set local references to mesh data
-    // (note that we should not use fvMesh sliced fields at this point yet
-    //  since this causes a loop when generating weighting factors in
-    //  coupledFvPatchField evaluation phase)
+    // Note that we should not use fvMesh sliced fields at this point yet
+    // since this causes a loop when generating weighting factors in
+    // coupledFvPatchField evaluation phase
     const labelUList& owner = mesh_.owner();
     const labelUList& neighbour = mesh_.neighbour();
 
@@ -174,7 +173,7 @@ void Foam::surfaceInterpolation::makeWeights() const
     forAll(owner, facei)
     {
         // Note: mag in the dot-product.
-        // For all valid meshes, the non-orthogonality will be less that
+        // For all valid meshes, the non-orthogonality will be less than
         // 90 deg and the dot-product will be positive.  For invalid
         // meshes (d & s <= 0), this will stabilise the calculation
         // but the result will be poor.
@@ -183,8 +182,7 @@ void Foam::surfaceInterpolation::makeWeights() const
         w[facei] = SfdNei/(SfdOwn + SfdNei);
     }
 
-    surfaceScalarField::Boundary& wBf =
-        weights.boundaryFieldRef();
+    surfaceScalarField::Boundary& wBf = weights.boundaryFieldRef();
 
     forAll(mesh_.boundary(), patchi)
     {
@@ -228,6 +226,7 @@ void Foam::surfaceInterpolation::makeDeltaCoeffs() const
         dimless/dimLength
     );
     surfaceScalarField& deltaCoeffs = *deltaCoeffs_;
+    deltaCoeffs.oriented().oriented() = true;
 
 
     // Set local references to mesh data
@@ -278,6 +277,7 @@ void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
         dimless/dimLength
     );
     surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
+    nonOrthDeltaCoeffs.oriented().oriented() = true;
 
 
     // Set local references to mesh data
@@ -342,6 +342,7 @@ void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
         dimless
     );
     surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
+    corrVecs.oriented().oriented() = true;
 
     // Set local references to mesh data
     const volVectorField& C = mesh_.C();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index 8d7074174fb..ab8f3fcc8c1 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -311,6 +311,8 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
 
     tlambdas.clear();
 
+//    tsf.ref().oriented() = Sf.oriented();
+
     return tsf;
 }
 
@@ -363,6 +365,8 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
         >
     > tsf = dotInterpolate(Sf, vf, weights(vf));
 
+    tsf.ref().oriented() = Sf.oriented();
+
     if (corrected())
     {
         tsf.ref() += Sf & correction(vf);
@@ -397,6 +401,7 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
             surfaceMesh
         >
     > tSfDotinterpVf = dotInterpolate(Sf, tvf());
+
     tvf.clear();
     return tSfDotinterpVf;
 }
-- 
GitLab