Automatic derivative, initial version

applications/utilities/miscellaneous/autoDerivative/testFad/Make/files

@@ -0,0 +1,3 @@
+testFad.C
+
+EXE = $(FOAM_APPBIN)/testFad

applications/utilities/miscellaneous/autoDerivative/testFad/Make/options

@@ -0,0 +1,4 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS =

applications/utilities/miscellaneous/autoDerivative/testFad/testFad.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    testFad
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "FadOne.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    typedef FadOne<1> fadScalar;
+
+    fadScalar x(1);
+    x.deriv(0) = 1;
+
+    fadScalar f1 = 2.0*Foam::pow(x, 2) + 3.0*x + 4.0;
+    Info << "fadScalar of f1 is: "
+         << f1 << " It should be: [9 7]" << endl;
+
+    fadScalar y(1);
+    y.deriv(0) = 1;
+
+    fadScalar f2 = 3*Foam::pow(y, 2)*Foam::sin(Foam::pow(y, 2));
+    Info << "fadScalar of f2 is: " << f2
+         << " It should be: [2.5244129 8.2906397]" << endl;
+
+    fadScalar z(4);
+    z.deriv(0) = 1;
+
+    fadScalar f3 = Foam::sqrt(Foam::pow(z, 3));
+    Info << "fadScalar of f3 is: " << f3
+         << " It should be: [8 3]" << endl;
+
+    fadScalar f4 = Foam::atan(x);
+    Info << "fadScalar of f4 is: " << f4
+         << " It should be: [?, ?]" << endl;
+
+    fadScalar f5 = Foam::log(x);
+    Info << "fadScalar of f4 is: " << f5
+         << " It should be: [?, ?]" << endl;
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/miscellaneous/autoDerivative/testFadField/Make/files

@@ -0,0 +1,4 @@
+testFadField.C
+fadOneFields.C
+
+EXE = $(FOAM_APPBIN)/testFadField

applications/utilities/miscellaneous/autoDerivative/testFadField/Make/options

@@ -0,0 +1,4 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS =

applications/utilities/miscellaneous/autoDerivative/testFadField/fadOneFields.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Specialisation of Field<T> for given FadOne<nVars>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fadOneFields.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.C"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+BINARY_TYPE_OPERATOR(fadScalar, fadScalar, fadScalar, +, add)
+BINARY_TYPE_OPERATOR(fadScalar, fadScalar, fadScalar, -, subtract)
+
+BINARY_OPERATOR(fadScalar, fadScalar, fadScalar, *, multiply)
+BINARY_OPERATOR(fadScalar, fadScalar, fadScalar, /, divide)
+
+BINARY_TYPE_OPERATOR_SF(fadScalar, fadScalar, fadScalar, /, divide)
+
+BINARY_FUNCTION(fadScalar, fadScalar, fadScalar, pow)
+BINARY_TYPE_FUNCTION(fadScalar, fadScalar, fadScalar, pow)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+UNARY_FUNCTION(fadScalar, fadScalar, pos)
+UNARY_FUNCTION(fadScalar, fadScalar, neg)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<fadScalarField>
+fadScalarField::component(const direction) const
+{
+    return *this;
+}
+
+template<>
+void component
+(
+    fadScalarField& lf,
+    const UList<fadScalar>& f,
+    const direction
+)
+{
+    lf = f;
+}
+
+template<>
+void fadScalarField::replace
+(
+    const direction,
+    const UList<fadScalar>& lf
+)
+{
+    *this = lf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// ************************************************************************* //

applications/utilities/miscellaneous/autoDerivative/testFadField/fadOneFields.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    labelField
+
+Description
+    Specialisation of Field<T> for given FadOne<nVars>.
+
+SourceFiles
+    fadOneFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fadOneFields_H
+#define fadOneFields_H
+
+#include "FadOneField.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef FadOne<2> fadScalar;
+typedef Field<fadScalar> fadScalarField;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+BINARY_TYPE_OPERATOR(fadScalar, fadScalar, fadScalar, +, add)
+BINARY_TYPE_OPERATOR(fadScalar, fadScalar, fadScalar, -, subtract)
+
+BINARY_OPERATOR(fadScalar, fadScalar, fadScalar, *, multiply)
+BINARY_OPERATOR(fadScalar, fadScalar, fadScalar, /, divide)
+
+BINARY_TYPE_OPERATOR_SF(fadScalar, fadScalar, fadScalar, /, divide)
+
+BINARY_FUNCTION(fadScalar, fadScalar, fadScalar, pow)
+BINARY_TYPE_FUNCTION(fadScalar, fadScalar, fadScalar, pow)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+UNARY_FUNCTION(fadScalar, fadScalar, pos)
+UNARY_FUNCTION(fadScalar, fadScalar, neg)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<fadScalarField>
+fadScalarField::component
+(
+    const direction
+) const;
+
+
+template<>
+void component
+(
+    fadScalarField& lf,
+    const UList<fadScalar>& f,
+    const direction
+);
+
+template<>
+void fadScalarField::replace
+(
+    const direction,
+    const UList<fadScalar>& lf
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    testFadField
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fadOneFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    fadScalarField tf(3, fadScalar(3.3));
+
+    Info << "name: " << fadScalar::typeName << endl;
+
+    Info << "tf value: " << FadOneValue(tf) << endl;
+    Info << "tf deriv: " << FadOneDeriv(tf, 0) << endl;
+
+    FadOneSetDeriv(tf, 0, Field<double>(3, 1));
+
+    fadScalar ten(10.0);
+    fadScalar bs = tf[0] - ten;
+    fadScalar cs = tf[0] - 10;
+
+    fadScalarField a = sqr(tf);
+    fadScalarField b = tf - ten;
+    fadScalarField c = ten - tf;
+    fadScalarField d = sqr(tf - ten);
+    fadScalarField e = (sqr(tf - bs + 10) + 22.3*sqr(cs - bs + 10));
+    fadScalarField f = d/e;
+//     fadScalarField g = Foam::atan(f);    // HJ, not prepared
+//     fadScalarField h = Foam::log(f);     // HJ, not prepared
+
+
+    Info << "new tf: " << tf << endl;
+
+    Field<Vector<fadScalar> > tfVector(3, Vector<fadScalar>(1, 2, 3));
+
+    Info << "tfVector: " << tfVector << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOne.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOne
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "FadOne.H"
+
+#include <sstream>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<int nVars>
+const char* const FadOne<nVars>::typeName = "FadOne";
+
+template<int nVars>
+const FadOne<nVars> FadOne<nVars>::zero = FadOne<nVars>(0);
+
+template<int nVars>
+const FadOne<nVars> FadOne<nVars>::one = FadOne<nVars>(1);
+
+
+template<int nVars>
+word name(const FadOne<nVars>& u)
+{
+    std::ostringstream osBuffer;
+    osBuffer << u;
+    return osBuffer.str();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOne.H

@@ -0,0 +1,345 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOne
+
+Description
+    Forward automatic derivative of the first order.  Supports multi-variable
+    derivative, where the number of independent variables is specified at
+    compile-time.  Replaces scalar in calculations
+
+SourceFiles
+    FadOneI.H
+    FadOne.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FadOne_H
+#define FadOne_H
+
+#include "error.H"
+#include "IOstream.H"
+#include "token.H"
+
+#include "products.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * Forward declaration of template friend fuctions * * * * * * * //
+
+template<int nVars> class FadOne;
+
+template<int nVars>
+inline FadOne<nVars> operator+(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator+(const double, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator+(const FadOne<nVars>&, const double);
+
+template<int nVars>
+inline FadOne<nVars> operator-(const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator-(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator-(const double, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator-(const FadOne<nVars>&, const double);
+
+template<int nVars>
+inline FadOne<nVars> operator*(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator*(const double, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator*(const FadOne<nVars>&, const double);
+
+template<int nVars>
+inline FadOne<nVars> operator/(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator/(const double, const FadOne<nVars>&);
+
+template<int nVars>
+inline FadOne<nVars> operator/(const FadOne<nVars>&, const double);
+
+template<int nVars>
+inline bool operator<(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline bool operator<=(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline bool operator>(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline bool operator>=(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline bool operator==(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+inline bool operator!=(const FadOne<nVars>&, const FadOne<nVars>&);
+
+template<int nVars>
+Istream& operator>>(Istream&, FadOne<nVars>&);
+
+template<int nVars>
+Ostream& operator<<(Ostream&, const FadOne<nVars>&);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class FadOne Declaration
+\*---------------------------------------------------------------------------*/
+
+template<int nVars>
+class FadOne
+{
+    // Private data
+
+        //- Length of data list
+        static const label Length = nVars + 1;
+
+        //- Value and derivatives
+        double data_[Length];
+
+
+    // Private Member Functions
+
+
+public:
+
+    //- Component type
+    typedef FadOne<nVars> cmptType;
+
+
+    // Member constants
+
+        enum
+        {
+            dim = 3,         // Dimensionality of space
+            rank = 0,        // Rank od Scalar is 0
+            nComponents = 1  // Number of components in FadOne is 1
+        };
+
+    // Static data members
+
+        static const char* const typeName;
+        static const char* componentNames[];
+        static const FadOne<nVars> zero;
+        static const FadOne<nVars> one;
+
+
+    // Constructors
+
+        //- Construct null
+        inline FadOne();
+
+        //- Construct from value
+        inline FadOne(const double& v);
+
+        //- Construct from Istream
+        inline FadOne(Istream&);
+
+        //- Construct as copy
+        inline FadOne(const FadOne<nVars>&);
+
+
+    // Destructor - default
+
+
+    // Member Functions
+
+        //- Return value
+        inline const double& value() const;
+        inline double& value();
+
+        //- Return derivative
+        inline const double& deriv(const label i) const;
+        inline double& deriv(const label i);
+
+
+    // Member Operators
+
+        inline void operator=(const double&);
+        inline void operator=(const FadOne<nVars>&);
+
+        inline void operator+=(const double&);
+        inline void operator+=(const FadOne<nVars>&);
+        inline void operator-=(const double&);
+        inline void operator-=(const FadOne<nVars>&);
+        inline void operator*=(const double&);
+        inline void operator*=(const FadOne<nVars>&);
+        inline void operator/=(const double&);
+        inline void operator/=(const FadOne<nVars>&);
+
+
+    // Friend Functions
+
+
+    // Friend operators
+
+        friend FadOne<nVars> operator+ <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator+ <nVars>
+        (
+            const double,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator+ <nVars>
+        (
+            const FadOne<nVars>&,
+            const double
+        );
+
+        friend FadOne<nVars> operator- <nVars>(const FadOne<nVars>&);
+
+        friend FadOne<nVars> operator- <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator- <nVars>
+        (
+            const double,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator- <nVars>
+        (
+            const FadOne<nVars>&,
+            const double
+        );
+
+        friend FadOne<nVars> operator* <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator* <nVars>
+        (
+            const double,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator* <nVars>
+        (
+            const FadOne<nVars>&,
+            const double
+        );
+
+        friend FadOne<nVars> operator/ <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator/ <nVars>
+        (
+            const double,
+            const FadOne<nVars>&
+        );
+
+        friend FadOne<nVars> operator/ <nVars>
+        (
+            const FadOne<nVars>&,
+            const double
+        );
+
+        friend bool operator< <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend bool operator<= <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend bool operator> <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend bool operator>= <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend bool operator== <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+        friend bool operator!= <nVars>
+        (
+            const FadOne<nVars>&,
+            const FadOne<nVars>&
+        );
+
+
+    // IOstream Operators
+
+        friend Istream& operator>> <nVars>(Istream&, FadOne<nVars>&);
+        friend Ostream& operator<< <nVars>(Ostream&, const FadOne<nVars>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "FadOneI.H"
+#include "FadOne.C"
+#include "FadOneFunctions.H"
+#include "FadOneFunctionsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOneField.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOneFields
+
+\*----------------------------------------------------------------------------*/
+
+#include "FadOneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<int nVars>
+Foam::tmp<Foam::Field<Foam::FadOne<nVars> > >
+Foam::ValueFadOneField(const UList<scalar>& u)
+{
+    tmp<Foam::Field<FadOne<nVars> > > tr
+    (
+        new Foam::Field<FadOne<nVars> > (u.size())
+    );
+
+    Field<FadOne<nVars> >& r = tr();
+
+    forAll (r, i)
+    {
+        r[i].value() = u[i];
+    }
+
+    return tr;
+}
+
+
+template<int nVars>
+Foam::tmp<Foam::scalarField> Foam::FadOneValue(const Field<FadOne<nVars> >& u)
+{
+    tmp<scalarField> tr(new scalarField(u.size()));
+    scalarField& r = tr();
+
+    forAll (r, i)
+    {
+        r[i] = u[i].value();
+    }
+
+    return tr;
+}
+
+
+template<int nVars>
+void Foam::FadOneSetValue
+(
+    Field<FadOne<nVars> >& u,
+    const scalarField& val
+)
+{
+    forAll (u, i)
+    {
+        u[i].value() = val[i];
+    }
+}
+
+
+template<int nVars>
+Foam::tmp<Foam::scalarField> Foam::FadOneDeriv
+(
+    const Field<FadOne<nVars> >& u,
+    const direction d
+)
+{
+    tmp<scalarField> tr(new scalarField(u.size()));
+    scalarField& r = tr();
+
+    forAll (r, i)
+    {
+        r[i] = u[i].deriv(d);
+    }
+
+    return tr;
+}
+
+
+template<int nVars>
+void Foam::FadOneSetDeriv
+(
+    Field<FadOne<nVars> >& u,
+    const direction d,
+    const scalarField& der
+)
+{
+    forAll (u, i)
+    {
+        u[i].deriv(d) = der[i];
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOneField.H

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOneField
+
+Description
+    Specialisation of Field<T> for FadOne
+
+SourceFiles
+    FadOneFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FadOneField_H
+#define FadOneField_H
+
+#include "FadOne.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<int nVars>
+tmp<Field<FadOne<nVars> > > ValueFadOneField(const UList<scalar>&);
+
+
+template<int nVars>
+tmp<scalarField> FadOneValue(const Field<FadOne<nVars> >&);
+
+template<int nVars>
+void FadOneSetValue
+(
+    Field<FadOne<nVars> >&,
+    const scalarField&
+);
+
+template<int nVars>
+tmp<scalarField> FadOneDeriv
+(
+    const Field<FadOne<nVars> >&,
+    const direction
+);
+
+template<int nVars>
+void FadOneSetDeriv
+(
+    Field<FadOne<nVars> >&,
+    const direction,
+    const scalarField&
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+#ifdef NoRepository
+#   include "FadOneField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOneFunctions.H

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Generic field type.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FadOne.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * */
+
+// Return a string representation of a complex
+template<int nVars>
+word name(const FadOne<nVars>&);
+
+
+// Allow binary writing for arrays of complex
+template<int nVars>
+inline bool writeBinary(const FadOne<nVars>*)
+{
+    return true;
+}
+
+
+template<int nVars>
+inline const FadOne<nVars>& max
+(
+    const FadOne<nVars>&,
+    const FadOne<nVars>&
+);
+
+template<int nVars>
+inline const FadOne<nVars>& min
+(
+    const FadOne<nVars>&,
+    const FadOne<nVars>&
+);
+
+template<int nVars>
+inline FadOne<nVars> limit
+(
+    const FadOne<nVars>&,
+    const FadOne<nVars>&
+);
+
+
+template<int nVars>
+inline FadOne<nVars>& setComponent(FadOne<nVars>& u, const direction)
+{
+    return u;
+}
+
+template<int nVars>
+inline FadOne<nVars> component(const FadOne<nVars>& u, const direction)
+{
+    return u;
+}
+
+template<int nVars>
+inline FadOne<nVars> mag(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> sin(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> asin(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> cos(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> acos(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> tan(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> atan(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> tanh(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> exp(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> log(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> sqr(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> sqrt(const FadOne<nVars>& u);
+
+template<int nVars>
+inline FadOne<nVars> pow
+(
+    const FadOne<nVars>&,
+    const FadOne<nVars>&
+);
+
+template<int nVars>
+inline FadOne<nVars> pow
+(
+    const FadOne<nVars>&,
+    const double
+);
+
+template<int nVars>
+inline FadOne<nVars> pow
+(
+    const double&,
+    const FadOne<nVars>&
+);
+
+template<int nVars>
+inline label sign(const FadOne<nVars>& u);
+
+template<int nVars>
+inline label pos(const FadOne<nVars>& u);
+
+template<int nVars>
+inline label neg(const FadOne<nVars>& u);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/primitives/FadOne/FadOneFunctionsI.H

@@ -0,0 +1,311 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOne
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::mag(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(u.value());
+
+    for (int i = 0; i < nVars; i++)
+    {
+        if (u.deriv(i) == 0)
+        {
+            r.deriv(i) = 0;
+        }
+        else
+        {
+            if (u.value() > 0)
+            {
+                r.deriv(i) = u.deriv(i);
+            }
+            else
+            {
+                r.deriv(i) = -u.deriv(i);
+            }
+        }
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::sin(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::sin(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = u.deriv(i)*Foam::cos(u.value());
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::asin(const FadOne<nVars>& u)
+{
+    const double& x = u.value();
+
+    FadOne<nVars> r(Foam::asin(x));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        const double& d = u.deriv(i);
+        r.deriv(i) = d/Foam::sqrt(1 - Foam::sqr(x));
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::cos(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::cos(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = -u.deriv(i)*Foam::sin(u.value());
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::acos(const FadOne<nVars>& u)
+{
+    const double& x = u.value();
+
+    FadOne<nVars> r(Foam::acos(x));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        const double& d = u.deriv(i);
+        r.deriv(i) = -d/Foam::sqrt(1 - Foam::sqr(x));
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::tan(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::tan(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = u.deriv(i)/(Foam::sqr(Foam::cos(u.value())));
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::atan(const FadOne<nVars>& u)
+{
+    const double& x = u.value();
+
+    FadOne<nVars> r(Foam::atan(x));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        const double& d = u.deriv(i);
+        r.deriv(i) = d/Foam::sqrt(1 + Foam::sqr(x));
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::tanh(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::tanh(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = u.deriv(i)*(1 - Foam::sqr(Foam::tanh(u.value())));
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::exp(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::exp(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = u.deriv(i)*Foam::exp(u.value());
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::log(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::log(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = u.deriv(i)/u.value();
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::sqr(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::sqr(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = 2*u.deriv(i)*u.value();
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::sqrt(const FadOne<nVars>& u)
+{
+    FadOne<nVars> r(Foam::sqrt(u.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        if (u.deriv(i) > 0)
+        {
+            r.deriv(i) = 0.5*u.deriv(i)/Foam::sqrt(u.value());
+        }
+        else
+        {
+            r.deriv(i) = 0;
+        }
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::pow
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(Foam::pow(u.value(), v.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) =
+            v.deriv(i)*Foam::log(u.value())*Foam::pow(u.value(), v.value())
+          + v.value()*u.deriv(i)*Foam::pow(u.value(), v.value() - 1);
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::pow
+(
+    const FadOne<nVars>& u,
+    const double v
+)
+{
+    FadOne<nVars> r(Foam::pow(u.value(), v));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = v*u.deriv(i)*Foam::pow(u.value(), v - 1);
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::pow
+(
+    const double& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(Foam::pow(u, v.value()));
+
+    for (int i = 0; i < nVars; i++)
+    {
+        r.deriv(i) = v.deriv(i)*Foam::log(u)*Foam::pow(u, v.value());
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::label Foam::sign(const FadOne<nVars>& u)
+{
+    return (u.value() >= 0) ? 1: -1;
+}
+
+
+template<int nVars>
+inline Foam::label Foam::pos(const FadOne<nVars>& u)
+{
+    return (u.value() >= 0) ? 1: 0;
+}
+
+
+template<int nVars>
+inline Foam::label Foam::neg(const FadOne<nVars>& u)
+{
+    return (u.value() < 0) ? 1: 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/FadOne/FadOneI.H

@@ -0,0 +1,531 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    FadOne
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<int nVars>
+inline Foam::FadOne<nVars>::FadOne()
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] = 0;
+    }
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars>::FadOne(const double& v)
+{
+    this->value() = v;
+
+    for (int i = 0; i < nVars; i++)
+    {
+        this->deriv(i) = 0;
+    }
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars>::FadOne(Istream& is)
+{
+    is >> *this;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars>::FadOne(const FadOne<nVars>& f)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] = f.data_[i];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<int nVars>
+inline const double& Foam::FadOne<nVars>::value() const
+{
+    return data_[0];
+}
+
+
+template<int nVars>
+inline double& Foam::FadOne<nVars>::value()
+{
+    return data_[0];
+}
+
+
+template<int nVars>
+inline const double& Foam::FadOne<nVars>::deriv(const label i) const
+{
+#   ifdef FULLDEBUG
+    if (i >= nVars)
+    {
+        FatalErrorIn
+        (
+            "double& FadOne<nVars>::deriv(const label i)"
+        )   << "index out of range"
+            << abort(FatalError);
+    }
+#   endif
+
+    return data_[i + 1];
+}
+
+
+template<int nVars>
+inline double& Foam::FadOne<nVars>::deriv(const label i)
+{
+#   ifdef FULLDEBUG
+    if (i >= nVars)
+    {
+        FatalErrorIn
+        (
+            "double& FadOne<nVars>::deriv(const label i)"
+        )   << "index out of range"
+            << abort(FatalError);
+    }
+#   endif
+
+    return data_[i + 1];
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator=(const double& d)
+{
+    for (int i = 0; i < nVars; i++)
+    {
+        this->deriv(i) = 0;
+    }
+
+    this->value() = d;
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator=(const FadOne<nVars>& v)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] = v.data_[i];
+    }
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator+=(const double& d)
+{
+    this->value() += d;
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator+=(const FadOne<nVars>& v)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] += v.data_[i];
+    }
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator-=(const double& d)
+{
+    this->value() -= d;
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator-=(const FadOne<nVars>& v)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] -= v.data_[i];
+    }
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator*=(const double& d)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] *= d;
+    }
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator*=(const FadOne<nVars>& v)
+{
+    for (int i = 0; i < nVars; i++)
+    {
+        this->deriv(i) = this->deriv(i)*v.value() + this->value()*v.deriv(i);
+    }
+
+    this->value() *= v.value();
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator/=(const double& d)
+{
+    for (int i = 0; i < Length; i++)
+    {
+        data_[i] /= d;
+    }
+}
+
+
+template<int nVars>
+inline void Foam::FadOne<nVars>::operator/=(const FadOne<nVars>& v)
+{
+    double& lv = this->value();
+    const double& vv = v.value();
+
+    for (int i = 0; i < nVars; i++)
+    {
+        this->deriv(i) = (this->deriv(i)*vv - lv*v.deriv(i))/(vv*vv);
+    }
+
+    lv /= vv;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator+
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(u);
+    r += v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator+
+(
+    const FadOne<nVars>& u,
+    const double v
+)
+{
+    FadOne<nVars> r(u);
+    r += v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator+
+(
+    const double u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(v);
+    r += u;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator-
+(
+    const FadOne<nVars>& u
+)
+{
+    FadOne<nVars> r;
+    for (int i = 0; i < FadOne<nVars>::Length; i++)
+    {
+        r.data_[i] = -u.data_[i];
+    }
+
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator-
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(u);
+    r -= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator-
+(
+    const FadOne<nVars>& u,
+    const double v
+)
+{
+    FadOne<nVars> r(u);
+    r -= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator-
+(
+    const double u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(v);
+    r -= u;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator*
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(u);
+    r *= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator*
+(
+    const FadOne<nVars>& u,
+    const double v
+)
+{
+    FadOne<nVars> r(u);
+    r *= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator*
+(
+    const double u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(v);
+    r *= u;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator/
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(u);
+    r /= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator/
+(
+    const FadOne<nVars>& u,
+    const double v
+)
+{
+    FadOne<nVars> r(u);
+    r /= v;
+    return r;
+}
+
+
+template<int nVars>
+inline Foam::FadOne<nVars> Foam::operator/
+(
+    const double u,
+    const FadOne<nVars>& v
+)
+{
+    FadOne<nVars> r(v);
+    r /= u;
+    return r;
+}
+
+
+template<int nVars>
+inline bool Foam::operator<
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() < v.value();
+}
+
+
+template<int nVars>
+inline bool Foam::operator<=
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() <= v.value();
+}
+
+
+template<int nVars>
+inline bool Foam::operator>
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() > v.value();
+}
+
+
+template<int nVars>
+inline bool Foam::operator>=
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() >= v.value();
+}
+
+
+template<int nVars>
+inline bool Foam::operator==
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() == v.value();
+}
+
+
+template<int nVars>
+inline bool Foam::operator!=
+(
+    const FadOne<nVars>& u,
+    const FadOne<nVars>& v
+)
+{
+    return u.value() != v.value();
+}
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+
+template<int nVars>
+Foam::Istream& Foam::operator>>(Istream& is, FadOne<nVars>& f)
+{
+    // Read begining of FadOne
+    if (token(is) != token::BEGIN_SQR)
+    {
+        Info<< "expected a " << token::BEGIN_SQR << " in FadOne<nVars>"
+            << endl << "in stream " << is.info() << endl;
+    }
+
+    for (int i = 0; i < FadOne<nVars>::Length ; i++)
+    {
+        is >> f.data_[i];
+    }
+
+    // Check end of FadOne
+    if (token(is) != token::END_SQR)
+    {
+        Info<< "expected a " << token::END_SQR << " in FadOne<nVars>"
+            << endl << "in stream " << is.info() << endl;
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, FadOne<nVars>&)");
+
+    return is;
+}
+
+
+template<int nVars>
+Foam::Ostream& Foam::operator<<(Ostream& os, const FadOne<nVars>& f)
+{
+    os << token::BEGIN_SQR;
+
+    for (int i = 0; i < FadOne<nVars>::Length - 1; i++)
+    {
+        os << f.data_[i] << token::SPACE;
+    }
+    os << f.data_[FadOne<nVars>::Length -1 ] << token::END_SQR;
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const FadOne<nVars>&)");
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //