ENH: Add triad primitive type and Field/IOField

src/OpenFOAM/Make/files

@@ -63,6 +63,7 @@ primitives/globalIndexAndTransform/globalIndexAndTransform.C
 primitives/globalIndexAndTransform/vectorTensorTransform/vectorTensorTransform.C
 primitives/quaternion/quaternion.C
 primitives/septernion/septernion.C
+primitives/triad/triad.C
 
 /* functions, data entries */
 primitives/functions/DataEntry/makeDataEntries.C
@@ -524,6 +525,7 @@ $(Fields)/sphericalTensorField/sphericalTensorField.C
 $(Fields)/diagTensorField/diagTensorField.C
 $(Fields)/symmTensorField/symmTensorField.C
 $(Fields)/tensorField/tensorField.C
+$(Fields)/triadField/triadField.C
 $(Fields)/complexFields/complexFields.C
 
 $(Fields)/labelField/labelIOField.C
@@ -542,6 +544,7 @@ $(Fields)/symmTensorField/symmTensorIOField.C
 $(Fields)/symmTensorField/symmTensorFieldIOField.C
 $(Fields)/tensorField/tensorIOField.C
 $(Fields)/tensorField/tensorFieldIOField.C
+$(Fields)/triadField/triadIOField.C
 $(Fields)/transformField/transformField.C
 
 pointPatchFields = fields/pointPatchFields

src/OpenFOAM/fields/Fields/fieldTypes.H

@@ -36,6 +36,7 @@ Description
 #include "sphericalTensor.H"
 #include "symmTensor.H"
 #include "tensor.H"
+#include "triad.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/fields/Fields/triadField/triadField.C

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "triadField.H"
+#include "transformField.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.C"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/triadField/triadField.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Typedef
+    Foam::triadField
+
+Description
+    Specialisation of Field\<T\> for triad.
+
+SourceFiles
+    triadField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triadField_H
+#define triadField_H
+
+#include "Field.H"
+#include "triad.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef Field<triad> triadField;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/triadField/triadIOField.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Description
+    triadField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triadIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName(triadIOField, "triadField", 0);
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/triadField/triadIOField.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Typedef
+    Foam::triadIOField
+
+Description
+    triadField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triadIOField_H
+#define triadIOField_H
+
+#include "triadField.H"
+#include "IOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<triad> triadIOField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/triad/triad.C

@@ -0,0 +1,338 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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 "triad.H"
+#include "transform.H"
+#include "tensor.H"
+#include "quaternion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+const char* const triad::Vector<vector>::typeName = "triad";
+
+template<>
+const char* triad::Vector<vector>::componentNames[] = {"x", "y", "z"};
+
+const triad triad::zero(vector::zero, vector::zero, vector::zero);
+
+const triad triad::one(vector::one, vector::one, vector::one);
+
+const triad triad::max(vector::max, vector::max, vector::max);
+
+const triad triad::min(vector::min, vector::min, vector::min);
+
+const triad triad::unset(triad::max);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::triad::triad(const quaternion& q)
+{
+    tensor Rt(q.R().T());
+    x() = Rt.x();
+    y() = Rt.y();
+    z() = Rt.z();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::triad::orthogonalize()
+{
+    // Hack for 2D z-slab cases
+    // if (!set(2))
+    // {
+    //     operator[](2) = vector(0, 0, 1);
+    // }
+
+    // If only two of the axes are set, set the third
+    if (set(0) && set(1) && !set(2))
+    {
+        operator[](2) = orthogonal(operator[](0), operator[](1));
+    }
+    else if (set(0) && set(2) && !set(1))
+    {
+        operator[](1) = orthogonal(operator[](0), operator[](2));
+    }
+    else if (set(1) && set(2) && !set(0))
+    {
+        operator[](0) = orthogonal(operator[](1), operator[](2));
+    }
+
+    // If all the axes are set
+    if (set())
+    {
+        for (int i=0; i<2; i++)
+        {
+        scalar o01 = mag(operator[](0) & operator[](1));
+        scalar o02 = mag(operator[](0) & operator[](2));
+        scalar o12 = mag(operator[](1) & operator[](2));
+
+        if (o01 < o02 && o01 < o12)
+        {
+            operator[](2) = orthogonal(operator[](0), operator[](1));
+
+            // if (o02 < o12)
+            // {
+            //     operator[](1) = orthogonal(operator[](0), operator[](2));
+            // }
+            // else
+            // {
+            //     operator[](0) = orthogonal(operator[](1), operator[](2));
+            // }
+        }
+        else if (o02 < o12)
+        {
+            operator[](1) = orthogonal(operator[](0), operator[](2));
+
+            // if (o01 < o12)
+            // {
+            //     operator[](2) = orthogonal(operator[](0), operator[](1));
+            // }
+            // else
+            // {
+            //     operator[](0) = orthogonal(operator[](1), operator[](2));
+            // }
+        }
+        else
+        {
+            operator[](0) = orthogonal(operator[](1), operator[](2));
+
+            // if (o02 < o01)
+            // {
+            //     operator[](1) = orthogonal(operator[](0), operator[](2));
+            // }
+            // else
+            // {
+            //     operator[](2) = orthogonal(operator[](0), operator[](1));
+            // }
+        }
+        }
+    }
+}
+
+
+void Foam::triad::operator+=(const triad& t2)
+{
+    if (t2.set(0) && !set(0))
+    {
+        operator[](0) = t2.operator[](0);
+    }
+
+    if (t2.set(1) && !set(1))
+    {
+        operator[](1) = t2.operator[](1);
+    }
+
+    if (t2.set(2) && !set(2))
+    {
+        operator[](2) = t2.operator[](2);
+    }
+
+    if (set() && t2.set())
+    {
+        direction correspondance[3];
+        short signd[3];
+
+        for (direction i=0; i<3; i++)
+        {
+            scalar mostAligned = -1;
+            for (direction j=0; j<3; j++)
+            {
+                bool set = false;
+                for (direction k=0; k<i; k++)
+                {
+                    if (correspondance[k] == j)
+                    {
+                        set = true;
+                        break;
+                    }
+                }
+
+                if (!set)
+                {
+                    scalar a = operator[](i) & t2.operator[](j);
+                    scalar maga = mag(a);
+
+                    if (maga > mostAligned)
+                    {
+                        correspondance[i] = j;
+                        mostAligned = maga;
+                        signd[i] = sign(a);
+                    }
+                }
+            }
+        }
+
+        for (direction i=0; i<3; i++)
+        {
+            operator[](i) += signd[i]*t2.operator[](correspondance[i]);
+        }
+    }
+}
+
+
+void Foam::triad::align(const vector& v)
+{
+    if (set())
+    {
+        vector mostAligned
+        (
+            mag(v & operator[](0)),
+            mag(v & operator[](1)),
+            mag(v & operator[](2))
+        );
+
+        scalar mav;
+
+        if
+        (
+            mostAligned.x() > mostAligned.y()
+         && mostAligned.x() > mostAligned.z()
+        )
+        {
+            mav = mostAligned.x();
+            mostAligned = operator[](0);
+        }
+        else if (mostAligned.y() > mostAligned.z())
+        {
+            mav = mostAligned.y();
+            mostAligned = operator[](1);
+        }
+        else
+        {
+            mav = mostAligned.z();
+            mostAligned = operator[](2);
+        }
+
+        if (mav < 0.99)
+        {
+            tensor R(rotationTensor(mostAligned, v));
+
+            operator[](0) = transform(R, operator[](0));
+            operator[](1) = transform(R, operator[](1));
+            operator[](2) = transform(R, operator[](2));
+        }
+    }
+}
+
+
+Foam::triad Foam::triad::sortxyz() const
+{
+    triad t;
+
+    if
+    (
+        mag(operator[](0).x()) > mag(operator[](1).x())
+     && mag(operator[](0).x()) > mag(operator[](2).x())
+    )
+    {
+        t[0] = operator[](0);
+
+        if (mag(operator[](1).y()) > mag(operator[](2).y()))
+        {
+            t[1] = operator[](1);
+            t[2] = operator[](2);
+        }
+        else
+        {
+            t[1] = operator[](2);
+            t[2] = operator[](1);
+        }
+    }
+    else if
+    (
+        mag(operator[](1).x()) > mag(operator[](2).x())
+    )
+    {
+        t[0] = operator[](1);
+
+        if (mag(operator[](0).y()) > mag(operator[](2).y()))
+        {
+            t[1] = operator[](0);
+            t[2] = operator[](2);
+        }
+        else
+        {
+            t[1] = operator[](2);
+            t[2] = operator[](0);
+        }
+    }
+    else
+    {
+        t[0] = operator[](2);
+
+        if (mag(operator[](0).y()) > mag(operator[](1).y()))
+        {
+            t[1] = operator[](0);
+            t[2] = operator[](1);
+        }
+        else
+        {
+            t[1] = operator[](1);
+            t[2] = operator[](0);
+        }
+    }
+
+    if (t[0].x() < 0) t[0] *= -1;
+    if (t[1].y() < 0) t[1] *= -1;
+    if (t[2].z() < 0) t[2] *= -1;
+
+    return t;
+}
+
+
+
+Foam::triad::operator quaternion() const
+{
+    tensor R;
+
+    R.xx() = x().x();
+    R.xy() = y().x();
+    R.xz() = z().x();
+
+    R.yx() = x().y();
+    R.yy() = y().y();
+    R.yz() = z().y();
+
+    R.zx() = x().z();
+    R.zy() = y().z();
+    R.zz() = z().z();
+
+    return quaternion(R);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/primitives/triad/triad.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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::triad
+
+Description
+    Representation of a 3D Cartesian coordinate system as a Vector of vectors.
+
+See Also
+    Foam::quaternion
+
+SourceFiles
+    triadI.H
+    triad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triad_H
+#define triad_H
+
+#include "vector.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class triad;
+Istream& operator>>(Istream&, triad&);
+Ostream& operator<<(Ostream&, const triad&);
+
+class quaternion;
+
+
+/*---------------------------------------------------------------------------*\
+                           Class triad Declaration
+\*---------------------------------------------------------------------------*/
+
+class triad
+:
+    public Vector<vector>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline triad();
+
+        //- Construct from components
+        inline triad(const Vector<vector>& vv);
+
+        //- Construct from coordinate axes
+        inline triad(const vector& x, const vector& y, const vector& z);
+
+        //- Construct from a primary axis with the other two unset
+        inline triad(const vector& pa);
+
+        //- Construct from a quaternion
+        triad(const quaternion& q);
+
+        //- Construct from Istream
+        inline triad(Istream&);
+
+
+    // Static data members
+
+        static const triad zero;
+        static const triad one;
+        static const triad max;
+        static const triad min;
+        static const triad unset;
+
+
+    // Member Functions
+
+        //- Is the vector in the direction d set
+        inline bool set(const direction d) const;
+
+        //- Are all the vector set
+        inline bool set() const;
+
+        //- Return the primary direction of the vector v
+        static inline direction primaryDirection(const vector& v);
+
+        //- Return the vector orthogonal to the two provided
+        static inline vector orthogonal(const vector& v1, const vector& v2);
+
+        //- Orthogonalize this triad so that it is ortho-normal
+        void orthogonalize();
+
+        //- Normalize each set axis vector to have a unit magnitude
+        void normalize();
+
+        //- Align this triad with the given vector v
+        //  by rotating the most aligned axis to be coincident with v
+        void align(const vector& v);
+
+        //- Sort the axes such that they are closest to the x, y and z axes
+        triad sortxyz() const;
+
+        //- convert to a quaternion
+        operator quaternion() const;
+
+
+    // Member Operators
+
+        inline void operator=(const Vector<vector>&);
+
+        //- Add the triad t2 to this triad
+        //  without normalizing or orthogonalizing
+        void operator+=(const triad& t2);
+
+
+    // IOstream Operators
+
+        inline friend Istream& operator>>(Istream&, triad&);
+        inline friend Ostream& operator<<(Ostream&, const triad&);
+};
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+//- Data associated with quaternion type are contiguous
+template<>
+inline bool contiguous<triad>() {return true;}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "triadI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/triad/triadI.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::triad::triad()
+:
+    Vector<vector>(triad::unset)
+{}
+
+
+inline Foam::triad::triad(const Vector<vector>& vv)
+:
+    Vector<vector>(vv)
+{}
+
+
+inline Foam::triad::triad(const vector& x, const vector& y, const vector& z)
+:
+    Vector<vector>(x, y, z)
+{}
+
+
+inline Foam::triad::triad(const vector& pa)
+{
+    operator=(triad::unset);
+    operator[](primaryDirection(pa)) = pa;
+}
+
+
+inline Foam::triad::triad(Istream& is)
+:
+    Vector<vector>(is)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::triad::set(const direction d) const
+{
+    return operator[](d)[0] < GREAT; // vector::zero;
+}
+
+
+inline bool Foam::triad::set() const
+{
+    return set(0) && set(1) && set(2);
+}
+
+
+inline Foam::direction Foam::triad::primaryDirection(const vector& v)
+{
+    if (mag(v.x()) > mag(v.y()) && mag(v.x()) > mag(v.z()))
+    {
+        return triad::X;
+    }
+    else if (mag(v.y()) > mag(v.z()))
+    {
+        return triad::Y;
+    }
+    else
+    {
+        return triad::Z;
+    }
+}
+
+
+inline Foam::vector Foam::triad::orthogonal
+(
+    const vector& v1,
+    const vector& v2
+)
+{
+    vector v3 = v1 ^ v2;
+
+    scalar magV3 = mag(v3);
+
+    if (magV3 > 0.5)
+    {
+        return v3/magV3;
+    }
+    else
+    {
+        return triad::unset[0];
+    }
+}
+
+
+inline void Foam::triad::normalize()
+{
+    if (set(0)) operator[](0) /= mag(operator[](0));
+    if (set(1)) operator[](1) /= mag(operator[](1));
+    if (set(2)) operator[](2) /= mag(operator[](2));
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::triad::operator=(const Vector<vector>& vv)
+{
+    Vector<vector>::operator=(vv);
+}
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+inline Foam::Istream& Foam::operator>>(Istream& is, triad& t)
+{
+    is >> static_cast<Vector<vector>&>(t);
+    return is;
+}
+
+
+inline Foam::Ostream& Foam::operator<<(Ostream& os, const triad& t)
+{
+    os << static_cast<const Vector<vector>&>(t);
+    return os;
+}
+
+
+// ************************************************************************* //