From 26f2c24d84d653169811d752b5ff3884c2c9c87a Mon Sep 17 00:00:00 2001 From: henry Date: Sat, 27 Sep 2008 09:25:30 +0100 Subject: [PATCH] Added support for quadraticFit interpolation scheme. Added the quadraticFit interpolation scheme. --- src/OpenFOAM/Make/files | 1 + .../matrices/DiagonalMatrix/DiagonalMatrix.C | 99 ++++ .../matrices/DiagonalMatrix/DiagonalMatrix.H | 102 ++++ .../LUscalarMatrix/LUscalarMatrixTemplates.C | 14 +- .../{containers => matrices}/Matrix/Matrix.C | 174 +++--- .../{containers => matrices}/Matrix/Matrix.H | 83 ++- .../{containers => matrices}/Matrix/MatrixI.H | 52 +- .../Matrix/MatrixIO.C | 44 +- src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C | 400 +++++++++++++ src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H | 128 +++++ src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H | 75 +++ .../matrices/scalarMatrix/scalarMatrix.C | 159 +++++- .../matrices/scalarMatrix/scalarMatrix.H | 46 +- .../scalarMatrix/scalarMatrixTemplates.C | 26 +- .../matrices/simpleMatrix/simpleMatrix.C | 68 +-- .../matrices/simpleMatrix/simpleMatrix.H | 66 +-- src/finiteVolume/Make/files | 5 + .../fvMesh/extendedStencil/extendedStencil.C | 525 ++++++++++++++++++ .../fvMesh/extendedStencil/extendedStencil.H | 151 +++++ .../extendedStencilTemplates.C | 129 +++++ .../schemes/quadraticFit/quadraticFit.C | 36 ++ .../schemes/quadraticFit/quadraticFit.H | 138 +++++ .../schemes/quadraticFit/quadraticFitData.C | 310 +++++++++++ .../schemes/quadraticFit/quadraticFitData.H | 140 +++++ 24 files changed, 2722 insertions(+), 249 deletions(-) create mode 100644 src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C create mode 100644 src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H rename src/OpenFOAM/{containers => matrices}/Matrix/Matrix.C (62%) rename src/OpenFOAM/{containers => matrices}/Matrix/Matrix.H (70%) rename src/OpenFOAM/{containers => matrices}/Matrix/MatrixI.H (66%) rename src/OpenFOAM/{containers => matrices}/Matrix/MatrixIO.C (84%) create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 7b137dedc0..57dc7f8dc7 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -179,6 +179,7 @@ matrices/solution/solution.C scalarMatrix = matrices/scalarMatrix $(scalarMatrix)/scalarMatrix.C +$(scalarMatrix)/SVD/SVD.C LUscalarMatrix = matrices/LUscalarMatrix $(LUscalarMatrix)/LUscalarMatrix.C diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C new file mode 100644 index 0000000000..c023e59c5a --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -0,0 +1,99 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 + +\*---------------------------------------------------------------------------*/ + +#include "DiagonalMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +: + List(min(a.n(), a.m())) +{ + forAll(*this, i) + { + this->operator[](i) = a[i][i]; + } +} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size) +: + List(size) +{} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size, const Type& val) +: + List(size, val) +{} + + +template +Foam::DiagonalMatrix& Foam::DiagonalMatrix::invert() +{ + forAll(*this, i) + { + Type x = this->operator[](i); + if (mag(x) < VSMALL) + { + this->operator[](i) = Type(0); + } + else + { + this->operator[](i) = Type(1)/x; + } + } + + return this; +} + + +template +Foam::DiagonalMatrix Foam::inv(const DiagonalMatrix& A) +{ + DiagonalMatrix Ainv = A; + + forAll(A, i) + { + Type x = A[i]; + if (mag(x) < VSMALL) + { + Ainv[i] = Type(0); + } + else + { + Ainv[i] = Type(1)/x; + } + } + + return Ainv; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H new file mode 100644 index 0000000000..36138fc25d --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 + Foam::DiagonalMatrix + +Description + DiagonalMatrix is a 2D diagonal matrix of objects + of type Type, size nxn + +SourceFiles + DiagonalMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef DiagonalMatrix_H +#define DiagonalMatrix_H + +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // + +template class Matrix; + +/*---------------------------------------------------------------------------*\ + Class DiagonalMatrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class DiagonalMatrix +: + public List +{ +public: + + // Constructors + + //- Construct from diagonal component of a Matrix + DiagonalMatrix(const Matrix&); + + //- Construct empty from size + DiagonalMatrix(const label size); + + //- Construct from size and a value + DiagonalMatrix(const label, const Type&); + + + // Member functions + + //- Invert the diaganol matrix and return itself + DiagonalMatrix& invert(); +}; + + +// Global functions + +//- Return the diagonal Matrix inverse +template +DiagonalMatrix inv(const DiagonalMatrix&); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "DiagonalMatrix.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C index cbb59e2c86..ac0bc17207 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C @@ -28,16 +28,16 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::LUscalarMatrix::solve(Field& sourceSol) const +template +void Foam::LUscalarMatrix::solve(Field& sourceSol) const { if (Pstream::parRun()) { - Field completeSourceSol(n()); + Field completeSourceSol(n()); if (Pstream::master()) { - typename Field::subField + typename Field::subField ( completeSourceSol, sourceSol.size() @@ -58,7 +58,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } @@ -77,7 +77,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const { LUBacksubstitute(*this, pivotIndices_, completeSourceSol); - sourceSol = typename Field::subField + sourceSol = typename Field::subField ( completeSourceSol, sourceSol.size() @@ -98,7 +98,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } diff --git a/src/OpenFOAM/containers/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C similarity index 62% rename from src/OpenFOAM/containers/Matrix/Matrix.C rename to src/OpenFOAM/matrices/Matrix/Matrix.C index 733e04034e..4fbef71c2b 100644 --- a/src/OpenFOAM/containers/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -28,13 +28,13 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -template -void Foam::Matrix::allocate() +template +void Foam::Matrix::allocate() { if (n_ && m_) { - v_ = new T*[n_]; - v_[0] = new T[n_*m_]; + v_ = new Type*[n_]; + v_[0] = new Type[n_*m_]; for (register label i=1; i::allocate() // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // -template -Foam::Matrix::~Matrix() +template +Foam::Matrix::~Matrix() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } } @@ -59,16 +59,16 @@ Foam::Matrix::~Matrix() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -const Foam::Matrix& Foam::Matrix::null() +template +const Foam::Matrix& Foam::Matrix::null() { - Matrix* nullPtr = reinterpret_cast*>(NULL); + Matrix* nullPtr = reinterpret_cast*>(NULL); return *nullPtr; } -template -Foam::Matrix::Matrix(const label n, const label m) +template +Foam::Matrix::Matrix(const label n, const label m) : n_(n), m_(m), @@ -76,7 +76,7 @@ Foam::Matrix::Matrix(const label n, const label m) { if (n_ < 0 || m_ < 0) { - FatalErrorIn("Matrix::Matrix(const label n, const label m)") + FatalErrorIn("Matrix::Matrix(const label n, const label m)") << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -85,8 +85,8 @@ Foam::Matrix::Matrix(const label n, const label m) } -template -Foam::Matrix::Matrix(const label n, const label m, const T& a) +template +Foam::Matrix::Matrix(const label n, const label m, const Type& a) : n_(n), m_(m), @@ -96,7 +96,7 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) { FatalErrorIn ( - "Matrix::Matrix(const label n, const label m, const T&)" + "Matrix::Matrix(const label n, const label m, const T&)" ) << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -105,7 +105,7 @@ Foam::Matrix::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 +117,8 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) } -template -Foam::Matrix::Matrix(const Matrix& a) +template +Foam::Matrix::Matrix(const Matrix& a) : n_(a.n_), m_(a.m_), @@ -127,8 +127,8 @@ Foam::Matrix::Matrix(const Matrix& a) if (a.v_) { allocate(); - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::Matrix(const Matrix& a) } } - -template -void Foam::Matrix::clear() + +template +void Foam::Matrix::clear() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } n_ = 0; @@ -153,8 +153,8 @@ void Foam::Matrix::clear() } -template -void Foam::Matrix::transfer(Matrix& a) +template +void Foam::Matrix::transfer(Matrix& a) { clear(); @@ -171,12 +171,30 @@ void Foam::Matrix::transfer(Matrix& a) // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -void Foam::Matrix::operator=(const T& t) +template +Foam::Matrix Foam::Matrix::T() const +{ + const Matrix& A = *this; + Matrix At(m(), n()); + + for (register label i=0; i +void Foam::Matrix::operator=(const Type& t) { if (v_) { - T* v = v_[0]; + Type* v = v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const T& t) // Assignment operator. Takes linear time. -template -void Foam::Matrix::operator=(const Matrix& a) +template +void Foam::Matrix::operator=(const Matrix& a) { if (this == &a) { - FatalErrorIn("Matrix::operator=(const Matrix&)") + FatalErrorIn("Matrix::operator=(const Matrix&)") << "attempted assignment to self" << abort(FatalError); } @@ -204,12 +222,12 @@ void Foam::Matrix::operator=(const Matrix& a) n_ = a.n_; m_ = a.m_; allocate(); - } + } if (v_) { - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const Matrix& a) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // -template -const T& Foam::max(const Matrix& a) +template +const Type& Foam::max(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMaxI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("max(const Matrix&)") + FatalErrorIn("max(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -254,15 +272,15 @@ const T& Foam::max(const Matrix& a) } -template -const T& Foam::min(const Matrix& a) +template +const Type& Foam::min(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMinI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("min(const Matrix&)") + FatalErrorIn("min(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -288,15 +306,15 @@ const T& Foam::min(const Matrix& a) // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // -template -Foam::Matrix Foam::operator-(const Matrix& a) +template +Foam::Matrix Foam::operator-(const Matrix& a) { - Matrix na(a.n_, a.m_); + Matrix na(a.n_, a.m_); if (a.n_ && a.m_) { - T* nav = na.v_[0]; - const T* av = a.v_[0]; + Type* nav = na.v_[0]; + const Type* av = a.v_[0]; label nm = a.n_*a.m_; for (register label i=0; i Foam::operator-(const Matrix& a) } -template -Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) +template +Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Matrix ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator+(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) +template +Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Matrix ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator-(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator*(const scalar s, const Matrix& a) +template +Foam::Matrix Foam::operator*(const scalar s, const Matrix& a) { - Matrix sa(a.n_, a.m_); + Matrix sa(a.n_, a.m_); if (a.n_ && a.m_) { - T* sav = sa.v_[0]; - const T* av = a.v_[0]; + Type* sav = sa.v_[0]; + const Type* av = a.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i class Matrix; +template class Matrix; -template const T& max(const Matrix&); -template const T& min(const Matrix&); +template const Type& max(const Matrix&); +template const Type& min(const Matrix&); -template Matrix operator-(const Matrix&); -template Matrix operator+(const Matrix&, const Matrix&); -template Matrix operator-(const Matrix&, const Matrix&); -template Matrix operator*(const scalar, const Matrix&); +template Matrix operator-(const Matrix&); +template Matrix operator+ +( + const Matrix&, + const Matrix& +); +template Matrix operator- +( + const Matrix&, + const Matrix& +); +template Matrix operator* +( + const scalar, + const Matrix& +); -template Istream& operator>>(Istream&, Matrix&); -template Ostream& operator<<(Ostream&, const Matrix&); +template Istream& operator>>(Istream&, Matrix&); +template Ostream& operator<<(Ostream&, const Matrix&); /*---------------------------------------------------------------------------*\ Class Matrix Declaration \*---------------------------------------------------------------------------*/ -template +template class Matrix { // Private data @@ -78,7 +90,7 @@ class Matrix label n_, m_; //- Row pointers - T** __restrict__ v_; + Type** __restrict__ v_; //- Allocate the storage for the row-pointers and the data // and set the row pointers @@ -97,16 +109,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&); + Matrix(const Matrix&); //- Construct from Istream. Matrix(Istream&); //- Clone - inline autoPtr > clone() const; + inline autoPtr > clone() const; // Destructor @@ -117,7 +129,7 @@ public: // Member functions //- Return a null Matrix - static const Matrix& null(); + static const Matrix& null(); // Access @@ -148,45 +160,60 @@ public: //- Transfer the contents of the argument Matrix into this Matrix // and annull the argument Matrix. - void transfer(Matrix&); + void transfer(Matrix&); // 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; + + //- Return the transpose of the matrix + Matrix T() const; //- Assignment operator. Takes linear time. - void operator=(const Matrix&); + void operator=(const Matrix&); //- Assignment of all entries to the given value - void operator=(const T&); + void operator=(const Type&); // Friend functions - friend const T& max(const Matrix&); - friend const T& min(const Matrix&); + friend const Type& max(const Matrix&); + friend const Type& min(const Matrix&); // Friend operators - friend Matrix operator-(const Matrix&); - friend Matrix operator+(const Matrix&, const Matrix&); - friend Matrix operator-(const Matrix&, const Matrix&); - friend Matrix operator*(const scalar, const Matrix&); + friend Matrix operator-(const Matrix&); + friend Matrix operator+ + ( + const Matrix&, + const Matrix& + ); + friend Matrix operator- + ( + const Matrix&, + const Matrix& + ); + friend Matrix operator* + ( + const scalar, + const Matrix& + ); // IOstream operators //- Read Matrix from Istream, discarding contents of existing Matrix. - friend Istream& operator>> (Istream&, Matrix&); + friend Istream& operator>> (Istream&, Matrix&); // Write Matrix to Ostream. - friend Ostream& operator<< (Ostream&, const Matrix&); + friend Ostream& operator<< (Ostream&, const Matrix&); }; diff --git a/src/OpenFOAM/containers/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H similarity index 66% rename from src/OpenFOAM/containers/Matrix/MatrixI.H rename to src/OpenFOAM/matrices/Matrix/MatrixI.H index d0e3fe8fdf..8323174d4f 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -26,8 +26,8 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -inline Foam::Matrix::Matrix() +template +inline Foam::Matrix::Matrix() : n_(0), m_(0), @@ -35,71 +35,67 @@ inline Foam::Matrix::Matrix() {} -template -inline Foam::autoPtr > Foam::Matrix::clone() const +template +inline Foam::autoPtr > Foam::Matrix::clone() const { - return autoPtr >(new Matrix(*this)); + return autoPtr >(new Matrix(*this)); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // //- Return the number of rows -template -inline Foam::label Foam::Matrix::n() const +template +inline Foam::label Foam::Matrix::n() const { return n_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::m() const +template +inline Foam::label Foam::Matrix::m() const { return m_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::size() const +template +inline Foam::label Foam::Matrix::size() const { return n_*m_; } -// Check index i is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checki(const label i) const +template +inline void Foam::Matrix::checki(const label i) const { if (!n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "attempt to access element from zero sized row" << abort(FatalError); } else if (i<0 || i>=n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "index " << i << " out of range 0 ... " << n_-1 << abort(FatalError); } } -// Check index j is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checkj(const label j) const +template +inline void Foam::Matrix::checkj(const label j) const { if (!m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "attempt to access element from zero sized column" << abort(FatalError); } else if (j<0 || j>=m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "index " << j << " out of range 0 ... " << m_-1 << abort(FatalError); } @@ -108,9 +104,8 @@ inline void Foam::Matrix::checkj(const label j) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -// Return subscript-checked element access -template -inline T* Foam::Matrix::operator[](const label i) +template +inline Type* Foam::Matrix::operator[](const label i) { # ifdef FULLDEBUG checki(i); @@ -119,9 +114,8 @@ inline T* Foam::Matrix::operator[](const label i) } -// Return subscript-checked const element access -template -inline const T* Foam::Matrix::operator[](const label i) const +template +inline const Type* Foam::Matrix::operator[](const label i) const { # ifdef FULLDEBUG checki(i); diff --git a/src/OpenFOAM/containers/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C similarity index 84% rename from src/OpenFOAM/containers/Matrix/MatrixIO.C rename to src/OpenFOAM/matrices/Matrix/MatrixIO.C index 710e42876e..9b5a922f0d 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixIO.C +++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C @@ -32,8 +32,8 @@ License // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // -template -Foam::Matrix::Matrix(Istream& is) +template +Foam::Matrix::Matrix(Istream& is) : n_(0), m_(0), @@ -43,17 +43,17 @@ Foam::Matrix::Matrix(Istream& is) } -template -Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) +template +Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) { // Anull matrix M.clear(); - is.fatalCheck("operator>>(Istream&, Matrix&)"); + is.fatalCheck("operator>>(Istream&, Matrix&)"); token firstToken(is); - is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); + is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); if (firstToken.isLabel()) { @@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) label nm = M.n_*M.m_; // Read list contents depending on data format - if (is.format() == IOstream::ASCII || !contiguous()) + if (is.format() == IOstream::ASCII || !contiguous()) { // Read beginning of contents char listDelimiter = is.readBeginList("Matrix"); @@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& 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& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading entry" ); } @@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - T element; + Type element; is >> element; is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the single entry" ); @@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) if (nm) { M.allocate(); - T* v = M.v_[0]; + Type* v = M.v_[0]; - is.read(reinterpret_cast(v), nm*sizeof(T)); + is.read(reinterpret_cast(v), nm*sizeof(Type)); is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the binary block" ); } @@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) + FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) << "incorrect first token, expected , found " << firstToken.info() << exit(FatalIOError); @@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } -template -Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& 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()) + if (os.format() == IOstream::ASCII || !contiguous()) { if (nm) { bool uniform = false; - const T* v = M.v_[0]; + const Type* v = M.v_[0]; - if (nm > 1 && contiguous()) + if (nm > 1 && contiguous()) { uniform = true; @@ -187,7 +187,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) // Write end of contents delimiter os << token::END_BLOCK; } - else if (nm < 10 && contiguous()) + else if (nm < 10 && contiguous()) { // Write size of list and start contents delimiter os << token::BEGIN_LIST; @@ -246,7 +246,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { if (nm) { - os.write(reinterpret_cast(M.v_[0]), nm*sizeof(T)); + os.write(reinterpret_cast(M.v_[0]), nm*sizeof(Type)); } } diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C new file mode 100644 index 0000000000..a1355c393d --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C @@ -0,0 +1,400 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +\*---------------------------------------------------------------------------*/ + +#include "SVD.H" +#include "scalarList.H" +#include "scalarMatrix.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * // + +Foam::SVD::SVD(const Matrix& A, const scalar minCondition) +: + U_(A), + V_(A.m(), A.m()), + S_(A.m()), + VSinvUt_(A.m(), A.n()), + nZeros_(0) +{ + // SVDcomp to find U_, V_ and S_ - the singular values + + const label Um = U_.m(); + const label Un = U_.n(); + + scalarList rv1(Um); + scalar g = 0; + scalar scale = 0; + scalar s = 0; + scalar anorm = 0; + label l = 0; + + for (label i = 0; i < Um; i++) + { + l = i+2; + rv1[i] = scale*g; + g = s = scale = 0; + + if (i < Un) + { + for (label k = i; k < Un; k++) + { + scale += mag(U_[k][i]); + } + + if (scale != 0) + { + for (label k = i; k < Un; k++) + { + U_[k][i] /= scale; + s += U_[k][i]*U_[k][i]; + } + + scalar f = U_[i][i]; + g = -sign(Foam::sqrt(s), f); + scalar h = f*g - s; + U_[i][i] = f - g; + + for (label j = l-1; j < Um; j++) + { + s = 0; + for (label k = i; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + f = s/h; + for (label k = i; k < A.n(); k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label k = i; k < Un; k++) + { + U_[k][i] *= scale; + } + } + } + + S_[i] = scale*g; + + g = s = scale = 0; + + if (i+1 <= Un && i != Um) + { + for (label k = l-1; k < Um; k++) + { + scale += mag(U_[i][k]); + } + + if (scale != 0) + { + for (label k=l-1; k < Um; k++) + { + U_[i][k] /= scale; + s += U_[i][k]*U_[i][k]; + } + + scalar f = U_[i][l-1]; + g = -sign(Foam::sqrt(s),f); + scalar h = f*g - s; + U_[i][l-1] = f - g; + + for (label k = l-1; k < Um; k++) + { + rv1[k] = U_[i][k]/h; + } + + for (label j = l-1; j < Un; j++) + { + s = 0; + for (label k = l-1; k < Um; k++) + { + s += U_[j][k]*U_[i][k]; + } + + for (label k = l-1; k < Um; k++) + { + U_[j][k] += s*rv1[k]; + } + } + for (label k = l-1; k < Um; k++) + { + U_[i][k] *= scale; + } + } + } + + anorm = max(anorm, mag(S_[i]) + mag(rv1[i])); + } + + for (label i = Um-1; i >= 0; i--) + { + if (i < Um-1) + { + if (g != 0) + { + for (label j = l; j < Um; j++) + { + V_[j][i] = (U_[i][j]/U_[i][l])/g; + } + + for (label j=l; j < Um; j++) + { + s = 0; + for (label k = l; k < Um; k++) + { + s += U_[i][k]*V_[k][j]; + } + + for (label k = l; k < Um; k++) + { + V_[k][j] += s*V_[k][i]; + } + } + } + + for (label j = l; j < Um;j++) + { + V_[i][j] = V_[j][i] = 0.0; + } + } + + V_[i][i] = 1; + g = rv1[i]; + l = i; + } + + for (label i = min(Um, Un) - 1; i >= 0; i--) + { + l = i+1; + g = S_[i]; + + for (label j = l; j < Um; j++) + { + U_[i][j] = 0.0; + } + + if (g != 0) + { + g = 1.0/g; + + for (label j = l; j < Um; j++) + { + s = 0; + for (label k = l; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + scalar f = (s/U_[i][i])*g; + + for (label k = i; k < Un; k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label j = i; j < Un; j++) + { + U_[j][i] *= g; + } + } + else + { + for (label j = i; j < Un; j++) + { + U_[j][i] = 0.0; + } + } + + ++U_[i][i]; + } + + for (label k = Um-1; k >= 0; k--) + { + for (label its = 0; its < 35; its++) + { + bool flag = true; + + label nm; + for (l = k; l >= 0; l--) + { + nm = l-1; + if (mag(rv1[l]) + anorm == anorm) + { + flag = false; + break; + } + if (mag(S_[nm]) + anorm == anorm) break; + } + + if (flag) + { + scalar c = 0.0; + s = 1.0; + for (label i = l-1; i < k+1; i++) + { + scalar f = s*rv1[i]; + rv1[i] = c*rv1[i]; + + if (mag(f) + anorm == anorm) break; + + g = S_[i]; + scalar h = sqrtSumSqr(f, g); + S_[i] = h; + h = 1.0/h; + c = g*h; + s = -f*h; + + for (label j = 0; j < Un; j++) + { + scalar y = U_[j][nm]; + scalar z = U_[j][i]; + U_[j][nm] = y*c + z*s; + U_[j][i] = z*c - y*s; + } + } + } + + scalar z = S_[k]; + + if (l == k) + { + if (z < 0.0) + { + S_[k] = -z; + + for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k]; + } + break; + } + if (its == 34) + { + WarningIn("SVD::SVD(Matrix& A)") + << "no convergence in 35 SVD iterations" + << endl; + } + + scalar x = S_[l]; + nm = k-1; + scalar y = S_[nm]; + g = rv1[nm]; + scalar h = rv1[k]; + scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y); + g = sqrtSumSqr(f, 1.0); + f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x; + scalar c = 1.0; + s = 1.0; + + for (label j = l; j <= nm; j++) + { + label i = j + 1; + g = rv1[i]; + y = S_[i]; + h = s*g; + g = c*g; + scalar z = sqrtSumSqr(f, h); + rv1[j] = z; + c = f/z; + s = h/z; + f = x*c + g*s; + g = g*c - x*s; + h = y*s; + y *= c; + + for (label jj = 0; jj < Um; jj++) + { + x = V_[jj][j]; + z = V_[jj][i]; + V_[jj][j] = x*c + z*s; + V_[jj][i] = z*c - x*s; + } + + z = sqrtSumSqr(f, h); + S_[j] = z; + if (z) + { + z = 1.0/z; + c = f*z; + s = h*z; + } + f = c*g + s*y; + x = c*y - s*g; + + for (label jj=0; jj < Un; jj++) + { + y = U_[jj][j]; + z = U_[jj][i]; + U_[jj][j] = y*c + z*s; + U_[jj][i] = z*c - y*s; + } + } + rv1[l] = 0.0; + rv1[k] = f; + S_[k] = x; + } + } + + // zero singular values that are less than minCondition*maxS + const scalar minS = minCondition*S_[findMax(S_)]; + for (label i = 0; i < S_.size(); i++) + { + if (S_[i] <= minS) + { + //Info << "Removing " << S_[i] << " < " << minS << endl; + S_[i] = 0; + nZeros_++; + } + } + + // now multiply out to find the pseudo inverse of A, VSinvUt_ + multiply(VSinvUt_, V_, inv(S_), U_.T()); + + // test SVD + /*Matrix SVDA(A.n(), A.m()); + multiply(SVDA, U_, S_, transpose(V_)); + scalar maxDiff = 0; + scalar diff = 0; + for(label i = 0; i < A.n(); i++) + { + for(label j = 0; j < A.m(); j++) + { + diff = mag(A[i][j] - SVDA[i][j]); + if (diff > maxDiff) maxDiff = diff; + } + } + Info << "Maximum discrepancy between A and svd(A) = " << maxDiff << endl; + + if (maxDiff > 4) + { + Info << "singular values " << S_ << endl; + } + */ +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H new file mode 100644 index 0000000000..b597027732 --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + SVD + +Description + Singular value decomposition of a rectangular matrix. + +SourceFiles + SVD.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SVD_H +#define SVD_H + +#include "DiagonalMatrix.H" +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +/*---------------------------------------------------------------------------*\ + Class SVD Declaration +\*---------------------------------------------------------------------------*/ + +class SVD +{ + // Private data + + //- Rectangular matrix with the same dimensions as the input + Matrix U_; + + //- square matrix V + Matrix V_; + + //- The singular values + DiagonalMatrix S_; + + //- The matrix product V S^(-1) U^T + Matrix VSinvUt_; + + //- The number of zero singular values + label nZeros_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + SVD(const SVD&); + + //- Disallow default bitwise assignment + void operator=(const SVD&); + + template + inline const T sign(const T& a, const T& b); + + +public: + + // Constructors + + //- Construct from a rectangular Matrix + SVD(const Matrix& A, const scalar minCondition = 0); + + + // Access functions + + //- Return U + inline const Matrix& U() const; + + //- Return the square matrix V + inline const Matrix& V() const; + + //- Return the singular values + inline const DiagonalMatrix& S() const; + + //- Return VSinvUt (the pseudo inverse) + inline const Matrix& VSinvUt() const; + + //- Return the number of zero singular values + inline label nZeros() const; + + //- Return the minimum non-zero singular value + inline scalar minNonZeroS() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "SVDI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H new file mode 100644 index 0000000000..337647edfb --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // + +template +inline const T Foam::SVD::sign(const T& a, const T& b) +{ + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::Matrix& Foam::SVD::U() const +{ + return U_; +} + +inline const Foam::Matrix& Foam::SVD::V() const +{ + return V_; +} + +inline const Foam::DiagonalMatrix& Foam::SVD::S() const +{ + return S_; +} + +inline const Foam::Matrix& Foam::SVD::VSinvUt() const +{ + return VSinvUt_; +} + +inline Foam::label Foam::SVD::nZeros() const +{ + return nZeros_; +} + +inline Foam::scalar Foam::SVD::minNonZeroS() const +{ + scalar minS = S_[0]; + for(label i = 1; i < S_.size(); i++) + { + scalar s = S_[i]; + if (s > VSMALL && s < minS) minS = s; + } + return minS; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C index 2ff9e39f18..b88e167b6f 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C +++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C @@ -25,6 +25,8 @@ License \*---------------------------------------------------------------------------*/ #include "scalarMatrix.H" +#include "SVD.H" + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -136,7 +138,7 @@ void Foam::scalarMatrix::LUDecompose { Swap(matrixj[k], matrixiMax[k]); } - + vv[iMax] = vv[j]; } @@ -158,4 +160,159 @@ void Foam::scalarMatrix::LUDecompose } +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +void Foam::multiply +( + Matrix& ans, // value changed in return + const Matrix& A, + const Matrix& B +) +{ + if (A.m() != B.n()) + { + FatalErrorIn + ( + "multiply(" + "Matrix& answer " + "const Matrix& A, " + "const Matrix& B)" + ) << "A and B must have identical inner dimensions but A.m = " + << A.m() << " and B.n = " << B.n() + << abort(FatalError); + } + + ans = Matrix(A.n(), B.m(), scalar(0)); + + for(register label i = 0; i < A.n(); i++) + { + for(register label j = 0; j < B.m(); j++) + { + for(register label l = 0; l < B.n(); l++) + { + ans[i][j] += A[i][l]*B[l][j]; + } + } + } +} + + +void Foam::multiply +( + Matrix& ans, // value changed in return + const Matrix& A, + const Matrix& B, + const Matrix& C +) +{ + if (A.m() != B.n()) + { + FatalErrorIn + ( + "multiply(" + "const Matrix& A, " + "const Matrix& B, " + "const Matrix& C, " + "Matrix& answer)" + ) << "A and B must have identical inner dimensions but A.m = " + << A.m() << " and B.n = " << B.n() + << abort(FatalError); + } + + if (B.m() != C.n()) + { + FatalErrorIn + ( + "multiply(" + "const Matrix& A, " + "const Matrix& B, " + "const Matrix& C, " + "Matrix& answer)" + ) << "B and C must have identical inner dimensions but B.m = " + << B.m() << " and C.n = " << C.n() + << abort(FatalError); + } + + ans = Matrix(A.n(), C.m(), scalar(0)); + + for(register label i = 0; i < A.n(); i++) + { + for(register label g = 0; g < C.m(); g++) + { + for(register label l = 0; l < C.n(); l++) + { + scalar ab = 0; + for(register label j = 0; j < A.m(); j++) + { + ab += A[i][j]*B[j][l]; + } + ans[i][g] += C[l][g] * ab; + } + } + } +} + + +void Foam::multiply +( + Matrix& ans, // value changed in return + const Matrix& A, + const DiagonalMatrix& B, + const Matrix& C +) +{ + if (A.m() != B.size()) + { + FatalErrorIn + ( + "multiply(" + "const Matrix& A, " + "const DiagonalMatrix& B, " + "const Matrix& C, " + "Matrix& answer)" + ) << "A and B must have identical inner dimensions but A.m = " + << A.m() << " and B.n = " << B.size() + << abort(FatalError); + } + + if (B.size() != C.n()) + { + FatalErrorIn + ( + "multiply(" + "const Matrix& A, " + "const DiagonalMatrix& B, " + "const Matrix& C, " + "Matrix& answer)" + ) << "B and C must have identical inner dimensions but B.m = " + << B.size() << " and C.n = " << C.n() + << abort(FatalError); + } + + ans = Matrix(A.n(), C.m(), scalar(0)); + + for(register label i = 0; i < A.n(); i++) + { + for(register label g = 0; g < C.m(); g++) + { + for(register label l = 0; l < C.n(); l++) + { + ans[i][g] += C[l][g] * A[i][l]*B[l]; + } + } + } +} + + +Foam::Matrix Foam::SVDinv +( + const Matrix& A, + scalar minCondition +) +{ + SVD svd(A, minCondition); + return svd.VSinvUt(); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H index f2b0764f1a..2d22e71d4e 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H +++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H @@ -37,6 +37,7 @@ SourceFiles #define scalarMatrix_H #include "Matrix.H" +#include "DiagonalMatrix.H" #include "scalarField.H" #include "labelList.H" @@ -75,13 +76,13 @@ public: //- Solve the matrix using Gaussian elimination with pivoting, // returning the solution in the source - template - static void solve(Matrix& matrix, Field& source); + template + static void solve(Matrix& matrix, Field& source); //- Solve the matrix using Gaussian elimination with pivoting // and return the solution - template - void solve(Field& psi, const Field& source) const; + template + void solve(Field& psi, const Field& source) const; //- LU decompose the matrix with pivoting @@ -93,21 +94,50 @@ public: //- LU back-substitution with given source, returning the solution // in the source - template + template static void LUBacksubstitute ( const Matrix& luMmatrix, const labelList& pivotIndices, - Field& source + Field& source ); //- Solve the matrix using LU decomposition with pivoting // returning the LU form of the matrix and the solution in the source - template - static void LUsolve(Matrix& matrix, Field& source); + template + static void LUsolve(Matrix& matrix, Field& source); }; +// Global functions + +void multiply +( + Matrix& answer, // value changed in return + const Matrix& A, + const Matrix& B +); + +void multiply +( + Matrix& answer, // value changed in return + const Matrix& A, + const Matrix& B, + const Matrix& C +); + +void multiply +( + Matrix& answer, // value changed in return + const Matrix& A, + const DiagonalMatrix& B, + const Matrix& C +); + +//- Return the inverse of matrix A using SVD +Matrix SVDinv(const Matrix& A, scalar minCondition = 0); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C index b729e32fcc..fb48a9b26d 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C @@ -29,11 +29,11 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template +template void Foam::scalarMatrix::solve ( Matrix& tmpMatrix, - Field& sourceSol + Field& sourceSol ) { label n = tmpMatrix.n(); @@ -89,7 +89,7 @@ void Foam::scalarMatrix::solve // Back-substitution for (register label j=n-1; j>=0; j--) { - T ntempvec = pTraits::zero; + Type ntempvec = pTraits::zero; for (register label k=j+1; k -void Foam::scalarMatrix::solve(Field& psi, const Field& source) const +template +void Foam::scalarMatrix::solve(Field& psi, const Field& source) const { Matrix tmpMatrix = *this; psi = source; @@ -110,12 +110,12 @@ void Foam::scalarMatrix::solve(Field& psi, const Field& source) const } -template +template void Foam::scalarMatrix::LUBacksubstitute ( const Matrix& luMatrix, const labelList& pivotIndices, - Field& sourceSol + Field& sourceSol ) { label n = luMatrix.n(); @@ -125,7 +125,7 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=0; i::zero) + else if (sum != pTraits::zero) { ii = i+1; } @@ -146,11 +146,11 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=n-1; i>=0; i--) { - T sum = sourceSol[i]; + Type sum = sourceSol[i]; const scalar* __restrict__ luMatrixi = luMatrix[i]; for (register label j=i+1; j +template void Foam::scalarMatrix::LUsolve ( Matrix& matrix, - Field& sourceSol + Field& sourceSol ) { labelList pivotIndices(matrix.n()); diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C index 22de758f3a..1265b87ee5 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C @@ -28,19 +28,19 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -Foam::simpleMatrix::simpleMatrix(const label mSize) +template +Foam::simpleMatrix::simpleMatrix(const label mSize) : scalarMatrix(mSize), - source_(mSize, pTraits::zero) + source_(mSize, pTraits::zero) {} -template -Foam::simpleMatrix::simpleMatrix +template +Foam::simpleMatrix::simpleMatrix ( const scalarMatrix& matrix, - const Field& source + const Field& source ) : scalarMatrix(matrix), @@ -48,8 +48,8 @@ Foam::simpleMatrix::simpleMatrix {} -template -Foam::simpleMatrix::simpleMatrix(Istream& is) +template +Foam::simpleMatrix::simpleMatrix(Istream& is) : scalarMatrix(is), source_(is) @@ -58,11 +58,11 @@ Foam::simpleMatrix::simpleMatrix(Istream& is) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -Foam::Field Foam::simpleMatrix::solve() const +template +Foam::Field Foam::simpleMatrix::solve() const { scalarMatrix tmpMatrix = *this; - Field sourceSol = source_; + Field sourceSol = source_; scalarMatrix::solve(tmpMatrix, sourceSol); @@ -70,11 +70,11 @@ Foam::Field Foam::simpleMatrix::solve() const } -template -Foam::Field Foam::simpleMatrix::LUsolve() const +template +Foam::Field Foam::simpleMatrix::LUsolve() const { scalarMatrix luMatrix = *this; - Field sourceSol = source_; + Field sourceSol = source_; scalarMatrix::LUsolve(luMatrix, sourceSol); @@ -84,26 +84,26 @@ Foam::Field Foam::simpleMatrix::LUsolve() const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -void Foam::simpleMatrix::operator=(const simpleMatrix& m) +template +void Foam::simpleMatrix::operator=(const simpleMatrix& m) { if (this == &m) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Attempted assignment to self" << abort(FatalError); } if (n() != m.n()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size matrices" << abort(FatalError); } if (source_.size() != m.source_.size()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size source vectors" << abort(FatalError); } @@ -115,14 +115,14 @@ void Foam::simpleMatrix::operator=(const simpleMatrix& m) // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // -template -Foam::simpleMatrix Foam::operator+ +template +Foam::simpleMatrix Foam::operator+ ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( static_cast(m1) + static_cast(m2), @@ -131,14 +131,14 @@ Foam::simpleMatrix Foam::operator+ } -template -Foam::simpleMatrix Foam::operator- +template +Foam::simpleMatrix Foam::operator- ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( static_cast(m1) - static_cast(m2), @@ -147,17 +147,17 @@ Foam::simpleMatrix Foam::operator- } -template -Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) +template +Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) { - return simpleMatrix(s*m.matrix_, s*m.source_); + return simpleMatrix(s*m.matrix_, s*m.source_); } // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -template -Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) { os << static_cast(m) << nl << m.source_; return os; diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H index d1c411452b..d7b1a027be 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H @@ -45,35 +45,35 @@ namespace Foam // Forward declaration of friend functions and operators -template +template class simpleMatrix; -template -simpleMatrix operator+ +template +simpleMatrix operator+ ( - const simpleMatrix&, - const simpleMatrix& + const simpleMatrix&, + const simpleMatrix& ); -template -simpleMatrix operator- +template +simpleMatrix operator- ( - const simpleMatrix&, - const simpleMatrix& + const simpleMatrix&, + const simpleMatrix& ); -template -simpleMatrix operator* +template +simpleMatrix operator* ( const scalar, - const simpleMatrix& + const simpleMatrix& ); -template +template Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); @@ -81,14 +81,14 @@ Ostream& operator<< Class simpleMatrix Declaration \*---------------------------------------------------------------------------*/ -template +template class simpleMatrix : public scalarMatrix { // Private data - Field source_; + Field source_; public: @@ -99,25 +99,25 @@ public: simpleMatrix(const label); //- Construct from components - simpleMatrix(const scalarMatrix&, const Field&); + simpleMatrix(const scalarMatrix&, const Field&); //- Construct from Istream simpleMatrix(Istream&); //- Construct as copy - simpleMatrix(const simpleMatrix&); + simpleMatrix(const simpleMatrix&); // Member Functions // Access - Field& source() + Field& source() { return source_; } - const Field& source() const + const Field& source() const { return source_; } @@ -125,45 +125,45 @@ public: //- Solve the matrix using Gaussian elimination with pivoting // and return the solution - Field solve() const; + Field solve() const; //- Solve the matrix using LU decomposition with pivoting // and return the solution - Field LUsolve() const; + Field LUsolve() const; // Member Operators - void operator=(const simpleMatrix&); + void operator=(const simpleMatrix&); // Friend Operators - friend simpleMatrix operator+ + friend simpleMatrix operator+ ( - const simpleMatrix&, - const simpleMatrix& + const simpleMatrix&, + const simpleMatrix& ); - friend simpleMatrix operator- + friend simpleMatrix operator- ( - const simpleMatrix&, - const simpleMatrix& + const simpleMatrix&, + const simpleMatrix& ); - friend simpleMatrix operator* + friend simpleMatrix operator* ( const scalar, - const simpleMatrix& + const simpleMatrix& ); // Ostream Operator - friend Ostream& operator<< + friend Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); }; diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 0c3fae1887..25137c58c7 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -38,6 +38,9 @@ fvMeshMapper = fvMesh/fvMeshMapper $(fvMeshMapper)/fvPatchMapper.C $(fvMeshMapper)/fvSurfaceMapper.C +extendedStencil = fvMesh/extendedStencil +$(extendedStencil)/extendedStencil.C + fvPatchFields = fields/fvPatchFields $(fvPatchFields)/fvPatchField/fvPatchFields.C @@ -165,6 +168,8 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C +$(schemes)/quadraticFit/quadraticFit.C +$(schemes)/quadraticFit/quadraticFitData.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C new file mode 100644 index 0000000000..bb4c6c1534 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C @@ -0,0 +1,525 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 + +\*---------------------------------------------------------------------------*/ + +#include "extendedStencil.H" +#include "globalIndex.H" +#include "syncTools.H" +#include "SortableList.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Calculates per face a list of global cell/face indices. +void Foam::extendedStencil::calcFaceStencils +( + const polyMesh& mesh, + const globalIndex& globalNumbering +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + const labelList& own = mesh.faceOwner(); + const labelList& nei = mesh.faceNeighbour(); + + + // Determine neighbouring global cell or boundary face + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelList neiGlobal(nBnd); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + // For coupled faces get the cell on the other side + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]); + faceI++; + } + } + else if (isA(pp)) + { + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = -1; + faceI++; + } + } + else + { + // For noncoupled faces get the boundary face. + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = + globalNumbering.toGlobal(mesh.nCells()+bFaceI); + faceI++; + } + } + } + syncTools::swapBoundaryFaceList(mesh, neiGlobal, false); + + + // Determine cellCells in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList globalCellCells(mesh.nCells()); + forAll(globalCellCells, cellI) + { + const cell& cFaces = mesh.cells()[cellI]; + + labelList& cCells = globalCellCells[cellI]; + + cCells.setSize(cFaces.size()); + + // Collect neighbouring cells/faces + label nNbr = 0; + forAll(cFaces, i) + { + label faceI = cFaces[i]; + + if (mesh.isInternalFace(faceI)) + { + label nbrCellI = own[faceI]; + if (nbrCellI == cellI) + { + nbrCellI = nei[faceI]; + } + cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI); + } + else + { + label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()]; + if (nbrCellI != -1) + { + cCells[nNbr++] = nbrCellI; + } + } + } + cCells.setSize(nNbr); + } + + + // Determine neighbouring global cell Cells + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList neiGlobalCellCells(nBnd); + for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + { + neiGlobalCellCells[faceI-mesh.nInternalFaces()] = + globalCellCells[own[faceI]]; + } + syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false); + + + + // Construct stencil in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + stencil_.setSize(mesh.nFaces()); + + labelHashSet faceStencil; + + for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + label globalNei = globalNumbering.toGlobal(nei[faceI]); + faceStencil.insert(globalNei); + const labelList& neiCCells = globalCellCells[nei[faceI]]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + } + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + // Get the coupled cell + label globalNei = neiGlobal[faceI-mesh.nInternalFaces()]; + faceStencil.insert(globalNei); + // And the neighbours of the coupled cell + const labelList& neiCCells = + neiGlobalCellCells[faceI-mesh.nInternalFaces()]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "coupledface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + else if (!isA(pp)) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "boundaryface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + } +} + + +// Calculates extended stencil. This is per face +// - owner +// - cellCells of owner +// - neighbour +// - cellCells of neighbour +// It comes in two parts: +// - a map which collects/distributes all necessary data in a compact array +// - the stencil (a labelList per face) which is a set of indices into this +// compact array. +// The compact array is laid out as follows: +// - first data for current processor (Pstream::myProcNo()) +// - all cells +// - all boundary faces +// - then per processor +// - all used cells and boundary faces +void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) +{ + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + + // Global numbering for cells and boundary faces + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + globalIndex globalNumbering(mesh.nCells()+nBnd); + + + // Calculate stencil in global cell indices + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + calcFaceStencils(mesh, globalNumbering); + + + // Convert stencil to schedule + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // We now know what information we need from other processors. This needs + // to be converted into what information I need to send as well + // (mapDistribute) + + + // 1. Construct per processor compact addressing of the global cells + // needed. The ones from the local processor are not included since + // these are always all needed. + List > globalToProc(Pstream::nProcs()); + { + const labelList& procPatchMap = mesh.globalData().procPatchMap(); + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Presize with (as estimate) size of patch to neighbour. + forAll(procPatchMap, procI) + { + if (procPatchMap[procI] != -1) + { + globalToProc[procI].resize + ( + patches[procPatchMap[procI]].size() + ); + } + } + + // Collect all (non-local) globalcells/faces needed. + forAll(stencil_, faceI) + { + const labelList& stencilCells = stencil_[faceI]; + + forAll(stencilCells, i) + { + label globalCellI = stencilCells[i]; + label procI = globalNumbering.whichProcID(stencilCells[i]); + + if (procI != Pstream::myProcNo()) + { + label nCompact = globalToProc[procI].size(); + globalToProc[procI].insert(globalCellI, nCompact); + } + } + } + // Sort global cells needed (not really necessary) + forAll(globalToProc, procI) + { + if (procI != Pstream::myProcNo()) + { + Map