Compilation, A

src/OpenFOAM/Make/files

@@ -414,6 +414,7 @@ $(primitiveMesh)/primitiveMeshFaceCentresAndAreas.C
 $(primitiveMesh)/primitiveMeshFindCell.C
 $(primitiveMesh)/primitiveMeshPointCells.C
 $(primitiveMesh)/primitiveMeshPointFaces.C
+$(primitiveMesh)/primitiveMeshPointEdges.C
 $(primitiveMesh)/primitiveMeshPointPoints.C
 $(primitiveMesh)/primitiveMeshCellPoints.C
 $(primitiveMesh)/primitiveMeshCalcCellShapes.C

src/OpenFOAM/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C

@@ -28,6 +28,7 @@ License
 
 #include "Switch.H"
 #include "mathematicalConstants.H"
+#include "dictionary.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/OpenFOAM/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C

@@ -28,6 +28,7 @@ License
 
 #include "Switch.H"
 #include "mathematicalConstants.H"
+#include "dictionary.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/OpenFOAM/coordinateSystems/coordinateRotation/axisCoordinateRotation.C

@@ -28,6 +28,7 @@ License
 
 #include "Switch.H"
 #include "mathematicalConstants.H"
+#include "dictionary.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/OpenFOAM/coordinateSystems/cylindricalCS.C

@@ -26,7 +26,6 @@ License
 
 #include "cylindricalCS.H"
 
-#include "one.H"
 #include "Switch.H"
 #include "mathematicalConstants.H"
 #include "addToRunTimeSelectionTable.H"

src/OpenFOAM/coordinateSystems/ellipticCylindricalCS.C

@@ -25,8 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ellipticCylindricalCS.H"
-#include "addToRunTimeSelectionTable.H"
+
+#include "Switch.H"
 #include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -59,6 +61,7 @@ Foam::ellipticCylindricalCS::ellipticCylindricalCS
 )
 :
     coordinateSystem(cs),
+    a_(0),
     inDegrees_(inDegrees)
 {}
 
@@ -71,6 +74,22 @@ Foam::ellipticCylindricalCS::ellipticCylindricalCS
 )
 :
     coordinateSystem(name, cs),
+    a_(0),
+    inDegrees_(inDegrees)
+{}
+
+
+Foam::ellipticCylindricalCS::ellipticCylindricalCS
+(
+    const word& name,
+    const point& origin,
+    const coordinateRotation& cr,
+    const scalar a,
+    const bool inDegrees
+)
+:
+    coordinateSystem(name, origin, cr),
+    a_(a),
     inDegrees_(inDegrees)
 {}
 
@@ -91,21 +110,6 @@ Foam::ellipticCylindricalCS::ellipticCylindricalCS
 {}
 
 
-Foam::ellipticCylindricalCS::ellipticCylindricalCS
-(
-    const word& name,
-    const point& origin,
-    const coordinateRotation& cr,
-    const scalar a,
-    const bool inDegrees
-)
-:
-    coordinateSystem(name, origin, cr),
-    a_(a),
-    inDegrees_(inDegrees)
-{}
-
-
 Foam::ellipticCylindricalCS::ellipticCylindricalCS
 (
     const word& name,
@@ -128,8 +132,7 @@ Foam::vector Foam::ellipticCylindricalCS::localToGlobal
 {
     // Notation: u = local.x() v = local.y() z = local.z();
     scalar theta =
-        local.y()
+        local.y()*( inDegrees_ ? mathematicalConstant::pi/180.0 : 1.0 );
-       *( inDegrees_ ? mathematicalConstant::pi/180.0 : 1.0 )
 
     return coordinateSystem::localToGlobal
     (
@@ -150,8 +153,8 @@ Foam::tmp<Foam::vectorField> Foam::ellipticCylindricalCS::localToGlobal
 ) const
 {
     scalarField theta =
-        local.component(vector::Y)
+        local.component(vector::Y)*
-       *( inDegrees_ ? mathematicalConstant::pi/180.0 : 1.0 )
+        ( inDegrees_ ? mathematicalConstant::pi/180.0 : 1.0 );
 
     vectorField lc(local.size());
     lc.replace

src/OpenFOAM/coordinateSystems/ellipticCylindricalCS.H

@@ -108,14 +108,14 @@ public:
         ellipticCylindricalCS(const bool inDegrees = true);
 
         //- Construct copy
-        sphericalCS
+        ellipticCylindricalCS
         (
             const coordinateSystem&,
             const bool inDegrees = true
         );
 
         //- Construct copy with a different name
-        sphericalCS
+        ellipticCylindricalCS
         (
             const word& name,
             const coordinateSystem&,
@@ -132,6 +132,17 @@ public:
             const bool inDegrees = true
         );
 
+        //- Construct from origin and 2 axes
+        ellipticCylindricalCS
+        (
+            const word& name,
+            const point& origin,
+            const vector& axis,
+            const vector& dirn,
+            const scalar a,
+            const bool inDegrees = true
+        );
+
         //- Construct from dictionary
         ellipticCylindricalCS
         (

src/OpenFOAM/db/IOobject/IOobjectI.H

@@ -33,16 +33,20 @@ inline Stream& Foam::IOobject::writeBanner(Stream& os, bool noHint)
 {
     static bool spacesSet = false;
     static char spaces[40];
-    static char spacesRev[40];
 
     if (!spacesSet)
     {
         memset(spaces, ' ', 40);
-        spaces[38 - strlen(Foam::FOAMversion)] = '\0';
-
-        memset(spacesRev, ' ', 40);
-        spacesRev[38 - strlen(Foam::FOAMDevRevisionNumber)] = '\0';
 
+        size_t len = strlen(Foam::FOAMversion);
+        if (len < 38)
+        {
+            spaces[38 - len] = '\0';
+        }
+        else
+        {
+            spaces[0] = '\0';
+        }
         spacesSet = true;
     }
 
@@ -59,10 +63,10 @@ inline Stream& Foam::IOobject::writeBanner(Stream& os, bool noHint)
 
     os  <<
         "| =========                 |                                                 |\n"
-        "| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\n"
+        "| \\\\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |\n"
         "|  \\\\    /   O peration     | Version:  " << FOAMversion << spaces << "|\n"
-        "|   \\\\  /    A nd           | Revision: " << FOAMDevRevisionNumber << spacesRev << "|\n"
+        "|   \\\\  /    A nd           | Web:      www.extend-project.de                 |\n"
-        "|    \\\\/     M anipulation  | Web:             www.OpenFOAM.org               |\n"
+        "|    \\\\/     M anipulation  |                                                 |\n"
         "\\*---------------------------------------------------------------------------*/\n";
 
     return os;

src/OpenFOAM/db/dlLibraryTable/dlLibraryTable.H

@@ -49,28 +49,6 @@ namespace Foam
                        Class dlLibraryTable Declaration
 \*---------------------------------------------------------------------------*/
 
-//- A means of hashing pointer addresses
-template<>
-class Hash<void*>
-{
-
-public:
-
-    Hash()
-    {}
-
-    long operator()(const void* const& p) const
-    {
-        return long(p);
-    }
-
-    label operator()(const void* const& p, const label tableSize) const
-    {
-        return abs(operator()(p)) % tableSize;
-    }
-};
-
-
 class dlLibraryTable
 :
     public HashTable<fileName, void*, Hash<void*> >

src/OpenFOAM/db/regIOobject/regIOobjectWrite.C

@@ -84,7 +84,15 @@ bool Foam::regIOobject::writeObject
 
     {
         // Try opening an OFstream for object
-        OFstream os(objectPath(), fmt, ver, cmp);
+        // Stream open for over-write.  HJ, 17/Aug/2010
+        OFstream os
+        (
+            objectPath(),
+            ios_base::out|ios_base::trunc,
+            fmt,
+            ver,
+            cmp
+        );
 
         // If any of these fail, return (leave error handling to Ostream class)
         if (!os.good())

src/OpenFOAM/fields/Fields/Field/SubFieldI.H

@@ -42,7 +42,7 @@ inline Foam::SubField<Type>::SubField
     const UList<Type>& list
 )
 :
-q    SubList<Type>(list, list.size())
+    SubList<Type>(list, list.size())
 {}
 
 

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -28,13 +28,13 @@ License
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-void Foam::Matrix<T>::allocate()
+void Foam::Matrix<Form, Type>::allocate()
 {
     if (n_ && m_)
     {
-        v_ = new T*[n_];
+        v_ = new Type*[n_];
-        v_[0] = new T[n_*m_];
+        v_[0] = new Type[n_*m_];
 
         for (register label i=1; i<n_; i++)
         {
@@ -46,8 +46,8 @@ void Foam::Matrix<T>::allocate()
 
 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T>::~Matrix()
+Foam::Matrix<Form, Type>::~Matrix()
 {
     if (v_)
     {
@@ -59,34 +59,8 @@ Foam::Matrix<T>::~Matrix()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-const Foam::Matrix<T>& Foam::Matrix<T>::null()
+Foam::Matrix<Form, Type>::Matrix(const label n, const label m)
-{
-    Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
-    return *nullPtr;
-}
-
-
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m)
-:
-    v_(NULL),
-    n_(n),
-    m_(m)
-{
-    if (n_ < 0 || m_ < 0)
-    {
-        FatalErrorIn("Matrix<T>::Matrix(const label n, const label m)")
-            << "bad n, m " << n_ << ", " << m_
-            << abort(FatalError);
-    }
-
-    allocate();
-}
-
-
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
 :
     v_(NULL),
     n_(n),
@@ -96,7 +70,28 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
     {
         FatalErrorIn
         (
-            "Matrix<T>::Matrix(const label n, const label m, const T&)"
+            "Matrix<Form, Type>::Matrix(const label n, const label m)"
+        )   << "bad n, m " << n_ << ", " << m_
+            << abort(FatalError);
+    }
+
+    allocate();
+}
+
+
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
+:
+    v_(NULL),
+    n_(n),
+    m_(m)
+{
+    if (n_ < 0 || m_ < 0)
+    {
+        FatalErrorIn
+        (
+            "Matrix<Form, Type>::Matrix"
+            "(const label n, const label m, const T&)"
         )   << "bad n, m " << n_ << ", " << m_
             << abort(FatalError);
     }
@@ -105,7 +100,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
 
     if (v_)
     {
-        T* v = v_[0];
+        Type* v = v_[0];
 
         label nm = n_*m_;
 
@@ -117,8 +112,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T>::Matrix(const Matrix<T>& a)
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
 :
     v_(NULL),
     n_(a.n_),
@@ -127,8 +122,8 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
     if (a.v_)
     {
         allocate();
-        T* v = v_[0];
+        Type* v = v_[0];
-        const T* av = a.v_[0];
+        const Type* av = a.v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -139,8 +134,8 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
 }
 
 
-template<class T>
+template<class Form, class Type>
-void Foam::Matrix<T>::clear()
+void Foam::Matrix<Form, Type>::clear()
 {
     if (v_)
     {
@@ -153,8 +148,8 @@ void Foam::Matrix<T>::clear()
 }
 
 
-template<class T>
+template<class Form, class Type>
-void Foam::Matrix<T>::transfer(Matrix<T>& a)
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
 {
     clear();
 
@@ -169,14 +164,32 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
 }
 
 
+template<class Form, class Type>
+Form Foam::Matrix<Form, Type>::T() const
+{
+    const Matrix<Form, Type>& A = *this;
+    Form At(m(), n());
+
+    for (register label i=0; i<n(); i++)
+    {
+        for (register label j=0; j<m(); j++)
+        {
+            At[j][i] = A[i][j];
+        }
+    }
+
+    return At;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-void Foam::Matrix<T>::operator=(const T& t)
+void Foam::Matrix<Form, Type>::operator=(const Type& t)
 {
     if (v_)
     {
-        T* v = v_[0];
+        Type* v = v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -188,12 +201,12 @@ void Foam::Matrix<T>::operator=(const T& t)
 
 
 // Assignment operator. Takes linear time.
-template<class T>
+template<class Form, class Type>
-void Foam::Matrix<T>::operator=(const Matrix<T>& a)
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 {
     if (this == &a)
     {
-        FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
+        FatalErrorIn("Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)")
             << "attempted assignment to self"
             << abort(FatalError);
     }
@@ -208,8 +221,8 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
 
     if (v_)
     {
-        T* v = v_[0];
+        Type* v = v_[0];
-        const T* av = a.v_[0];
+        const Type* av = a.v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -222,15 +235,15 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-const T& Foam::max(const Matrix<T>& a)
+const Type& Foam::max(const Matrix<Form, Type>& a)
 {
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
 
     if (nm)
     {
         label curMaxI = 0;
-        const T* v = a.v_[0];
+        const Type* v = a[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -244,7 +257,7 @@ const T& Foam::max(const Matrix<T>& a)
     }
     else
     {
-        FatalErrorIn("max(const Matrix<T>&)")
+        FatalErrorIn("max(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -254,15 +267,15 @@ const T& Foam::max(const Matrix<T>& a)
 }
 
 
-template<class T>
+template<class Form, class Type>
-const T& Foam::min(const Matrix<T>& a)
+const Type& Foam::min(const Matrix<Form, Type>& a)
 {
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
 
     if (nm)
     {
         label curMinI = 0;
-        const T* v = a.v_[0];
+        const Type* v = a[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -276,7 +289,7 @@ const T& Foam::min(const Matrix<T>& a)
     }
     else
     {
-        FatalErrorIn("min(const Matrix<T>&)")
+        FatalErrorIn("min(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -288,17 +301,17 @@ const T& Foam::min(const Matrix<T>& a)
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
+Form Foam::operator-(const Matrix<Form, Type>& a)
 {
-    Matrix<T> na(a.n_, a.m_);
+    Form na(a.n(), a.m());
 
-    if (a.n_ && a.m_)
+    if (a.n() && a.m())
     {
-        T* nav = na.v_[0];
+        Type* nav = na[0];
-        const T* av = a.v_[0];
+        const Type* av = a[0];
 
-        label nm = a.n_*a.m_;
+        label nm = a.n()*a.m();
         for (register label i=0; i<nm; i++)
         {
             nav[i] = -av[i];
@@ -309,32 +322,36 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
+Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
-    if (a.n_ != b.n_)
+    if (a.n() != b.n())
     {
-        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
+        FatalErrorIn
-            << "attempted add matrices with different number of rows: "
+        (
-            << a.n_ << ", " << b.n_
+            "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of rows: "
+            << a.n() << ", " << b.n()
             << abort(FatalError);
     }
 
-    if (a.m_ != b.m_)
+    if (a.m() != b.m())
     {
-        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
+        FatalErrorIn
-            << "attempted add matrices with different number of columns: "
+        (
-            << a.m_ << ", " << b.m_
+            "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of columns: "
+            << a.m() << ", " << b.m()
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Form ab(a.n(), a.m());
 
-    T* abv = ab.v_[0];
+    Type* abv = ab[0];
-    const T* av = a.v_[0];
+    const Type* av = a[0];
-    const T* bv = b.v_[0];
+    const Type* bv = b[0];
 
-    label nm = a.n_*a.m_;;
+    label nm = a.n()*a.m();
     for (register label i=0; i<nm; i++)
     {
         abv[i] = av[i] + bv[i];
@@ -344,32 +361,36 @@ Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
+Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
-    if (a.n_ != b.n_)
+    if (a.n() != b.n())
     {
-        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
+        FatalErrorIn
-            << "attempted add matrices with different number of rows: "
+        (
-            << a.n_ << ", " << b.n_
+            "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of rows: "
+            << a.n() << ", " << b.n()
             << abort(FatalError);
     }
 
-    if (a.m_ != b.m_)
+    if (a.m() != b.m())
     {
-        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
+        FatalErrorIn
-            << "attempted add matrices with different number of columns: "
+        (
-            << a.m_ << ", " << b.m_
+            "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of columns: "
+            << a.m() << ", " << b.m()
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Form ab(a.n(), a.m());
 
-    T* abv = ab.v_[0];
+    Type* abv = ab[0];
-    const T* av = a.v_[0];
+    const Type* av = a[0];
-    const T* bv = b.v_[0];
+    const Type* bv = b[0];
 
-    label nm = a.n_*a.m_;;
+    label nm = a.n()*a.m();
     for (register label i=0; i<nm; i++)
     {
         abv[i] = av[i] - bv[i];
@@ -379,17 +400,17 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T> Foam::operator*(const scalar s, const Matrix<T>& a)
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 {
-    Matrix<T> sa(a.n_, a.m_);
+    Form sa(a.n(), a.m());
 
-    if (a.n_ && a.m_)
+    if (a.n() && a.m())
     {
-        T* sav = sa.v_[0];
+        Type* sav = sa[0];
-        const T* av = a.v_[0];
+        const Type* av = a[0];
 
-        label nm = a.n_*a.m_;;
+        label nm = a.n()*a.m();
         for (register label i=0; i<nm; i++)
         {
             sav[i] = s*av[i];

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -51,31 +51,32 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class T> class Matrix;
+template<class Form, class Type> class Matrix;
 
-template<class T> const T& max(const Matrix<T>&);
+template<class Form, class Type> Istream& operator>>
-template<class T> const T& min(const Matrix<T>&);
+(
+    Istream&,
+    Matrix<Form, Type>&
+);
 
-template<class T> Matrix<T> operator-(const Matrix<T>&);
+template<class Form, class Type> Ostream& operator<<
-template<class T> Matrix<T> operator+(const Matrix<T>&, const Matrix<T>&);
+(
-template<class T> Matrix<T> operator-(const Matrix<T>&, const Matrix<T>&);
+    Ostream&,
-template<class T> Matrix<T> operator*(const scalar, const Matrix<T>&);
+    const Matrix<Form, Type>&
-
+);
-template<class T> Istream& operator>>(Istream&, Matrix<T>&);
-template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
 
 
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class T>
+template<class Form, class Type>
 class Matrix
 {
     // Private data
 
         //- Row pointers
-        T** __restrict__ v_;
+        Type** __restrict__ v_;
 
         //- Number of rows and columns in Matrix.
         label n_, m_;
@@ -87,6 +88,12 @@ class Matrix
 
 public:
 
+    // Static Member Functions
+
+        //- Return a null Matrix
+        inline static const Matrix<Form, Type>& null();
+
+
     // Constructors
 
         //- Null constructor.
@@ -97,16 +104,16 @@ public:
 
         //- Construct with given number of rows and columns
         //  and value for all elements.
-        Matrix(const label n, const label m, const T&);
+        Matrix(const label n, const label m, const Type&);
 
         //- Copy constructor.
-        Matrix(const Matrix<T>&);
+    Matrix(const Matrix<Form, Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);
 
         //- Clone
-        inline autoPtr<Matrix<T> > clone() const;
+        inline autoPtr<Matrix<Form, Type> > clone() const;
 
 
     // Destructor
@@ -116,10 +123,6 @@ public:
 
     // Member functions
 
-        //- Return a null Matrix
-        static const Matrix<T>& null();
-
-
         // Access
 
             //- Return the number of rows
@@ -148,48 +151,72 @@ public:
 
             //- Transfer the contents of the argument Matrix into this Matrix
             //  and annull the argument Matrix.
-            void transfer(Matrix<T>&);
+            void transfer(Matrix<Form, Type>&);
+
+
+        //- Return the transpose of the matrix
+        Form T() const;
 
 
     // Member operators
 
         //- Return subscript-checked element of Matrix.
-        inline T* operator[](const label);
+        inline Type* operator[](const label);
 
         //- Return subscript-checked element of constant Matrix.
-        inline const T* operator[](const label) const;
+        inline const Type* operator[](const label) const;
 
         //- Assignment operator. Takes linear time.
-        void operator=(const Matrix<T>&);
+        void operator=(const Matrix<Form, Type>&);
 
         //- Assignment of all entries to the given value
-        void operator=(const T&);
+        void operator=(const Type&);
-
-
-    // Friend functions
-
-        friend const T& max<T>(const Matrix<T>&);
-        friend const T& min<T>(const Matrix<T>&);
-
-
-    // Friend operators
-
-        friend Matrix<T> operator-<T>(const Matrix<T>&);
-        friend Matrix<T> operator+<T>(const Matrix<T>&, const Matrix<T>&);
-        friend Matrix<T> operator-<T>(const Matrix<T>&, const Matrix<T>&);
-        friend Matrix<T> operator*<T>(const scalar, const Matrix<T>&);
 
 
     // IOstream operators
 
         //- Read Matrix from Istream, discarding contents of existing Matrix.
-        friend Istream& operator>> <T>(Istream&, Matrix<T>&);
+        friend Istream& operator>> <Form, Type>
+        (
+            Istream&,
+            Matrix<Form, Type>&
+        );
 
         // Write Matrix to Ostream.
-        friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
+        friend Ostream& operator<< <Form, Type>
+        (
+            Ostream&,
+            const Matrix<Form, Type>&
+        );
 };
 
 
+// Global functions and operators
+
+template<class Form, class Type> const Type& max(const Matrix<Form, Type>&);
+template<class Form, class Type> const Type& min(const Matrix<Form, Type>&);
+
+template<class Form, class Type> Form operator-(const Matrix<Form, Type>&);
+
+template<class Form, class Type> Form operator+
+(
+    const Matrix<Form, Type>&,
+    const Matrix<Form, Type>&
+);
+
+template<class Form, class Type> Form operator-
+(
+    const Matrix<Form, Type>&,
+    const Matrix<Form, Type>&
+);
+
+template<class Form, class Type> Form operator*
+(
+    const scalar,
+    const Matrix<Form, Type>&
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -26,8 +26,8 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-inline Foam::Matrix<T>::Matrix()
+inline Foam::Matrix<Form, Type>::Matrix()
 :
     v_(NULL),
     n_(0),
@@ -35,71 +35,75 @@ inline Foam::Matrix<T>::Matrix()
 {}
 
 
-template<class T>
+template<class Form, class Type>
-inline Foam::autoPtr<Foam::Matrix<T> > Foam::Matrix<T>::clone() const
+inline Foam::autoPtr<Foam::Matrix<Form, Type> >
+Foam::Matrix<Form, Type>::clone() const
 {
-    return autoPtr<Matrix<T> >(new Matrix<T>(*this));
+    return autoPtr<Matrix<Form, Type> >(new Matrix<Form, Type>(*this));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Form, class Type>
+inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
+{
+    return *reinterpret_cast< Matrix<Form, Type>* >(0);
+}
+
+
 //- Return the number of rows
-template<class T>
+template<class Form, class Type>
-inline Foam::label Foam::Matrix<T>::n() const
+inline Foam::label Foam::Matrix<Form, Type>::n() const
 {
     return n_;
 }
 
 
-//- Return the number of columns
+template<class Form, class Type>
-template<class T>
+inline Foam::label Foam::Matrix<Form, Type>::m() const
-inline Foam::label Foam::Matrix<T>::m() const
 {
     return m_;
 }
 
 
-//- Return the number of columns
+template<class Form, class Type>
-template<class T>
+inline Foam::label Foam::Matrix<Form, Type>::size() const
-inline Foam::label Foam::Matrix<T>::size() const
 {
     return n_*m_;
 }
 
 
-// Check index i is within valid range (0 ... n-1).
+template<class Form, class Type>
-template<class T>
+inline void Foam::Matrix<Form, Type>::checki(const label i) const
-inline void Foam::Matrix<T>::checki(const label i) const
 {
     if (!n_)
     {
-        FatalErrorIn("Matrix<T>::checki(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checki(const label)")
             << "attempt to access element from zero sized row"
             << abort(FatalError);
     }
     else if (i<0 || i>=n_)
     {
-        FatalErrorIn("Matrix<T>::checki(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checki(const label)")
             << "index " << i << " out of range 0 ... " << n_-1
             << abort(FatalError);
     }
 }
 
 
-// Check index j is within valid range (0 ... n-1).
+template<class Form, class Type>
-template<class T>
+inline void Foam::Matrix<Form, Type>::checkj(const label j) const
-inline void Foam::Matrix<T>::checkj(const label j) const
 {
     if (!m_)
     {
-        FatalErrorIn("Matrix<T>::checkj(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
             << "attempt to access element from zero sized column"
             << abort(FatalError);
     }
     else if (j<0 || j>=m_)
     {
-        FatalErrorIn("Matrix<T>::checkj(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
             << "index " << j << " out of range 0 ... " << m_-1
             << abort(FatalError);
     }
@@ -108,9 +112,8 @@ inline void Foam::Matrix<T>::checkj(const label j) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-// Return subscript-checked element access
+template<class Form, class Type>
-template<class T>
+inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
-inline T* Foam::Matrix<T>::operator[](const label i)
 {
 #   ifdef FULLDEBUG
     checki(i);
@@ -119,9 +122,8 @@ inline T* Foam::Matrix<T>::operator[](const label i)
 }
 
 
-// Return subscript-checked const element access
+template<class Form, class Type>
-template<class T>
+inline const Type* Foam::Matrix<Form, Type>::operator[](const label i) const
-inline const T* Foam::Matrix<T>::operator[](const label i) const
 {
 #   ifdef FULLDEBUG
     checki(i);

src/OpenFOAM/matrices/Matrix/MatrixIO.C

@@ -32,8 +32,8 @@ License
 
 // * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
 
-template<class T>
+template<class Form, class Type>
-Foam::Matrix<T>::Matrix(Istream& is)
+Foam::Matrix<Form, Type>::Matrix(Istream& is)
 :
     v_(NULL),
     n_(0),
@@ -43,17 +43,17 @@ Foam::Matrix<T>::Matrix(Istream& is)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
 {
     // Anull matrix
     M.clear();
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&)");
 
     token firstToken(is);
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&) : reading first token");
 
     if (firstToken.isLabel())
     {
@@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
         label nm = M.n_*M.m_;
 
         // Read list contents depending on data format
-        if (is.format() == IOstream::ASCII || !contiguous<T>())
+        if (is.format() == IOstream::ASCII || !contiguous<Type>())
         {
             // Read beginning of contents
             char listDelimiter = is.readBeginList("Matrix");
@@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
             if (nm)
             {
                 M.allocate();
-                T* v = M.v_[0];
+                Type* v = M.v_[0];
 
                 if (listDelimiter == token::BEGIN_LIST)
                 {
@@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
 
                             is.fatalCheck
                             (
-                                "operator>>(Istream&, Matrix<T>&) : "
+                                "operator>>(Istream&, Matrix<Form, Type>&) : "
                                 "reading entry"
                             );
                         }
@@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
                 }
                 else
                 {
-                    T element;
+                    Type element;
                     is >> element;
 
                     is.fatalCheck
                     (
-                        "operator>>(Istream&, Matrix<T>&) : "
+                        "operator>>(Istream&, Matrix<Form, Type>&) : "
                         "reading the single entry"
                     );
 
@@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
             if (nm)
             {
                 M.allocate();
-                T* v = M.v_[0];
+                Type* v = M.v_[0];
 
-                is.read(reinterpret_cast<char*>(v), nm*sizeof(T));
+                is.read(reinterpret_cast<char*>(v), nm*sizeof(Type));
 
                 is.fatalCheck
                 (
-                    "operator>>(Istream&, Matrix<T>&) : "
+                    "operator>>(Istream&, Matrix<Form, Type>&) : "
                     "reading the binary block"
                 );
             }
@@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
     }
     else
     {
-        FatalIOErrorIn("operator>>(Istream&, Matrix<T>&)", is)
+        FatalIOErrorIn("operator>>(Istream&, Matrix<Form, Type>&)", is)
             << "incorrect first token, expected <int>, found "
             << firstToken.info()
             << exit(FatalIOError);
@@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
 }
 
 
-template<class T>
+template<class Form, class Type>
-Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
 {
     label nm = M.n_*M.m_;
 
     os  << M.n() << token::SPACE << M.m();
 
     // Write list contents depending on data format
-    if (os.format() == IOstream::ASCII || !contiguous<T>())
+    if (os.format() == IOstream::ASCII || !contiguous<Type>())
     {
         if (nm)
         {
             bool uniform = false;
 
-            const T* v = M.v_[0];
+            const Type* v = M.v_[0];
 
-            if (nm > 1 && contiguous<T>())
+            if (nm > 1 && contiguous<Type>())
             {
                 uniform = true;
 
@@ -189,7 +189,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
             }
             // Fix: matrices smaller than 20x20 will be written square.
             // HJ, 22/Jan/2009
-            else if (nm < 400 && contiguous<T>())
+            else if (nm < 400 && contiguous<Type>())
             {
                 // Write size of list and start contents delimiter
                 os  << token::BEGIN_LIST;
@@ -248,7 +248,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
     {
         if (nm)
         {
-            os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(T));
+            os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(Type));
         }
     }
 

src/OpenFOAM/matrices/Matrix/tools/DenseMatrixTools.C

@@ -25,7 +25,7 @@ License
 Class
     DenseMatrixTools
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "DenseMatrixTools.H"
 #include "Swap.H"
@@ -36,7 +36,7 @@ Class
 template<class T>
 static void Foam::DenseMatrixTools::solve
 (
-    Matrix<T>& A,
+    SquareMatrix<T>& A,
     List<T>& x,
     List<T>& b
 )
@@ -107,12 +107,12 @@ static void Foam::DenseMatrixTools::solve
 }
 
 
-template<class T>
+template<class Form, class Type>
 static void Foam::DenseMatrixTools::qrDecompose
 (
     const label nCols,
-    FieldField<Field, T>& A,
+    FieldField<Field, Type>& A,
-    Matrix<T>& R
+    Matrix<Form, T>& R
 )
 {
     // Note: consider Arnoldi algorithm for speed-up.  HJ, 14/Sep/2006
@@ -128,9 +128,9 @@ static void Foam::DenseMatrixTools::qrDecompose
 
         for (label k = j + 1; k < nCols; k++)
         {
-            const Field<T>& Aj = A[j];
+            const Field<Type>& Aj = A[j];
-            Field<T>& Ak = A[k];
+            Field<Type>& Ak = A[k];
-            T& Rjk = R[j][k];
+            Type& Rjk = R[j][k];
 
             Rjk = gSumProd(Aj, Ak);
 

src/OpenFOAM/matrices/Matrix/tools/DenseMatrixTools.H

@@ -33,7 +33,7 @@ SourceFiles
 #ifndef DenseMatrixTools_H
 #define DenseMatrixTools_H
 
-#include "Matrix.H"
+#include "SquareMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,16 +47,16 @@ namespace Foam
 namespace DenseMatrixTools
 {
     //- Solve using Gauss Elimination with pivoting.  Destroys matrix and b
-    template<class T>
+    template<class Type>
-    static void solve(Matrix<T>& A, List<T>& x, List<T>& b);
+    static void solve(SquareMatrix<Type>& A, List<Type>& x, List<Type>& b);
 
     //- Q-R decomposition
-    template<class T>
+    template<class Form, class T>
     static void qrDecompose
     (
         const label nCols,
         FieldField<Field, T>& A,
-        Matrix<T>& R
+        Matrix<class Form, class Type>& R
     );
 };
 

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -250,6 +250,11 @@ public:
             //- Maximum number of iterations
             label maxIter_;
 
+
+    protected:
+
+        // Protected data
+
             //- Matrix reference
             const lduMatrix& matrix_;
 
@@ -451,7 +456,9 @@ public:
     //- Abstract base-class for lduMatrix smoothers
     class smoother
     {
-        // Private data
+    protected:
+
+        // Protected data
 
             //- Matrix reference
             const lduMatrix& matrix_;
@@ -585,7 +592,9 @@ public:
     //- Abstract base-class for lduMatrix preconditioners
     class preconditioner
     {
-        // Private data
+    protected:
+
+        // Protected data
 
             //- Matrix reference
             const lduMatrix& matrix_;

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixPreconditioner.C

@@ -84,7 +84,7 @@ Foam::lduPreconditioner::New
 
     const dictionary& controls = e.isDict() ? e.dict() : dictionary::null;
 
-    if (sol.matrix().symmetric())
+    if (matrix.symmetric())
     {
         symMatrixConstructorTable::iterator constructorIter =
             symMatrixConstructorTablePtr_->find(preconName);
@@ -121,7 +121,7 @@ Foam::lduPreconditioner::New
             )
         );
     }
-    else if (sol.matrix().asymmetric())
+    else if (matrix.asymmetric())
     {
         asymMatrixConstructorTable::iterator constructorIter =
             asymMatrixConstructorTablePtr_->find(preconName);

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSmoother.C

@@ -172,6 +172,8 @@ Foam::autoPtr<Foam::lduMatrix::smoother> Foam::lduMatrix::smoother::New
         )   << "cannot solve incomplete matrix, "
                "no diagonal or off-diagonal coefficient"
             << exit(FatalIOError);
+
+        return autoPtr<lduSmoother>(NULL);
     }
 }
 

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

@@ -48,7 +48,7 @@ Foam::autoPtr<Foam::lduMatrix::solver> Foam::lduMatrix::solver::New
     const dictionary& dict
 )
 {
-    word name(dict.lookup("solver"));
+    word solverName(dict.lookup("solver"));
 
     if (matrix.diagonal())
     {
@@ -243,10 +243,11 @@ Foam::lduMatrix::solver::solver
 
 void Foam::lduMatrix::solver::readControls()
 {
+    tolerance_ = dict_.lookupOrDefault<scalar>("tolerance", 1e-6);
+    relTolerance_ = dict_.lookupOrDefault<scalar>("relTol", 0);
+
     minIter_ = dict_.lookupOrDefault<label>("minIter", 0);
     maxIter_ = dict_.lookupOrDefault<label>("maxIter", 1000);
-    tolerance_ = dict_.lookupOrDefault<scalar>("tolerance", 1e-6);
-    relTol_    = dict_.lookupOrDefault<scalar>("relTol", 0);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C

@@ -38,21 +38,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::DICPreconditioner::DICPreconditioner
-(
-    const lduMatrix::solver& sol,
-    const dictionary&
-)
-:
-    lduMatrix::preconditioner(sol),
-    rD_(sol.matrix().diag())
-{
-    calcReciprocalD(rD_, sol.matrix());
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::DICPreconditioner::calcReciprocalD

src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C

@@ -84,7 +84,7 @@ Foam::GAMGPreconditioner::~GAMGPreconditioner()
 void Foam::GAMGPreconditioner::readControls()
 {
     GAMG_.readControls();
-    nVcycles_ = controlDict_.lookupOrDefault<label>("nVcycles", 2);
+    nVcycles_ = GAMG_.dict().lookupOrDefault<label>("nVcycles", 2);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H

@@ -62,7 +62,7 @@ class GAMGPreconditioner
         //- Number of V-cycles to perform
         label nVcycles_;
 
-        //- Read the control parameters from the controlDict_
+        //- Read the control parameters from the dictionary
         virtual void readControls();
 
 public:

src/OpenFOAM/matrices/lduMatrix/solvers/BICCG/BICCG.H

@@ -94,7 +94,7 @@ public:
             const FieldField<Field, scalar>& coupleBouCoeffs,
             const FieldField<Field, scalar>& coupleIntCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const dictionary& solverControls
+            const dictionary& dict
         );
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C

@@ -108,7 +108,7 @@ Foam::GAMGSolver::GAMGSolver
     nFinestSweeps_(2),
     scaleCorrection_(matrix.symmetric()),
     directSolveCoarsest_(false),
-    agglomeration_(GAMGAgglomeration::New(matrix_, dict())),
+    agglomeration_(GAMGAgglomeration::New(matrix_, dict)),
 
     matrixLevels_(agglomeration_.size()),
     interfaceLevels_(agglomeration_.size()),

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -75,7 +75,7 @@ Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
         PtrList<scalarField> coarseB;
 
         // Create the smoothers for all levels
-        PtrList<lduMatrix::smoother> smoothers;
+        PtrList<lduSmoother> smoothers;
 
         // Initialise the above data structures
         initVcycle(coarseCorrX, coarseB, smoothers);
@@ -115,7 +115,7 @@ Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
 
 void Foam::GAMGSolver::Vcycle
 (
-    const PtrList<lduMatrix::smoother>& smoothers,
+    const PtrList<lduSmoother>& smoothers,
     scalarField& x,
     const scalarField& b,
     scalarField& Ax,
@@ -352,7 +352,7 @@ void Foam::GAMGSolver::initVcycle
 (
     PtrList<scalarField>& coarseCorrX,
     PtrList<scalarField>& coarseB,
-    PtrList<lduMatrix::smoother>& smoothers
+    PtrList<lduSmoother>& smoothers
 ) const
 {
     coarseCorrX.setSize(matrixLevels_.size());
@@ -363,13 +363,13 @@ void Foam::GAMGSolver::initVcycle
     smoothers.set
     (
         0,
-        lduMatrix::smoother::New
+        lduSmoother::New
         (
-            dict().lookup("smoother"),
             matrix_,
             coupleBouCoeffs_,
             coupleIntCoeffs_,
-            interfaces_
+            interfaces_,
+            dict()
         )
     );
 
@@ -396,13 +396,13 @@ void Foam::GAMGSolver::initVcycle
         smoothers.set
         (
             leveli + 1,
-            lduMatrix::smoother::New
+            lduSmoother::New
             (
-                dict().lookup("smoother"),
                 matrixLevels_[leveli],
                 coupleLevelsBouCoeffs_[leveli],
                 coupleLevelsIntCoeffs_[leveli],
-                interfaceLevels_[leveli]
+                interfaceLevels_[leveli],
+                dict()
             )
         );
     }

src/OpenFOAM/matrices/lduMatrix/solvers/ICCG/ICCG.C

@@ -72,7 +72,7 @@ Foam::ICCG::ICCG
     const FieldField<Field, scalar>& coupleBouCoeffs,
     const FieldField<Field, scalar>& coupleIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& solverControls
+    const dictionary& dict
 )
 :
     PCG
@@ -82,7 +82,7 @@ Foam::ICCG::ICCG
         coupleBouCoeffs,
         coupleIntCoeffs,
         interfaces,
-        solverControls
+        dict
     )
 {}
 

src/OpenFOAM/matrices/lduMatrix/solvers/ICCG/ICCG.H

@@ -94,7 +94,7 @@ public:
             const FieldField<Field, scalar>& coupleBouCoeffs,
             const FieldField<Field, scalar>& coupleIntCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const dictionary& solverControls
+            const dictionary& dict
         );
 
         //- Construct from matrix components and tolerances

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C

@@ -46,7 +46,7 @@ Foam::PBiCG::PBiCG
     const FieldField<Field, scalar>& coupleBouCoeffs,
     const FieldField<Field, scalar>& coupleIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& solverControls
+    const dictionary& dict
 )
 :
     lduSolver
@@ -56,7 +56,7 @@ Foam::PBiCG::PBiCG
         coupleBouCoeffs,
         coupleIntCoeffs,
         interfaces,
-        solverControls
+        dict
     )
 {}
 
@@ -71,10 +71,10 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
 ) const
 {
     // --- Setup class containing solver performance data
-    lduMatrix::solverPerformance solverPerf
+    lduSolverPerformance solverPerf
     (
-        lduMatrix::preconditioner::getName(controlDict_) + typeName,
+        lduMatrix::preconditioner::getName(dict()) + typeName,
-        fieldName_
+        fieldName()
     );
 
 
@@ -206,7 +206,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
 
             for (register label cell=0; cell<nCells; cell++)
             {
-                psiPtr[cell] += alpha*pAPtr[cell];
+                xPtr[cell] += alpha*pAPtr[cell];
                 rAPtr[cell] -= alpha*wAPtr[cell];
                 rTPtr[cell] -= alpha*wTPtr[cell];
             }

src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C

@@ -46,7 +46,7 @@ Foam::PCG::PCG
     const FieldField<Field, scalar>& coupleBouCoeffs,
     const FieldField<Field, scalar>& coupleIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& solverControls
+    const dictionary& dict
 )
 :
     lduSolver
@@ -56,7 +56,7 @@ Foam::PCG::PCG
         coupleBouCoeffs,
         coupleIntCoeffs,
         interfaces,
-        solverControls
+        dict
     )
 {}
 
@@ -71,10 +71,10 @@ Foam::lduSolverPerformance Foam::PCG::solve
 ) const
 {
     // --- Setup class containing solver performance data
-    lduMatrix::solverPerformance solverPerf
+    lduSolverPerformance solverPerf
     (
-        lduMatrix::preconditioner::getName(controlDict_) + typeName,
+        lduMatrix::preconditioner::getName(dict()) + typeName,
-        fieldName_
+        fieldName()
     );
 
     register label nCells = x.size();
@@ -194,7 +194,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
 
             for (register label cell=0; cell<nCells; cell++)
             {
-                psiPtr[cell] += alpha*pAPtr[cell];
+                xPtr[cell] += alpha*pAPtr[cell];
                 rAPtr[cell] -= alpha*wAPtr[cell];
             }
 

src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H

@@ -77,7 +77,7 @@ public:
             const FieldField<Field, scalar>& coupleBouCoeffs,
             const FieldField<Field, scalar>& coupleIntCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const dictionary& solverControls
+            const dictionary& dict
         );
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C

@@ -43,7 +43,7 @@ Foam::diagonalSolver::diagonalSolver
     const FieldField<Field, scalar>& coupleBouCoeffs,
     const FieldField<Field, scalar>& coupleIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& solverControls
+    const dictionary& dict
 )
 :
     lduMatrix::solver
@@ -53,7 +53,7 @@ Foam::diagonalSolver::diagonalSolver
         coupleBouCoeffs,
         coupleIntCoeffs,
         interfaces,
-        solverControls
+        dict
     )
 {}
 
@@ -72,7 +72,7 @@ Foam::lduMatrix::solverPerformance Foam::diagonalSolver::solve
     return lduSolverPerformance
     (
         typeName,
-        fieldName_,
+        fieldName(),
         0,
         0,
         0,

src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H

@@ -76,7 +76,7 @@ public:
             const FieldField<Field, scalar>& coupleBouCoeffs,
             const FieldField<Field, scalar>& coupleIntCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const dictionary& solverControls
+            const dictionary& dict
         );
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C

@@ -49,7 +49,7 @@ Foam::smoothSolver::smoothSolver
     const FieldField<Field, scalar>& coupleBouCoeffs,
     const FieldField<Field, scalar>& coupleIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& solverControls
+    const dictionary& dict
 )
 :
     lduSolver
@@ -59,7 +59,7 @@ Foam::smoothSolver::smoothSolver
         coupleBouCoeffs,
         coupleIntCoeffs,
         interfaces,
-        solverControls
+        dict
     )
 {
     readControls();
@@ -71,7 +71,7 @@ Foam::smoothSolver::smoothSolver
 void Foam::smoothSolver::readControls()
 {
     lduSolver::readControls();
-    nSweeps_ = controlDict_.lookupOrDefault<label>("nSweeps", 1);
+    nSweeps_ = dict().lookupOrDefault<label>("nSweeps", 1);
 }
 
 
@@ -93,6 +93,7 @@ Foam::lduSolverPerformance Foam::smoothSolver::solve
             matrix_,
             coupleBouCoeffs_,
             coupleIntCoeffs_,
+            interfaces_,
             dict()
         );
 
@@ -140,6 +141,7 @@ Foam::lduSolverPerformance Foam::smoothSolver::solve
                 matrix_,
                 coupleBouCoeffs_,
                 coupleIntCoeffs_,
+                interfaces_,
                 dict()
             );
 

src/OpenFOAM/meshes/pointMesh/pointMesh.H

@@ -106,7 +106,7 @@ public:
         //  HJ, 16/Oct/2008
         label nPoints() const
         {
-            return mesh_.nPoints();
+            return mesh.GeoMesh<polyMesh>::mesh_.nPoints();
         }
 
         //- Return number of points
@@ -115,7 +115,7 @@ public:
         //  HJ, 16/Oct/2008
         label nCells() const
         {
-            return mesh_.nCells();
+            return mesh.GeoMesh<polyMesh>::mesh_.nCells();
         }
 
         //- Return reference to boundary mesh

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C

@@ -26,11 +26,6 @@ License
 
 #include "Map.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template
@@ -41,7 +36,8 @@ template
     class PointType
 >
 const bool
-PrimitivePatch<Face, FaceList, PointField, PointType>::nSquaredProjection_
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+nSquaredProjection_
 (
     debug::optimisationSwitch("nSquaredProjection", 0) > 0
 );

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C

@@ -27,7 +27,6 @@ License
 #include "PrimitivePatch.H"
 #include "HashSet.H"
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchCheck.C

@@ -30,7 +30,7 @@ Description
 #include "PrimitivePatch.H"
 #include "Map.H"
 #include "ListOps.H"
-
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -497,8 +497,4 @@ writeVTKNormals
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C

@@ -27,7 +27,6 @@ License
 #include "PrimitivePatch.H"
 #include "demandDrivenData.H"
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C

@@ -112,10 +112,8 @@ calcMeshData() const
         markedPoints.find(pointPatch[pointi])() = pointi;
     }
 
-    //- Unsorted version:
+    forAll(*this, faceI)
-    //DynamicList<label> meshPoints(2*this->size());
+    {
-    //forAll(*this, facei)
-    //{
         const Face& curPoints = this->operator[](faceI);
 
         forAll (curPoints, pointI)
@@ -124,20 +122,6 @@ calcMeshData() const
         }
     }
 
-    // Create the storage and store the meshPoints.  Mesh points are
-    // the ones marked by the usage loop above
-    meshPointsPtr_ = new labelList(markedPoints.toc());
-    labelList& pointPatch = *meshPointsPtr_;
-
-    // Sort the list to preserve compatibility with the old ordering
-    sort(pointPatch);
-
-    // For every point in map give it its label in mesh points
-    forAll (pointPatch, pointI)
-    {
-        markedPoints.find(pointPatch[pointI])() = pointI;
-    }
-
     // Create local faces. Note that we start off from copy of original face
     // list (even though vertices are overwritten below). This is done so
     // additional data gets copied (e.g. region number of labelledTri)

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -125,85 +125,6 @@ Foam::primitiveMesh::~primitiveMesh()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::primitiveMesh::calcPointOrder
-(
-    label& nInternalPoints,
-    labelList& oldToNew,
-    const faceList& faces,
-    const label nInternalFaces,
-    const label nPoints
-)
-{
-    // Internal points are points that are not used by a boundary face.
-
-    // Map from old to new position
-    oldToNew.setSize(nPoints);
-    oldToNew = -1;
-
-
-    // 1. Create compact addressing for boundary points. Start off by indexing
-    // from 0 inside oldToNew. (shifted up later on)
-
-    label nBoundaryPoints = 0;
-    for (label faceI = nInternalFaces; faceI < faces.size(); faceI++)
-    {
-        const face& f = faces[faceI];
-
-        forAll(f, fp)
-        {
-            label pointI = f[fp];
-
-            if (oldToNew[pointI] == -1)
-            {
-                oldToNew[pointI] = nBoundaryPoints++;
-            }
-        }
-    }
-
-    // Now we know the number of boundary and internal points
-
-    nInternalPoints = nPoints - nBoundaryPoints;
-
-    // Move the boundary addressing up
-    forAll(oldToNew, pointI)
-    {
-        if (oldToNew[pointI] != -1)
-        {
-            oldToNew[pointI] += nInternalPoints;
-        }
-    }
-
-
-    // 2. Compact the internal points. Detect whether internal and boundary
-    // points are mixed.
-
-    label internalPointI = 0;
-
-    bool ordered = true;
-
-    for (label faceI = 0; faceI < nInternalFaces; faceI++)
-    {
-        const face& f = faces[faceI];
-
-        forAll(f, fp)
-        {
-            label pointI = f[fp];
-
-            if (oldToNew[pointI] == -1)
-            {
-                if (pointI >= nInternalPoints)
-                {
-                    ordered = false;
-                }
-                oldToNew[pointI] = internalPointI++;
-            }
-        }
-    }
-
-    return ordered;
-}
-
-
 void Foam::primitiveMesh::reset
 (
     const label nPoints,
@@ -216,40 +137,14 @@ void Foam::primitiveMesh::reset
 
     nPoints_ = nPoints;
     nEdges_ = -1;
-    nInternal0Edges_ = -1;
-    nInternal1Edges_ = -1;
-    nInternalEdges_ = -1;
 
     nInternalFaces_ = nInternalFaces;
     nFaces_ = nFaces;
     nCells_ = nCells;
 
-    // Check if points are ordered
-    label nInternalPoints;
-    labelList pointMap;
-
-    bool isOrdered = calcPointOrder
-    (
-        nInternalPoints,
-        pointMap,
-        faces(),
-        nInternalFaces_,
-        nPoints_
-    );
-
-    if (isOrdered)
-    {
-        nInternalPoints_ = nInternalPoints;
-    }
-    else
-    {
-        nInternalPoints_ = -1;
-    }
-
     if (debug)
     {
         Pout<< "primitiveMesh::reset : mesh reset to"
-            << " nInternalPoints:" << nInternalPoints_
             << " nPoints:" << nPoints_
             << " nEdges:" << nEdges_
             << " nInternalFaces:" << nInternalFaces_

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -200,12 +200,15 @@ class primitiveMesh
             //- Calculate point-point addressing
             void calcPointPoints() const;
 
-            //- Calculate edges, pointEdges and faceEdges (if doFaceEdges=true)
+            //- Calculate edges, pointEdges and faceEdges
-            //  During edge calculation, a larger set of data is assembled.
+            void calcEdges(const bool) const;
+
+            //- During edge calculation, a larger set of data is assembled.
             // Create and destroy as a set, using clearOutEdges()
-            void calcEdges(const bool doFaceEdges) const;
             void clearOutEdges();
-            //- Helper: return (after optional creation) edge between two points
+
+            //- Helper:
+            //  return (after optional creation) edge between two points
             static label getEdge
             (
                 List<DynamicList<label> >&,
@@ -432,18 +435,8 @@ public:
                     const label nCells = -1
                 );
 
-                //- Helper function to calculate point ordering. Returns true
+            // Removed calcPointOrder - garbage.  HJ, 27/Aug/2010
-                //  if points already ordered, false and fills pointMap (old to
+
-                //  new). Map splits points into those not used by any boundary
-                //  face and those that are.
-                static bool calcPointOrder
-                (
-                    label& nInternalPoints,
-                    labelList& pointMap,
-                    const faceList&,
-                    const label nInternalFaces,
-                    const label nPoints
-                );
 
             // Return mesh connectivity
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

@@ -71,386 +71,232 @@ Foam::label Foam::primitiveMesh::getEdge
 }
 
 
-void Foam::primitiveMesh::calcEdges(const bool doFaceEdges) const
+void Foam::primitiveMesh::calcEdges() const
 {
     if (debug)
     {
-        Pout<< "primitiveMesh::calcEdges(const bool) : "
+        Pout<< "primitiveMesh::calcEdges() : "
-            << "calculating edges, pointEdges and optionally faceEdges"
+            << "calculating edges, pointEdges and faceEdges"
             << endl;
     }
 
     // It is an error to attempt to recalculate edges
     // if the pointer is already set
-    if ((edgesPtr_ || pePtr_) || (doFaceEdges && fePtr_))
+    if (edgesPtr_ || fePtr_)
     {
         FatalErrorIn("primitiveMesh::calcEdges(const bool) const")
-            << "edges or pointEdges or faceEdges already calculated"
+            << "edges or faceEdges already calculated"
             << abort(FatalError);
     }
     else
     {
+        // Note: replaced with the original and correct algorithm
+        // which preserved correct ordering for edges list
+        // Version 1.5.x, 1.5-dev were wrong, with failures in parallel
+        // point-based data exchange.  Fixed in 1.6-ext and
+        // consecutive versions
+        // HJ, 27/Aug/2010
+
         // ALGORITHM:
-        // Go through the faces list. Search pointEdges for existing edge.
+        // Go through the pointFace list.  Go through the list of faces for that
-        // If not found create edge and add to pointEdges.
+        // point and ask for edges.  If the edge has got the point in question
-        // In second loop sort edges in order of increasing neighbouring
+        // AND the second point in the edge is larger than the first, add the
-        // point.
+        // edge to the list.  At the same time, add the edge label to the list
-        // This algorithm replaces the one using pointFaces which used more
+        // of edges for the current face (faceEdges) and log the other face as
-        // allocations but less memory and was on practical cases
+        // the neighbour of this face.
-        // quite a bit slower.
 
-        const faceList& fcs = faces();
+        const faceList& f = faces();
 
-        // Size up lists
+        const labelListList& pf = pointFaces();
-        // ~~~~~~~~~~~~~
 
-        // Estimate pointEdges storage
+        fePtr_ = new labelListList(nFaces());
-        List<DynamicList<label> > pe(nPoints());
+        labelListList& fe = *fePtr_;
-        forAll(pe, pointI)
+
+        // count the maximum number of edges
+        label maxEdges = 0;
+
+        // create a list of edges for each face and store for efficiency
+        edgeListList edgesOfFace(nFaces());
+
+        forAll (f, faceI)
         {
-            pe[pointI].setCapacity(primitiveMesh::edgesPerPoint_);
+            edgesOfFace[faceI] = f[faceI].edges();
-        }
 
-        // Estimate edges storage
+            maxEdges += f[faceI].nEdges();
-        DynamicList<edge> es(pe.size()*primitiveMesh::edgesPerPoint_/2);
 
-        // Estimate faceEdges storage
+            labelList& curFE = fe[faceI];
-        if (doFaceEdges)
+
+            curFE.setSize(f[faceI].nEdges());
+
+            forAll (curFE, curFEI)
             {
-            fePtr_ = new labelListList(fcs.size());
+                curFE[curFEI] = -1;
-            labelListList& faceEdges = *fePtr_;
-            forAll(fcs, faceI)
-            {
-                faceEdges[faceI].setSize(fcs[faceI].size());
             }
         }
 
+        // EDGE CALCULATION
 
-        // Find consecutive face points in edge list
+        edgesPtr_ = new edgeList(maxEdges);
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        edgeList& e = *edgesPtr_;
+        label nEdges = 0;
 
-        // Edges on boundary faces
+        forAll (pf, pointI)
-        label nExtEdges = 0;
-        // Edges using no boundary point
-        nInternal0Edges_ = 0;
-        // Edges using 1 boundary point
-        label nInt1Edges = 0;
-        // Edges using two boundary points but not on boundary face:
-        // edges.size()-nExtEdges-nInternal0Edges_-nInt1Edges
-
-        // Ordering is different if points are ordered.
-        if (nInternalPoints_ == -1)
         {
-            // No ordering. No distinction between types.
+            const labelList& curFaces = pf[pointI];
-            forAll(fcs, faceI)
-            {
-                const face& f = fcs[faceI];
 
-                forAll(f, fp)
+            // create a list of labels to keep the neighbours that
-                {
+            // have already been added
-                    label pointI = f[fp];
+            DynamicList<label, edgesPerPoint_> addedNeighbours;
-                    label nextPointI = f[f.fcIndex(fp)];
+            DynamicList<DynamicList<label, edgesPerPoint_> >
+                faceGivingNeighbour;
+            DynamicList<DynamicList<label, edgesPerPoint_> >
+                edgeOfFaceGivingNeighbour;
 
-                    label edgeI = getEdge(pe, es, pointI, nextPointI);
+            forAll (curFaces, faceI)
-
-                    if (doFaceEdges)
             {
-                        (*fePtr_)[faceI][fp] = edgeI;
+                // get the edges
+                const edgeList& fEdges = edgesOfFace[curFaces[faceI]];
+
+                // for every edge
+                forAll(fEdges, edgeI)
+                {
+                    const edge& ends = fEdges[edgeI];
+
+                    // does the edge has got the point in question
+                    bool found = false;
+                    label secondPoint = -1;
+
+                    if (ends.start() == pointI)
+                    {
+                        found = true;
+                        secondPoint = ends.end();
                     }
-                }
-            }
-            // Assume all edges internal
-            nExtEdges = 0;
-            nInternal0Edges_ = es.size();
-            nInt1Edges = 0;
-        }
-        else
-        {
-            // 1. Do external faces first. This creates external edges.
-            for (label faceI = nInternalFaces_; faceI < fcs.size(); faceI++)
-            {
-                const face& f = fcs[faceI];
 
-                forAll(f, fp)
+                    if (ends.end() == pointI)
                     {
-                    label pointI = f[fp];
+                        found = true;
-                    label nextPointI = f[f.fcIndex(fp)];
+                        secondPoint = ends.start();
-
-                    label oldNEdges = es.size();
-                    label edgeI = getEdge(pe, es, pointI, nextPointI);
-
-                    if (es.size() > oldNEdges)
-                    {
-                        nExtEdges++;
                     }
-                    if (doFaceEdges)
+
+                    // if the edge has got the point and second label is larger
+                    // than first, it is a candidate for adding
+                    if (found && (secondPoint > pointI))
                     {
-                        (*fePtr_)[faceI][fp] = edgeI;
+                        // check if the edge has already been added
-                    }
+                        bool added = false;
+
+                        forAll (addedNeighbours, eopI)
+                        {
+                            if (secondPoint == addedNeighbours[eopI])
+                            {
+                                // Edge is already added. New face sharing it
+                                added = true;
+
+                                // Remember the face and edge giving neighbour
+                                faceGivingNeighbour[eopI].append
+                                    (curFaces[faceI]);
+
+                                edgeOfFaceGivingNeighbour[eopI].append(edgeI);
+
+                                break;
                             }
                         }
 
-            // 2. Do internal faces. This creates internal edges.
+                        // If not added, add the edge to the list
-            for (label faceI = 0; faceI < nInternalFaces_; faceI++)
+                        if (!added)
                         {
-                const face& f = fcs[faceI];
+                            addedNeighbours.append(secondPoint);
 
-                forAll(f, fp)
+                            // Remember the face and subShape giving neighbour
-                {
+                            faceGivingNeighbour(addedNeighbours.size() - 1)
-                    label pointI = f[fp];
+                                .append(curFaces[faceI]);
-                    label nextPointI = f[f.fcIndex(fp)];
+                            edgeOfFaceGivingNeighbour
-
+                                (addedNeighbours.size() - 1).append(edgeI);
-                    label oldNEdges = es.size();
-                    label edgeI = getEdge(pe, es, pointI, nextPointI);
-
-                    if (es.size() > oldNEdges)
-                    {
-                        if (pointI < nInternalPoints_)
-                        {
-                            if (nextPointI < nInternalPoints_)
-                            {
-                                nInternal0Edges_++;
-                            }
-                            else
-                            {
-                                nInt1Edges++;
-                            }
-                        }
-                        else
-                        {
-                            if (nextPointI < nInternalPoints_)
-                            {
-                                nInt1Edges++;
-                            }
-                            else
-                            {
-                                // Internal edge with two points on boundary
-                            }
-                        }
-                    }
-                    if (doFaceEdges)
-                    {
-                        (*fePtr_)[faceI][fp] = edgeI;
                         }
                     }
                 }
             }
 
+            // All edges for the current point found. Before adding them to the
+            // list, it is necessary to sort them in the increasing order of the
+            // neighbouring point.
 
-        // For unsorted meshes the edges will be in order of occurrence inside
+            // Make real list out of SLList to simplify the manipulation.
-        // the faces. For sorted meshes the first nExtEdges will be the external
+            // Also, make another list to "remember" how the original list was
-        // edges.
+            // reshuffled.
+            labelList shuffleList(addedNeighbours.size());
 
-        if (nInternalPoints_ != -1)
+            forAll (shuffleList, i)
             {
-            nInternalEdges_ = es.size()-nExtEdges;
+                shuffleList[i] = i;
-            nInternal1Edges_ = nInternal0Edges_+nInt1Edges;
-
-            //Pout<< "Edge overview:" << nl
-            //    << "    total number of edges           : " << es.size()
-            //    << nl
-            //    << "    boundary edges                  : " << nExtEdges
-            //    << nl
-            //    << "    edges using no boundary point   : "
-            //    << nInternal0Edges_
-            //    << nl
-            //    << "    edges using one boundary points : " << nInt1Edges
-            //   << nl
-            //    << "    edges using two boundary points : "
-            //    << es.size()-nExtEdges-nInternal0Edges_-nInt1Edges << endl;
             }
 
+            // Use a simple sort to sort the addedNeighbours list.
+            //  Other two lists mimic the same behaviour
+            label i, j, a, b;
 
-        // Like faces sort edges in order of increasing neigbouring point.
+            for (j = 1; j <= addedNeighbours.size() - 1; j++)
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // Automatically if points are sorted into internal and external points
-        // the edges will be sorted into
-        // - internal point to internal point
-        // - internal point to external point
-        // - external point to external point
-        // The problem is that the last one mixes external edges with internal
-        // edges with two points on the boundary.
-
-
-        // Map to sort into new upper-triangular order
-        labelList oldToNew(es.size());
-        if (debug)
             {
-            oldToNew = -1;
+                a = addedNeighbours[j];
+                b = shuffleList[j];
+
+                i = j - 1;
+
+                while (i >= 0 && addedNeighbours[i] > a)
+                {
+                    addedNeighbours[i + 1] = addedNeighbours[i];
+                    shuffleList[i + 1] = shuffleList[i];
+                    i--;
                 }
 
-        // start of edges with 0 boundary points
+                addedNeighbours[i + 1] = a;
-        label internal0EdgeI = 0;
+                shuffleList[i + 1] = b;
+            }
 
-        // start of edges with 1 boundary points
+            labelList reshuffleList(shuffleList.size());
-        label internal1EdgeI = nInternal0Edges_;
 
-        // start of edges with 2 boundary points
+            forAll(shuffleList, i)
-        label internal2EdgeI = nInternal1Edges_;
+            {
+                reshuffleList[shuffleList[i]] = i;
+            }
 
-        // start of external edges
+            // Reshuffle other lists
-        label externalEdgeI = nInternalEdges_;
 
+            labelListList fgn(faceGivingNeighbour.size());
 
-        // To sort neighbouring points in increasing order. Defined outside
+            forAll (faceGivingNeighbour, i)
-        // for optimisation reasons: if all connectivity size same will need
+            {
-        // no reallocations
+                fgn[reshuffleList[i]].transfer(faceGivingNeighbour[i].shrink());
-        SortableList<label> nbrPoints(primitiveMesh::edgesPerPoint_);
+            }
 
-        forAll(pe, pointI)
+            labelListList eofgn(edgeOfFaceGivingNeighbour.size());
-        {
-            const DynamicList<label>& pEdges = pe[pointI];
 
-            nbrPoints.setSize(pEdges.size());
+            forAll (edgeOfFaceGivingNeighbour, i)
+            {
+                eofgn[reshuffleList[i]].transfer
+                (
+                    edgeOfFaceGivingNeighbour[i].shrink()
+                );
+            }
 
-            forAll(pEdges, i)
+            // adding the edges
+            forAll(addedNeighbours, edgeI)
             {
-                const edge& e = es[pEdges[i]];
+                const labelList& curFgn = fgn[edgeI];
+                const labelList& curEofgn = eofgn[edgeI];
 
-                label nbrPointI = e.otherVertex(pointI);
+                forAll (curFgn, fgnI)
+                {
+                    fe[curFgn[fgnI]][curEofgn[fgnI]] = nEdges;
+                }
 
-                if (nbrPointI < pointI)
+                e[nEdges] = edge(pointI, addedNeighbours[edgeI]);
-                {
+                nEdges++;
-                    nbrPoints[i] = -1;
-                }
-                else
-                {
-                    nbrPoints[i] = nbrPointI;
-                }
-            }
-            nbrPoints.sort();
-
-
-            if (nInternalPoints_ == -1)
-            {
-                // Sort all upper-triangular
-                forAll(nbrPoints, i)
-                {
-                    if (nbrPoints[i] != -1)
-                    {
-                        label edgeI = pEdges[nbrPoints.indices()[i]];
-
-                        oldToNew[edgeI] = internal0EdgeI++;
-                    }
-                }
-            }
-            else
-            {
-                if (pointI < nInternalPoints_)
-                {
-                    forAll(nbrPoints, i)
-                    {
-                        label nbrPointI = nbrPoints[i];
-
-                        label edgeI = pEdges[nbrPoints.indices()[i]];
-
-                        if (nbrPointI != -1)
-                        {
-                            if (edgeI < nExtEdges)
-                            {
-                                // External edge
-                                oldToNew[edgeI] = externalEdgeI++;
-                            }
-                            else if (nbrPointI < nInternalPoints_)
-                            {
-                                // Both points inside
-                                oldToNew[edgeI] = internal0EdgeI++;
-                            }
-                            else
-                            {
-                                // One points inside, one outside
-                                oldToNew[edgeI] = internal1EdgeI++;
-                            }
-                        }
-                    }
-                }
-                else
-                {
-                    forAll(nbrPoints, i)
-                    {
-                        label nbrPointI = nbrPoints[i];
-
-                        label edgeI = pEdges[nbrPoints.indices()[i]];
-
-                        if (nbrPointI != -1)
-                        {
-                            if (edgeI < nExtEdges)
-                            {
-                                // External edge
-                                oldToNew[edgeI] = externalEdgeI++;
-                            }
-                            else if (nbrPointI < nInternalPoints_)
-                            {
-                                // Not possible!
-                                FatalErrorIn("primitiveMesh::calcEdges(..)")
-                                    << abort(FatalError);
-                            }
-                            else
-                            {
-                                // Both points outside
-                                oldToNew[edgeI] = internal2EdgeI++;
-                            }
-                        }
-                    }
-                }
             }
         }
 
-
+        // reset the size
-        if (debug)
+        e.setSize(nEdges);
-        {
-            label edgeI = findIndex(oldToNew, -1);
-
-            if (edgeI != -1)
-            {
-                const edge& e = es[edgeI];
-
-                FatalErrorIn("primitiveMesh::calcEdges(..)")
-                    << "Did not sort edge " << edgeI << " points:" << e
-                    << " coords:" << points()[e[0]] << points()[e[1]]
-                    << endl
-                    << "Current buckets:" << endl
-                    << "    internal0EdgeI:" << internal0EdgeI << endl
-                    << "    internal1EdgeI:" << internal1EdgeI << endl
-                    << "    internal2EdgeI:" << internal2EdgeI << endl
-                    << "    externalEdgeI:" << externalEdgeI << endl
-                    << abort(FatalError);
-            }
-        }
-
-
-
-        // Renumber and transfer.
-
-        // Edges
-        edgesPtr_ = new edgeList(es.size());
-        edgeList& edges = *edgesPtr_;
-        forAll(es, edgeI)
-        {
-            edges[oldToNew[edgeI]] = es[edgeI];
-        }
-
-        // pointEdges
-        pePtr_ = new labelListList(nPoints());
-        labelListList& pointEdges = *pePtr_;
-        forAll(pe, pointI)
-        {
-            DynamicList<label>& pEdges = pe[pointI];
-            pEdges.shrink();
-            inplaceRenumber(oldToNew, pEdges);
-            pointEdges[pointI].transfer(pEdges);
-            Foam::sort(pointEdges[pointI]);
-        }
-
-        // faceEdges
-        if (doFaceEdges)
-        {
-            labelListList& faceEdges = *fePtr_;
-            forAll(faceEdges, faceI)
-            {
-                inplaceRenumber(oldToNew, faceEdges[faceI]);
-            }
-        }
     }
 }
 
@@ -505,70 +351,22 @@ const Foam::edgeList& Foam::primitiveMesh::edges() const
 {
     if (!edgesPtr_)
     {
-        //calcEdges(true);
+        // 1.6.x merge: reverting to correct code
-        calcEdges(false);
+        // HJ, 17/Aug/2010
+        calcEdges();
     }
 
     return *edgesPtr_;
 }
 
-const Foam::labelListList& Foam::primitiveMesh::pointEdges() const
-{
-    if (!pePtr_)
-    {
-        //calcEdges(true);
-        calcEdges(false);
-    }
-
-    return *pePtr_;
-}
-
 
 const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
 {
     if (!fePtr_)
     {
-        if (debug)
+        // 1.6.x merge: reverting to correct code
-        {
+        // HJ, 17/Aug/2010
-            Pout<< "primitiveMesh::faceEdges() : "
+        calcEdges();
-                << "calculating faceEdges" << endl;
-        }
-
-        //calcEdges(true);
-        const faceList& fcs = faces();
-        const labelListList& pe = pointEdges();
-        const edgeList& es = edges();
-
-        fePtr_ = new labelListList(fcs.size());
-        labelListList& faceEdges = *fePtr_;
-
-        forAll(fcs, faceI)
-        {
-            const face& f = fcs[faceI];
-
-            labelList& fEdges = faceEdges[faceI];
-            fEdges.setSize(f.size());
-
-            forAll(f, fp)
-            {
-                label pointI = f[fp];
-                label nextPointI = f[f.fcIndex(fp)];
-
-                // Find edge between pointI, nextPontI
-                const labelList& pEdges = pe[pointI];
-
-                forAll(pEdges, i)
-                {
-                    label edgeI = pEdges[i];
-
-                    if (es[edgeI].otherVertex(pointI) == nextPointI)
-                    {
-                        fEdges[fp] = edgeI;
-                        break;
-                    }
-                }
-            }
-        }
     }
 
     return *fePtr_;
@@ -578,7 +376,6 @@ const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
 void Foam::primitiveMesh::clearOutEdges()
 {
     deleteDemandDrivenData(edgesPtr_);
-    deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(fePtr_);
     labels_.clear();
     labelSet_.clear();
@@ -678,6 +475,4 @@ const Foam::labelList& Foam::primitiveMesh::cellEdges(const label cellI) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H

@@ -33,42 +33,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label primitiveMesh::nInternalPoints() const
-{
-    return nInternalPoints_;
-}
-
-
 inline label primitiveMesh::nPoints() const
 {
     return nPoints_;
 }
 
 
-inline label primitiveMesh::nInternal0Edges() const
-{
-    // Force edge calculation
-    nEdges();
-    return nInternal0Edges_;
-}
-
-
-inline label primitiveMesh::nInternal1Edges() const
-{
-    // Force edge calculation
-    nEdges();
-    return nInternal1Edges_;
-}
-
-
-inline label primitiveMesh::nInternalEdges() const
-{
-    // Force edge calculation
-    nEdges();
-    return nInternalEdges_;
-}
-
-
 inline label primitiveMesh::nEdges() const
 {
     if (nEdges_ < 0)

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -624,7 +624,6 @@ inline pointHit triangle<Point, PointRef>::fastIntersection
             // Ray parallel to triangle. Return miss
             return intersection;
         }
-    }
 
         const scalar inv_det = 1/det;
 

src/OpenFOAM/primitives/Scalar/Scalar.H

@@ -208,6 +208,11 @@ inline Scalar cmptDivide(const Scalar s1, const Scalar s2)
     return s1/s2;
 }
 
+inline Scalar cmptSumMultiply(const Scalar s1, const Scalar s2)
+{
+    return s1*s2;
+}
+
 inline Scalar cmptMax(const Scalar s)
 {
     return s;
@@ -275,6 +280,24 @@ inline Scalar stabilise(const Scalar s, const Scalar small)
 }
 
 
+inline Scalar cmptStabilise
+(
+    const Scalar s,
+    const Scalar small,
+    const Scalar value
+)
+{
+    if (mag(s) < small)
+    {
+        return sign(s)*value;
+    }
+    else
+    {
+        return s;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Scalar readScalar(Istream&);