From 26f2c24d84d653169811d752b5ff3884c2c9c87a Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
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 | 89 +--
.../{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, 2725 insertions(+), 252 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 7b137dedc03..57dc7f8dc70 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 00000000000..c023e59c5a2
--- /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<class Type>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Type>& a)
+:
+ List<Type>(min(a.n(), a.m()))
+{
+ forAll(*this, i)
+ {
+ this->operator[](i) = a[i][i];
+ }
+}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size)
+:
+ List<Type>(size)
+{}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size, const Type& val)
+:
+ List<Type>(size, val)
+{}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::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<class Type>
+Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& A)
+{
+ DiagonalMatrix<Type> 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 00000000000..36138fc25de
--- /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<Type>
+
+Description
+ DiagonalMatrix<Type> 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 Type> class Matrix;
+
+/*---------------------------------------------------------------------------*\
+ Class DiagonalMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class DiagonalMatrix
+:
+ public List<Type>
+{
+public:
+
+ // Constructors
+
+ //- Construct from diagonal component of a Matrix
+ DiagonalMatrix<Type>(const Matrix<Type>&);
+
+ //- Construct empty from size
+ DiagonalMatrix<Type>(const label size);
+
+ //- Construct from size and a value
+ DiagonalMatrix<Type>(const label, const Type&);
+
+
+ // Member functions
+
+ //- Invert the diaganol matrix and return itself
+ DiagonalMatrix<Type>& invert();
+};
+
+
+// Global functions
+
+//- Return the diagonal Matrix inverse
+template<class Type>
+DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>&);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // 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 cbb59e2c864..ac0bc17207e 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
@@ -28,16 +28,16 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template<class T>
-void Foam::LUscalarMatrix::solve(Field<T>& sourceSol) const
+template<class Type>
+void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
{
if (Pstream::parRun())
{
- Field<T> completeSourceSol(n());
+ Field<Type> completeSourceSol(n());
if (Pstream::master())
{
- typename Field<T>::subField
+ typename Field<Type>::subField
(
completeSourceSol,
sourceSol.size()
@@ -58,7 +58,7 @@ void Foam::LUscalarMatrix::solve(Field<T>& 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<T>& sourceSol) const
{
LUBacksubstitute(*this, pivotIndices_, completeSourceSol);
- sourceSol = typename Field<T>::subField
+ sourceSol = typename Field<Type>::subField
(
completeSourceSol,
sourceSol.size()
@@ -98,7 +98,7 @@ void Foam::LUscalarMatrix::solve(Field<T>& 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 733e04034ef..4fbef71c2b4 100644
--- a/src/OpenFOAM/containers/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -28,13 +28,13 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-template<class T>
-void Foam::Matrix<T>::allocate()
+template<class Type>
+void Foam::Matrix<Type>::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<n_; i++)
{
@@ -46,12 +46,12 @@ void Foam::Matrix<T>::allocate()
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
-template<class T>
-Foam::Matrix<T>::~Matrix()
+template<class Type>
+Foam::Matrix<Type>::~Matrix()
{
if (v_)
{
- delete[] (v_[0]);
+ delete[] (v_[0]);
delete[] v_;
}
}
@@ -59,16 +59,16 @@ Foam::Matrix<T>::~Matrix()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template<class T>
-const Foam::Matrix<T>& Foam::Matrix<T>::null()
+template<class Type>
+const Foam::Matrix<Type>& Foam::Matrix<Type>::null()
{
- Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
+ Matrix<Type>* nullPtr = reinterpret_cast<Matrix<Type>*>(NULL);
return *nullPtr;
}
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m)
+template<class Type>
+Foam::Matrix<Type>::Matrix(const label n, const label m)
:
n_(n),
m_(m),
@@ -76,7 +76,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
{
if (n_ < 0 || m_ < 0)
{
- FatalErrorIn("Matrix<T>::Matrix(const label n, const label m)")
+ FatalErrorIn("Matrix<Type>::Matrix(const label n, const label m)")
<< "bad n, m " << n_ << ", " << m_
<< abort(FatalError);
}
@@ -85,8 +85,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
}
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
+template<class Type>
+Foam::Matrix<Type>::Matrix(const label n, const label m, const Type& a)
:
n_(n),
m_(m),
@@ -96,7 +96,7 @@ 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<Type>::Matrix(const label n, const label m, const T&)"
) << "bad n, m " << n_ << ", " << m_
<< abort(FatalError);
}
@@ -105,7 +105,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 +117,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
}
-template<class T>
-Foam::Matrix<T>::Matrix(const Matrix<T>& a)
+template<class Type>
+Foam::Matrix<Type>::Matrix(const Matrix<Type>& a)
:
n_(a.n_),
m_(a.m_),
@@ -127,8 +127,8 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& 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<nm; i++)
@@ -138,13 +138,13 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
}
}
-
-template<class T>
-void Foam::Matrix<T>::clear()
+
+template<class Type>
+void Foam::Matrix<Type>::clear()
{
if (v_)
{
- delete[] (v_[0]);
+ delete[] (v_[0]);
delete[] v_;
}
n_ = 0;
@@ -153,8 +153,8 @@ void Foam::Matrix<T>::clear()
}
-template<class T>
-void Foam::Matrix<T>::transfer(Matrix<T>& a)
+template<class Type>
+void Foam::Matrix<Type>::transfer(Matrix<Type>& a)
{
clear();
@@ -171,12 +171,30 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-template<class T>
-void Foam::Matrix<T>::operator=(const T& t)
+template<class Type>
+Foam::Matrix<Type> Foam::Matrix<Type>::T() const
+{
+ const Matrix<Type>& A = *this;
+ Matrix<Type> 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;
+}
+
+
+template<class Type>
+void Foam::Matrix<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 +206,12 @@ void Foam::Matrix<T>::operator=(const T& t)
// Assignment operator. Takes linear time.
-template<class T>
-void Foam::Matrix<T>::operator=(const Matrix<T>& a)
+template<class Type>
+void Foam::Matrix<Type>::operator=(const Matrix<Type>& a)
{
if (this == &a)
{
- FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
+ FatalErrorIn("Matrix<Type>::operator=(const Matrix<Type>&)")
<< "attempted assignment to self"
<< abort(FatalError);
}
@@ -204,12 +222,12 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& 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<nm; i++)
@@ -222,15 +240,15 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
-template<class T>
-const T& Foam::max(const Matrix<T>& a)
+template<class Type>
+const Type& Foam::max(const Matrix<Type>& 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<nm; i++)
{
@@ -244,7 +262,7 @@ const T& Foam::max(const Matrix<T>& a)
}
else
{
- FatalErrorIn("max(const Matrix<T>&)")
+ FatalErrorIn("max(const Matrix<Type>&)")
<< "matrix is empty"
<< abort(FatalError);
@@ -254,15 +272,15 @@ const T& Foam::max(const Matrix<T>& a)
}
-template<class T>
-const T& Foam::min(const Matrix<T>& a)
+template<class Type>
+const Type& Foam::min(const Matrix<Type>& 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<nm; i++)
{
@@ -276,7 +294,7 @@ const T& Foam::min(const Matrix<T>& a)
}
else
{
- FatalErrorIn("min(const Matrix<T>&)")
+ FatalErrorIn("min(const Matrix<Type>&)")
<< "matrix is empty"
<< abort(FatalError);
@@ -288,15 +306,15 @@ const T& Foam::min(const Matrix<T>& a)
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
-template<class T>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
+template<class Type>
+Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a)
{
- Matrix<T> na(a.n_, a.m_);
+ Matrix<Type> 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<nm; i++)
@@ -309,30 +327,34 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
}
-template<class T>
-Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
+template<class Type>
+Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<Type>& b)
{
if (a.n_ != b.n_)
{
- FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
- << "attempted add matrices with different number of rows: "
+ FatalErrorIn
+ (
+ "Matrix<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
+ ) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
- FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
- << "attempted add matrices with different number of columns: "
+ FatalErrorIn
+ (
+ "Matrix<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
+ ) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_
<< abort(FatalError);
}
- Matrix<T> ab(a.n_, a.m_);
+ Matrix<Type> 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<nm; i++)
@@ -344,30 +366,34 @@ Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
}
-template<class T>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
+template<class Type>
+Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<Type>& b)
{
if (a.n_ != b.n_)
{
- FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
- << "attempted add matrices with different number of rows: "
+ FatalErrorIn
+ (
+ "Matrix<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
+ ) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
- FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
- << "attempted add matrices with different number of columns: "
+ FatalErrorIn
+ (
+ "Matrix<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
+ ) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_
<< abort(FatalError);
}
- Matrix<T> ab(a.n_, a.m_);
+ Matrix<Type> 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<nm; i++)
@@ -379,15 +405,15 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
}
-template<class T>
-Foam::Matrix<T> Foam::operator*(const scalar s, const Matrix<T>& a)
+template<class Type>
+Foam::Matrix<Type> Foam::operator*(const scalar s, const Matrix<Type>& a)
{
- Matrix<T> sa(a.n_, a.m_);
+ Matrix<Type> 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<nm; i++)
diff --git a/src/OpenFOAM/containers/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
similarity index 70%
rename from src/OpenFOAM/containers/Matrix/Matrix.H
rename to src/OpenFOAM/matrices/Matrix/Matrix.H
index ab68ce5badb..6b4da39d4f2 100644
--- a/src/OpenFOAM/containers/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -51,25 +51,37 @@ namespace Foam
// Forward declaration of friend functions and operators
-template<class T> class Matrix;
-
-template<class T> const T& max(const Matrix<T>&);
-template<class T> const T& min(const Matrix<T>&);
-
-template<class T> Matrix<T> operator-(const Matrix<T>&);
-template<class T> Matrix<T> operator+(const Matrix<T>&, const Matrix<T>&);
-template<class T> Matrix<T> operator-(const Matrix<T>&, const Matrix<T>&);
-template<class T> Matrix<T> operator*(const scalar, const Matrix<T>&);
-
-template<class T> Istream& operator>>(Istream&, Matrix<T>&);
-template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
+template<class Type> class Matrix;
+
+template<class Type> const Type& max(const Matrix<Type>&);
+template<class Type> const Type& min(const Matrix<Type>&);
+
+template<class Type> Matrix<Type> operator-(const Matrix<Type>&);
+template<class Type> Matrix<Type> operator+
+(
+ const Matrix<Type>&,
+ const Matrix<Type>&
+);
+template<class Type> Matrix<Type> operator-
+(
+ const Matrix<Type>&,
+ const Matrix<Type>&
+);
+template<class Type> Matrix<Type> operator*
+(
+ const scalar,
+ const Matrix<Type>&
+);
+
+template<class Type> Istream& operator>>(Istream&, Matrix<Type>&);
+template<class Type> Ostream& operator<<(Ostream&, const Matrix<Type>&);
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
\*---------------------------------------------------------------------------*/
-template<class T>
+template<class Type>
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<T>&);
+ Matrix(const Matrix<Type>&);
//- Construct from Istream.
Matrix(Istream&);
//- Clone
- inline autoPtr<Matrix<T> > clone() const;
+ inline autoPtr<Matrix<Type> > clone() const;
// Destructor
@@ -117,7 +129,7 @@ public:
// Member functions
//- Return a null Matrix
- static const Matrix<T>& null();
+ static const Matrix<Type>& 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<T>&);
+ void transfer(Matrix<Type>&);
// 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<Type> T() const;
//- Assignment operator. Takes linear time.
- void operator=(const Matrix<T>&);
+ void operator=(const Matrix<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 const Type& max<Type>(const Matrix<Type>&);
+ friend const Type& min<Type>(const Matrix<Type>&);
// 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>&);
+ friend Matrix<Type> operator-<Type>(const Matrix<Type>&);
+ friend Matrix<Type> operator+<Type>
+ (
+ const Matrix<Type>&,
+ const Matrix<Type>&
+ );
+ friend Matrix<Type> operator-<Type>
+ (
+ const Matrix<Type>&,
+ const Matrix<Type>&
+ );
+ friend Matrix<Type> operator*<Type>
+ (
+ const scalar,
+ const Matrix<Type>&
+ );
// IOstream operators
//- Read Matrix from Istream, discarding contents of existing Matrix.
- friend Istream& operator>> <T>(Istream&, Matrix<T>&);
+ friend Istream& operator>> <Type>(Istream&, Matrix<Type>&);
// Write Matrix to Ostream.
- friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
+ friend Ostream& operator<< <Type>(Ostream&, const Matrix<Type>&);
};
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 d0e3fe8fdf3..8323174d4f9 100644
--- a/src/OpenFOAM/containers/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -26,8 +26,8 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-template<class T>
-inline Foam::Matrix<T>::Matrix()
+template<class Type>
+inline Foam::Matrix<Type>::Matrix()
:
n_(0),
m_(0),
@@ -35,71 +35,67 @@ inline Foam::Matrix<T>::Matrix()
{}
-template<class T>
-inline Foam::autoPtr<Foam::Matrix<T> > Foam::Matrix<T>::clone() const
+template<class Type>
+inline Foam::autoPtr<Foam::Matrix<Type> > Foam::Matrix<Type>::clone() const
{
- return autoPtr<Matrix<T> >(new Matrix<T>(*this));
+ return autoPtr<Matrix<Type> >(new Matrix<Type>(*this));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return the number of rows
-template<class T>
-inline Foam::label Foam::Matrix<T>::n() const
+template<class Type>
+inline Foam::label Foam::Matrix<Type>::n() const
{
return n_;
}
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::m() const
+template<class Type>
+inline Foam::label Foam::Matrix<Type>::m() const
{
return m_;
}
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::size() const
+template<class Type>
+inline Foam::label Foam::Matrix<Type>::size() const
{
return n_*m_;
}
-// Check index i is within valid range (0 ... n-1).
-template<class T>
-inline void Foam::Matrix<T>::checki(const label i) const
+template<class Type>
+inline void Foam::Matrix<Type>::checki(const label i) const
{
if (!n_)
{
- FatalErrorIn("Matrix<T>::checki(const label)")
+ FatalErrorIn("Matrix<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<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 T>
-inline void Foam::Matrix<T>::checkj(const label j) const
+template<class Type>
+inline void Foam::Matrix<Type>::checkj(const label j) const
{
if (!m_)
{
- FatalErrorIn("Matrix<T>::checkj(const label)")
+ FatalErrorIn("Matrix<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<Type>::checkj(const label)")
<< "index " << j << " out of range 0 ... " << m_-1
<< abort(FatalError);
}
@@ -108,9 +104,8 @@ inline void Foam::Matrix<T>::checkj(const label j) const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-// Return subscript-checked element access
-template<class T>
-inline T* Foam::Matrix<T>::operator[](const label i)
+template<class Type>
+inline Type* Foam::Matrix<Type>::operator[](const label i)
{
# ifdef FULLDEBUG
checki(i);
@@ -119,9 +114,8 @@ inline T* Foam::Matrix<T>::operator[](const label i)
}
-// Return subscript-checked const element access
-template<class T>
-inline const T* Foam::Matrix<T>::operator[](const label i) const
+template<class Type>
+inline const Type* Foam::Matrix<Type>::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 710e42876e6..9b5a922f0d5 100644
--- a/src/OpenFOAM/containers/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -32,8 +32,8 @@ License
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
-template<class T>
-Foam::Matrix<T>::Matrix(Istream& is)
+template<class Type>
+Foam::Matrix<Type>::Matrix(Istream& is)
:
n_(0),
m_(0),
@@ -43,17 +43,17 @@ Foam::Matrix<T>::Matrix(Istream& is)
}
-template<class T>
-Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
+template<class Type>
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
{
// Anull matrix
M.clear();
- is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
+ is.fatalCheck("operator>>(Istream&, Matrix<Type>&)");
token firstToken(is);
- is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
+ is.fatalCheck("operator>>(Istream&, Matrix<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<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<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<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<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>
-Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
+template<class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<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;
@@ -187,7 +187,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
// Write end of contents delimiter
os << token::END_BLOCK;
}
- else if (nm < 10 && contiguous<T>())
+ else if (nm < 10 && contiguous<Type>())
{
// Write size of list and start contents delimiter
os << token::BEGIN_LIST;
@@ -246,7 +246,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));
}
}
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
new file mode 100644
index 00000000000..a1355c393d1
--- /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<scalar>& 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<scalar>& 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<scalar> 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 00000000000..b597027732d
--- /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<scalar> U_;
+
+ //- square matrix V
+ Matrix<scalar> V_;
+
+ //- The singular values
+ DiagonalMatrix<scalar> S_;
+
+ //- The matrix product V S^(-1) U^T
+ Matrix<scalar> 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<class T>
+ inline const T sign(const T& a, const T& b);
+
+
+public:
+
+ // Constructors
+
+ //- Construct from a rectangular Matrix
+ SVD(const Matrix<scalar>& A, const scalar minCondition = 0);
+
+
+ // Access functions
+
+ //- Return U
+ inline const Matrix<scalar>& U() const;
+
+ //- Return the square matrix V
+ inline const Matrix<scalar>& V() const;
+
+ //- Return the singular values
+ inline const DiagonalMatrix<scalar>& S() const;
+
+ //- Return VSinvUt (the pseudo inverse)
+ inline const Matrix<scalar>& 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 00000000000..337647edfb5
--- /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<class T>
+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::scalar>& Foam::SVD::U() const
+{
+ return U_;
+}
+
+inline const Foam::Matrix<Foam::scalar>& Foam::SVD::V() const
+{
+ return V_;
+}
+
+inline const Foam::DiagonalMatrix<Foam::scalar>& Foam::SVD::S() const
+{
+ return S_;
+}
+
+inline const Foam::Matrix<Foam::scalar>& 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 2ff9e39f18a..b88e167b6fb 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<scalar>& ans, // value changed in return
+ const Matrix<scalar>& A,
+ const Matrix<scalar>& B
+)
+{
+ if (A.m() != B.n())
+ {
+ FatalErrorIn
+ (
+ "multiply("
+ "Matrix<scalar>& answer "
+ "const Matrix<scalar>& A, "
+ "const Matrix<scalar>& B)"
+ ) << "A and B must have identical inner dimensions but A.m = "
+ << A.m() << " and B.n = " << B.n()
+ << abort(FatalError);
+ }
+
+ ans = Matrix<scalar>(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<scalar>& ans, // value changed in return
+ const Matrix<scalar>& A,
+ const Matrix<scalar>& B,
+ const Matrix<scalar>& C
+)
+{
+ if (A.m() != B.n())
+ {
+ FatalErrorIn
+ (
+ "multiply("
+ "const Matrix<scalar>& A, "
+ "const Matrix<scalar>& B, "
+ "const Matrix<scalar>& C, "
+ "Matrix<scalar>& 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<scalar>& A, "
+ "const Matrix<scalar>& B, "
+ "const Matrix<scalar>& C, "
+ "Matrix<scalar>& answer)"
+ ) << "B and C must have identical inner dimensions but B.m = "
+ << B.m() << " and C.n = " << C.n()
+ << abort(FatalError);
+ }
+
+ ans = Matrix<scalar>(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<scalar>& ans, // value changed in return
+ const Matrix<scalar>& A,
+ const DiagonalMatrix<scalar>& B,
+ const Matrix<scalar>& C
+)
+{
+ if (A.m() != B.size())
+ {
+ FatalErrorIn
+ (
+ "multiply("
+ "const Matrix<scalar>& A, "
+ "const DiagonalMatrix<scalar>& B, "
+ "const Matrix<scalar>& C, "
+ "Matrix<scalar>& 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<scalar>& A, "
+ "const DiagonalMatrix<scalar>& B, "
+ "const Matrix<scalar>& C, "
+ "Matrix<scalar>& answer)"
+ ) << "B and C must have identical inner dimensions but B.m = "
+ << B.size() << " and C.n = " << C.n()
+ << abort(FatalError);
+ }
+
+ ans = Matrix<scalar>(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::scalar> Foam::SVDinv
+(
+ const Matrix<scalar>& 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 f2b0764f1a5..2d22e71d4ea 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<class T>
- static void solve(Matrix<scalar>& matrix, Field<T>& source);
+ template<class Type>
+ static void solve(Matrix<scalar>& matrix, Field<Type>& source);
//- Solve the matrix using Gaussian elimination with pivoting
// and return the solution
- template<class T>
- void solve(Field<T>& psi, const Field<T>& source) const;
+ template<class Type>
+ void solve(Field<Type>& psi, const Field<Type>& 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<class T>
+ template<class Type>
static void LUBacksubstitute
(
const Matrix<scalar>& luMmatrix,
const labelList& pivotIndices,
- Field<T>& source
+ Field<Type>& source
);
//- Solve the matrix using LU decomposition with pivoting
// returning the LU form of the matrix and the solution in the source
- template<class T>
- static void LUsolve(Matrix<scalar>& matrix, Field<T>& source);
+ template<class Type>
+ static void LUsolve(Matrix<scalar>& matrix, Field<Type>& source);
};
+// Global functions
+
+void multiply
+(
+ Matrix<scalar>& answer, // value changed in return
+ const Matrix<scalar>& A,
+ const Matrix<scalar>& B
+);
+
+void multiply
+(
+ Matrix<scalar>& answer, // value changed in return
+ const Matrix<scalar>& A,
+ const Matrix<scalar>& B,
+ const Matrix<scalar>& C
+);
+
+void multiply
+(
+ Matrix<scalar>& answer, // value changed in return
+ const Matrix<scalar>& A,
+ const DiagonalMatrix<scalar>& B,
+ const Matrix<scalar>& C
+);
+
+//- Return the inverse of matrix A using SVD
+Matrix<scalar> SVDinv(const Matrix<scalar>& A, scalar minCondition = 0);
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
index b729e32fcc7..fb48a9b26d5 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
@@ -29,11 +29,11 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template<class T>
+template<class Type>
void Foam::scalarMatrix::solve
(
Matrix<scalar>& tmpMatrix,
- Field<T>& sourceSol
+ Field<Type>& 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<T>::zero;
+ Type ntempvec = pTraits<Type>::zero;
for (register label k=j+1; k<n; k++)
{
@@ -101,8 +101,8 @@ void Foam::scalarMatrix::solve
}
-template<class T>
-void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
+template<class Type>
+void Foam::scalarMatrix::solve(Field<Type>& psi, const Field<Type>& source) const
{
Matrix<scalar> tmpMatrix = *this;
psi = source;
@@ -110,12 +110,12 @@ void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
}
-template<class T>
+template<class Type>
void Foam::scalarMatrix::LUBacksubstitute
(
const Matrix<scalar>& luMatrix,
const labelList& pivotIndices,
- Field<T>& sourceSol
+ Field<Type>& sourceSol
)
{
label n = luMatrix.n();
@@ -125,7 +125,7 @@ void Foam::scalarMatrix::LUBacksubstitute
for (register label i=0; i<n; i++)
{
label ip = pivotIndices[i];
- T sum = sourceSol[ip];
+ Type sum = sourceSol[ip];
sourceSol[ip] = sourceSol[i];
const scalar* __restrict__ luMatrixi = luMatrix[i];
@@ -136,7 +136,7 @@ void Foam::scalarMatrix::LUBacksubstitute
sum -= luMatrixi[j]*sourceSol[j];
}
}
- else if (sum != pTraits<T>::zero)
+ else if (sum != pTraits<Type>::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<n; j++)
- {
+ {
sum -= luMatrixi[j]*sourceSol[j];
}
@@ -159,11 +159,11 @@ void Foam::scalarMatrix::LUBacksubstitute
}
-template<class T>
+template<class Type>
void Foam::scalarMatrix::LUsolve
(
Matrix<scalar>& matrix,
- Field<T>& sourceSol
+ Field<Type>& sourceSol
)
{
labelList pivotIndices(matrix.n());
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 22de758f3ac..1265b87ee57 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -28,19 +28,19 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix(const label mSize)
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
:
scalarMatrix(mSize),
- source_(mSize, pTraits<T>::zero)
+ source_(mSize, pTraits<Type>::zero)
{}
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix
(
const scalarMatrix& matrix,
- const Field<T>& source
+ const Field<Type>& source
)
:
scalarMatrix(matrix),
@@ -48,8 +48,8 @@ Foam::simpleMatrix<T>::simpleMatrix
{}
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix(Istream& is)
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
:
scalarMatrix(is),
source_(is)
@@ -58,11 +58,11 @@ Foam::simpleMatrix<T>::simpleMatrix(Istream& is)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template<class T>
-Foam::Field<T> Foam::simpleMatrix<T>::solve() const
+template<class Type>
+Foam::Field<Type> Foam::simpleMatrix<Type>::solve() const
{
scalarMatrix tmpMatrix = *this;
- Field<T> sourceSol = source_;
+ Field<Type> sourceSol = source_;
scalarMatrix::solve(tmpMatrix, sourceSol);
@@ -70,11 +70,11 @@ Foam::Field<T> Foam::simpleMatrix<T>::solve() const
}
-template<class T>
-Foam::Field<T> Foam::simpleMatrix<T>::LUsolve() const
+template<class Type>
+Foam::Field<Type> Foam::simpleMatrix<Type>::LUsolve() const
{
scalarMatrix luMatrix = *this;
- Field<T> sourceSol = source_;
+ Field<Type> sourceSol = source_;
scalarMatrix::LUsolve(luMatrix, sourceSol);
@@ -84,26 +84,26 @@ Foam::Field<T> Foam::simpleMatrix<T>::LUsolve() const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-template<class T>
-void Foam::simpleMatrix<T>::operator=(const simpleMatrix<T>& m)
+template<class Type>
+void Foam::simpleMatrix<Type>::operator=(const simpleMatrix<Type>& m)
{
if (this == &m)
{
- FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+ FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
<< "Attempted assignment to self"
<< abort(FatalError);
}
if (n() != m.n())
{
- FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+ FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
<< "Different size matrices"
<< abort(FatalError);
}
if (source_.size() != m.source_.size())
{
- FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+ FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
<< "Different size source vectors"
<< abort(FatalError);
}
@@ -115,14 +115,14 @@ void Foam::simpleMatrix<T>::operator=(const simpleMatrix<T>& m)
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
-template<class T>
-Foam::simpleMatrix<T> Foam::operator+
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator+
(
- const simpleMatrix<T>& m1,
- const simpleMatrix<T>& m2
+ const simpleMatrix<Type>& m1,
+ const simpleMatrix<Type>& m2
)
{
- return simpleMatrix<T>
+ return simpleMatrix<Type>
(
static_cast<const scalarMatrix&>(m1)
+ static_cast<const scalarMatrix&>(m2),
@@ -131,14 +131,14 @@ Foam::simpleMatrix<T> Foam::operator+
}
-template<class T>
-Foam::simpleMatrix<T> Foam::operator-
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator-
(
- const simpleMatrix<T>& m1,
- const simpleMatrix<T>& m2
+ const simpleMatrix<Type>& m1,
+ const simpleMatrix<Type>& m2
)
{
- return simpleMatrix<T>
+ return simpleMatrix<Type>
(
static_cast<const scalarMatrix&>(m1)
- static_cast<const scalarMatrix&>(m2),
@@ -147,17 +147,17 @@ Foam::simpleMatrix<T> Foam::operator-
}
-template<class T>
-Foam::simpleMatrix<T> Foam::operator*(const scalar s, const simpleMatrix<T>& m)
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator*(const scalar s, const simpleMatrix<Type>& m)
{
- return simpleMatrix<T>(s*m.matrix_, s*m.source_);
+ return simpleMatrix<Type>(s*m.matrix_, s*m.source_);
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
-template<class T>
-Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix<T>& m)
+template<class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix<Type>& m)
{
os << static_cast<const scalarMatrix&>(m) << nl << m.source_;
return os;
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
index d1c411452b8..d7b1a027bec 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<class T>
+template<class Type>
class simpleMatrix;
-template<class T>
-simpleMatrix<T> operator+
+template<class Type>
+simpleMatrix<Type> operator+
(
- const simpleMatrix<T>&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&,
+ const simpleMatrix<Type>&
);
-template<class T>
-simpleMatrix<T> operator-
+template<class Type>
+simpleMatrix<Type> operator-
(
- const simpleMatrix<T>&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&,
+ const simpleMatrix<Type>&
);
-template<class T>
-simpleMatrix<T> operator*
+template<class Type>
+simpleMatrix<Type> operator*
(
const scalar,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&
);
-template<class T>
+template<class Type>
Ostream& operator<<
(
Ostream&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&
);
@@ -81,14 +81,14 @@ Ostream& operator<<
Class simpleMatrix Declaration
\*---------------------------------------------------------------------------*/
-template<class T>
+template<class Type>
class simpleMatrix
:
public scalarMatrix
{
// Private data
- Field<T> source_;
+ Field<Type> source_;
public:
@@ -99,25 +99,25 @@ public:
simpleMatrix(const label);
//- Construct from components
- simpleMatrix(const scalarMatrix&, const Field<T>&);
+ simpleMatrix(const scalarMatrix&, const Field<Type>&);
//- Construct from Istream
simpleMatrix(Istream&);
//- Construct as copy
- simpleMatrix(const simpleMatrix<T>&);
+ simpleMatrix(const simpleMatrix<Type>&);
// Member Functions
// Access
- Field<T>& source()
+ Field<Type>& source()
{
return source_;
}
- const Field<T>& source() const
+ const Field<Type>& source() const
{
return source_;
}
@@ -125,45 +125,45 @@ public:
//- Solve the matrix using Gaussian elimination with pivoting
// and return the solution
- Field<T> solve() const;
+ Field<Type> solve() const;
//- Solve the matrix using LU decomposition with pivoting
// and return the solution
- Field<T> LUsolve() const;
+ Field<Type> LUsolve() const;
// Member Operators
- void operator=(const simpleMatrix<T>&);
+ void operator=(const simpleMatrix<Type>&);
// Friend Operators
- friend simpleMatrix<T> operator+ <T>
+ friend simpleMatrix<Type> operator+ <Type>
(
- const simpleMatrix<T>&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&,
+ const simpleMatrix<Type>&
);
- friend simpleMatrix<T> operator- <T>
+ friend simpleMatrix<Type> operator- <Type>
(
- const simpleMatrix<T>&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&,
+ const simpleMatrix<Type>&
);
- friend simpleMatrix<T> operator* <T>
+ friend simpleMatrix<Type> operator* <Type>
(
const scalar,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&
);
// Ostream Operator
- friend Ostream& operator<< <T>
+ friend Ostream& operator<< <Type>
(
Ostream&,
- const simpleMatrix<T>&
+ const simpleMatrix<Type>&
);
};
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 0c3fae18872..25137c58c7d 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 00000000000..bb4c6c1534c
--- /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<emptyPolyPatch>(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<emptyPolyPatch>(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<Map<label> > 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<label>& globalMap = globalToProc[procI];
+
+ SortableList<label> sorted(globalMap.toc());
+
+ forAll(sorted, i)
+ {
+ Map<label>::iterator iter = globalMap.find(sorted[i]);
+ iter() = i;
+ }
+ }
+ }
+
+
+ // forAll(globalToProc, procI)
+ // {
+ // Pout<< "From processor:" << procI << " want cells/faces:" << endl;
+ // forAllConstIter(Map<label>, globalToProc[procI], iter)
+ // {
+ // Pout<< " global:" << iter.key()
+ // << " local:" << globalNumbering.toLocal(procI, iter.key())
+ // << endl;
+ // }
+ // Pout<< endl;
+ // }
+ }
+
+
+ // 2. The overall compact addressing is
+ // - myProcNo first
+ // - all other processors consecutively
+
+ labelList compactStart(Pstream::nProcs());
+ compactStart[Pstream::myProcNo()] = 0;
+ label nCompact = mesh.nCells()+nBnd;
+ forAll(compactStart, procI)
+ {
+ if (procI != Pstream::myProcNo())
+ {
+ compactStart[procI] = nCompact;
+ nCompact += globalToProc[procI].size();
+
+ // Pout<< "Data wanted from " << procI << " starts at "
+ // << compactStart[procI] << endl;
+ }
+ }
+ // Pout<< "Overall cells needed:" << nCompact << endl;
+
+
+ // 3. Find out what to receive/send in compact addressing.
+ labelListList recvCompact(Pstream::nProcs());
+ for (label procI = 0; procI < Pstream::nProcs(); procI++)
+ {
+ if (procI != Pstream::myProcNo())
+ {
+ labelList wantedGlobals(globalToProc[procI].size());
+ recvCompact[procI].setSize(globalToProc[procI].size());
+
+ label i = 0;
+ forAllConstIter(Map<label>, globalToProc[procI], iter)
+ {
+ wantedGlobals[i] = iter.key();
+ recvCompact[procI][i] = compactStart[procI]+iter();
+ i++;
+ }
+
+ // Pout<< "From proc:" << procI
+ // << " I need (globalcells):" << wantedGlobals
+ // << " which are my compact:" << recvCompact[procI]
+ // << endl;
+
+ // Send the global cell numbers I need from procI
+ OPstream str(Pstream::blocking, procI);
+ str << wantedGlobals;
+ }
+ else
+ {
+ recvCompact[procI] =
+ compactStart[procI]
+ + identity(mesh.nCells()+nBnd);
+ }
+ }
+ labelListList sendCompact(Pstream::nProcs());
+ for (label procI = 0; procI < Pstream::nProcs(); procI++)
+ {
+ if (procI != Pstream::myProcNo())
+ {
+ // See what neighbour wants to receive (= what I need to send)
+
+ IPstream str(Pstream::blocking, procI);
+ labelList globalCells(str);
+
+ labelList& procCompact = sendCompact[procI];
+ procCompact.setSize(globalCells.size());
+
+ // Convert from globalCells (all on my processor!) into compact
+ // addressing
+ forAll(globalCells, i)
+ {
+ label cellI = globalNumbering.toLocal(globalCells[i]);
+ procCompact[i] = compactStart[Pstream::myProcNo()]+cellI;
+ }
+ }
+ else
+ {
+ sendCompact[procI] = recvCompact[procI];
+ }
+ }
+
+ // Convert stencil to compact numbering
+ forAll(stencil_, faceI)
+ {
+ labelList& stencilCells = stencil_[faceI];
+
+ forAll(stencilCells, i)
+ {
+ label globalCellI = stencilCells[i];
+ label procI = globalNumbering.whichProcID(globalCellI);
+ if (procI != Pstream::myProcNo())
+ {
+ label localCompact = globalToProc[procI][globalCellI];
+ stencilCells[i] = compactStart[procI]+localCompact;
+ }
+ else
+ {
+ label localCompact = globalNumbering.toLocal(globalCellI);
+ stencilCells[i] = compactStart[procI]+localCompact;
+ }
+
+ }
+ }
+ // Pout<< "***stencil_:" << stencil_ << endl;
+
+ // Constuct map for distribution of compact data.
+ mapPtr_.reset
+ (
+ new mapDistribute
+ (
+ nCompact,
+ sendCompact,
+ recvCompact,
+ true // reuse send/recv maps.
+ )
+ );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::extendedStencil::extendedStencil
+(
+ const mapDistribute& map,
+ const labelListList& stencil
+)
+:
+ mapPtr_
+ (
+ autoPtr<mapDistribute>
+ (
+ new mapDistribute
+ (
+ map.constructSize(),
+ map.subMap(),
+ map.constructMap()
+ )
+ )
+ ),
+ stencil_(stencil)
+{}
+
+
+Foam::extendedStencil::extendedStencil(const polyMesh& mesh)
+{
+ calcExtendedFaceStencil(mesh);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
new file mode 100644
index 00000000000..1d00972403b
--- /dev/null
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::extendedStencil
+
+Description
+ Calculates/constains the extended face stencil.
+
+ The stencil is a list of indices into either cells or boundary faces
+ in a compact way. (element 0 is owner, 1 is neighbour). The index numbering
+ is
+ - cells first
+ - then all (non-empty patch) boundary faces
+
+ When used in evaluation is a two stage process:
+ - collect the data (cell data and non-empty boundaries) into a
+ single field
+ - (parallel) distribute the field
+ - sum the weights*field.
+
+SourceFiles
+ extendedStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedStencil_H
+#define extendedStencil_H
+
+#include "mapDistribute.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndex;
+
+/*---------------------------------------------------------------------------*\
+ Class extendedStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedStencil
+{
+ // Private data
+
+ //- Swap map for getting neigbouring data
+ autoPtr<mapDistribute> mapPtr_;
+
+ //- Per face the stencil.
+ labelListList stencil_;
+
+
+ // Private Member Functions
+
+ void calcFaceStencils(const polyMesh&, const globalIndex&);
+
+ //- Calculate the stencil (but not weights)
+ void calcExtendedFaceStencil(const polyMesh&);
+
+
+ //- Disallow default bitwise copy construct
+ extendedStencil(const extendedStencil&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const extendedStencil&);
+
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ extendedStencil(const mapDistribute& map, const labelListList&);
+
+ //- Construct from all cells and boundary faces
+ extendedStencil(const polyMesh&);
+
+
+
+ // Member Functions
+
+ //- Return reference to the parallel distribution map
+ const mapDistribute& map() const
+ {
+ return mapPtr_();
+ }
+
+ //- Return reference to the stencil
+ const labelListList& stencil() const
+ {
+ return stencil_;
+ }
+
+ //- Use map to get the data into stencil order
+ template<class T>
+ void collectData
+ (
+ const GeometricField<T, fvPatchField, volMesh>& fld,
+ List<List<T> >& stencilFld
+ ) const;
+
+ //- Given weights interpolate vol field
+ template<class Type>
+ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
+ (
+ const GeometricField<Type, fvPatchField, volMesh>& fld,
+ const List<List<scalar> >& stencilWeights
+ ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "extendedStencilTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
new file mode 100644
index 00000000000..a408621dc4c
--- /dev/null
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::extendedStencil::collectData
+(
+ const GeometricField<Type, fvPatchField, volMesh>& fld,
+ List<List<Type> >& stencilFld
+) const
+{
+ // 1. Construct cell data in compact addressing
+ List<Type> compactFld(map().constructSize(), pTraits<Type>::zero);
+
+ // Insert my internal values
+ forAll(fld, cellI)
+ {
+ compactFld[cellI] = fld[cellI];
+ }
+ // Insert my boundary values
+ label nCompact = fld.size();
+ forAll(fld.boundaryField(), patchI)
+ {
+ const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];
+
+ forAll(pfld, i)
+ {
+ compactFld[nCompact++] = pfld[i];
+ }
+ }
+
+ // Do all swapping
+ map().distribute(compactFld);
+
+ // 2. Pull to stencil
+ stencilFld.setSize(stencil_.size());
+
+ forAll(stencil_, faceI)
+ {
+ const labelList& compactCells = stencil_[faceI];
+
+ stencilFld[faceI].setSize(compactCells.size());
+
+ forAll(compactCells, i)
+ {
+ stencilFld[faceI][i] = compactFld[compactCells[i]];
+ }
+ }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::extendedStencil::interpolate
+(
+ const GeometricField<Type, fvPatchField, volMesh>& fld,
+ const List<List<scalar> >& stencilWeights
+) const
+{
+ const fvMesh& mesh = fld.mesh();
+
+ // Collect internal and boundary values
+ List<List<Type> > stencilFld;
+ collectData(fld, stencilFld);
+
+ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+ (
+ new GeometricField<Type, fvsPatchField, surfaceMesh>
+ (
+ IOobject
+ (
+ fld.name(),
+ mesh.time().timeName(),
+ mesh
+ ),
+ mesh,
+ dimensioned<Type>
+ (
+ fld.name(),
+ fld.dimensions(),
+ pTraits<Type>::zero
+ )
+ )
+ );
+ GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsfCorr();
+
+ for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+ {
+ const List<Type>& stField = stencilFld[faceI];
+ const List<scalar>& stWeight = stencilWeights[faceI];
+
+ forAll(stField, i)
+ {
+ sf[faceI] += stField[i]*stWeight[i];
+ }
+ }
+ // And what for boundaries?
+
+ return tsfCorr;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
new file mode 100644
index 00000000000..a215576330b
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "quadraticFit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ makeSurfaceInterpolationScheme(quadraticFit);
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
new file mode 100644
index 00000000000..412b454015d
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ quadraticFit
+
+Description
+ Quadratic fit interpolation scheme which applies an explicit correction to
+ linear.
+
+SourceFiles
+ quadraticFit.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFit_H
+#define quadraticFit_H
+
+#include "linear.H"
+#include "quadraticFitData.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class quadraticFit Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class quadraticFit
+:
+ public linear<Type>
+{
+ // Private Data
+ const scalar centralWeight_;
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ quadraticFit(const quadraticFit&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const quadraticFit&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("quadraticFit");
+
+
+ // Constructors
+
+ //- Construct from mesh and Istream
+ quadraticFit(const fvMesh& mesh, Istream& is)
+ :
+ linear<Type>(mesh),
+ centralWeight_(readScalar(is))
+ {}
+
+
+ //- Construct from mesh, faceFlux and Istream
+ quadraticFit
+ (
+ const fvMesh& mesh,
+ const surfaceScalarField& faceFlux,
+ Istream& is
+ )
+ :
+ linear<Type>(mesh),
+ centralWeight_(readScalar(is))
+ {}
+
+
+ // Member Functions
+
+ //- Return true if this scheme uses an explicit correction
+ virtual bool corrected() const
+ {
+ return true;
+ }
+
+ //- Return the explicit correction to the face-interpolate
+ virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+ correction
+ (
+ const GeometricField<Type, fvPatchField, volMesh>& vf
+ ) const
+ {
+ const fvMesh& mesh = this->mesh();
+
+ const quadraticFitData& cfd = quadraticFitData::New
+ (
+ mesh,
+ centralWeight_
+ );
+
+ const extendedStencil& stencil = cfd.stencil();
+ const List<scalarList>& f = cfd.fit();
+
+ return stencil.interpolate(vf, f);
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
new file mode 100644
index 00000000000..fd0e55fcdd6
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
@@ -0,0 +1,310 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "quadraticFitData.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "SVD.H"
+#include "syncTools.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(quadraticFitData, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::quadraticFitData::quadraticFitData
+(
+ const fvMesh& mesh,
+ const scalar cWeight
+)
+:
+ MeshObject<fvMesh, quadraticFitData>(mesh),
+ centralWeight_(cWeight),
+#ifdef SPHERICAL_GEOMETRY
+ dim_(2),
+#else
+ dim_(mesh.nGeometricD()),
+#endif
+ minSize_
+ (
+ dim_ == 1 ? 3 :
+ dim_ == 2 ? 6 :
+ dim_ == 3 ? 9 : 0
+ ),
+ stencil_(mesh),
+ fit_(mesh.nInternalFaces())
+{
+ if (debug)
+ {
+ Info << "Contructing quadraticFitData" << endl;
+ }
+
+ // check input
+ if (centralWeight_ < 1 - SMALL)
+ {
+ FatalErrorIn("quadraticFitData::quadraticFitData")
+ << "centralWeight requested = " << centralWeight_
+ << " should not be less than one"
+ << exit(FatalError);
+ }
+
+ if (minSize_ == 0)
+ {
+ FatalErrorIn("quadraticFitSnGradData")
+ << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
+ }
+
+ // store the polynomial size for each cell to write out
+ surfaceScalarField interpPolySize
+ (
+ IOobject
+ (
+ "quadraticFitInterpPolySize",
+ "constant",
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("quadraticFitInterpPolySize", dimless, scalar(0))
+ );
+
+ // Get the cell/face centres in stencil order.
+ // Centred face stencils no good for triangles of tets. Need bigger stencils
+ List<List<point> > stencilPoints(stencil_.stencil().size());
+ stencil_.collectData
+ (
+ mesh.C(),
+ stencilPoints
+ );
+
+ // find the fit coefficients for every face in the mesh
+
+ for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
+ {
+ interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
+ }
+ interpPolySize.write();
+ if (debug)
+ {
+ Info<< "quadraticFitData::quadraticFitData() :"
+ << "Finished constructing polynomialFit data"
+ << endl;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::quadraticFitData::findFaceDirs
+(
+ vector& idir, // value changed in return
+ vector& jdir, // value changed in return
+ vector& kdir, // value changed in return
+ const fvMesh& mesh,
+ const label faci
+)
+{
+ idir = mesh.Sf()[faci];
+ idir /= mag(idir);
+
+# ifndef SPHERICAL_GEOMETRY
+ if (mesh.nGeometricD() <= 2) // find the normal direcion
+ {
+ if (mesh.directions()[0] == -1)
+ {
+ kdir = vector(1, 0, 0);
+ }
+ else if (mesh.directions()[1] == -1)
+ {
+ kdir = vector(0, 1, 0);
+ }
+ else
+ {
+ kdir = vector(0, 0, 1);
+ }
+ }
+ else // 3D so find a direction in the place of the face
+ {
+ const face& f = mesh.faces()[faci];
+ kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
+ }
+# else
+ // Spherical geometry so kdir is the radial direction
+ kdir = mesh.Cf()[faci];
+# endif
+
+ if (mesh.nGeometricD() == 3)
+ {
+ // Remove the idir component from kdir and normalise
+ kdir -= (idir & kdir)*idir;
+
+ scalar magk = mag(kdir);
+
+ if (magk < SMALL)
+ {
+ FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+ << exit(FatalError);
+ }
+ else
+ {
+ kdir /= magk;
+ }
+ }
+
+ jdir = kdir ^ idir;
+}
+
+
+Foam::label Foam::quadraticFitData::calcFit
+(
+ const List<point>& C,
+ const label faci
+)
+{
+ vector idir(1,0,0);
+ vector jdir(0,1,0);
+ vector kdir(0,0,1);
+ findFaceDirs(idir, jdir, kdir, mesh(), faci);
+
+ scalarList wts(C.size(), scalar(1));
+ wts[0] = centralWeight_;
+ wts[1] = centralWeight_;
+
+ point p0 = mesh().faceCentres()[faci];
+ scalar scale = 0;
+
+ // calculate the matrix of the polynomial components
+ Matrix<scalar> B(C.size(), minSize_, scalar(0));
+
+ for(label ip = 0; ip < C.size(); ip++)
+ {
+ const point& p = C[ip];
+
+ scalar px = (p - p0)&idir;
+ scalar py = (p - p0)&jdir;
+#ifndef SPHERICAL_GEOMETRY
+ scalar pz = (p - p0)&kdir;
+#else
+ scalar pz = mag(p) - mag(p0);
+#endif
+
+ if (ip == 0)
+ {
+ scale = max(max(mag(px), mag(py)), mag(pz));
+ }
+
+ px /= scale;
+ py /= scale;
+ pz /= scale;
+
+ label is = 0;
+ B[ip][is++] = wts[ip];
+
+ B[ip][is++] = wts[ip]*px;
+ B[ip][is++] = wts[ip]*sqr(px);
+
+ if (dim_ >= 2)
+ {
+ B[ip][is++] = wts[ip]*py;
+ B[ip][is++] = wts[ip]*px*py;
+ B[ip][is++] = wts[ip]*sqr(py);
+ }
+ if (dim_ == 3)
+ {
+ B[ip][is++] = wts[ip]*pz;
+ B[ip][is++] = wts[ip]*px*pz;
+ //B[ip][is++] = wts[ip]*py*pz;
+ B[ip][is++] = wts[ip]*sqr(pz);
+ }
+ }
+
+ // Set the fit
+ label stencilSize = C.size();
+ fit_[faci].setSize(stencilSize);
+ scalarList singVals(minSize_);
+ label nSVDzeros = 0;
+
+ bool goodFit = false;
+ for(scalar tol = SMALL; tol < 0.1 && !goodFit; tol *= 10)
+ {
+ SVD svd(B, tol);
+
+ scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
+ scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
+
+ goodFit = sign(fit0) == sign(fit1);
+
+ if (goodFit)
+ {
+ fit_[faci][0] = fit0;
+ fit_[faci][1] = fit1;
+ for(label i = 2; i < stencilSize; i++)
+ {
+ fit_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
+ }
+
+ singVals = svd.S();
+ nSVDzeros = svd.nZeros();
+ }
+ }
+ if (!goodFit)
+ {
+ FatalErrorIn
+ (
+ "quadraticFitData::calcFit(const pointField&, const label)"
+ ) << "For face " << faci << endl
+ << "Fit not good even once tol >= 0.1"
+ << exit(FatalError);
+ }
+
+ const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
+ mesh().surfaceInterpolation::weights();
+
+ // remove the uncorrected linear coefficients
+
+ fit_[faci][0] -= w[faci];
+ fit_[faci][1] -= 1 - w[faci];
+
+ return minSize_ - nSVDzeros;
+}
+
+
+bool Foam::quadraticFitData::movePoints()
+{
+ notImplemented("quadraticFitData::movePoints()");
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H
new file mode 100644
index 00000000000..19859ea621c
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ quadraticFitData
+
+Description
+ Data for the quadratic fit correction interpolation scheme
+
+SourceFiles
+ quadraticFitData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFitData_H
+#define quadraticFitData_H
+
+#include "MeshObject.H"
+#include "fvMesh.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndex;
+
+/*---------------------------------------------------------------------------*\
+ Class quadraticFitData Declaration
+\*---------------------------------------------------------------------------*/
+
+class quadraticFitData
+:
+ public MeshObject<fvMesh, quadraticFitData>
+{
+ // Private data
+
+ //- weights for central stencil
+ const scalar centralWeight_;
+
+ //- dimensionality of the geometry
+ const label dim_;
+
+ //- minimum stencil size
+ const label minSize_;
+
+ //- Extended stencil addressing
+ extendedStencil stencil_;
+
+ //- For each cell in the mesh store the values which multiply the
+ // values of the stencil to obtain the gradient for each direction
+ List<scalarList> fit_;
+
+
+ // Private member functions
+
+ //- Find the normal direction and i, j and k directions for face faci
+ static void findFaceDirs
+ (
+ vector& idir, // value changed in return
+ vector& jdir, // value changed in return
+ vector& kdir, // value changed in return
+ const fvMesh& mesh,
+ const label faci
+ );
+
+ label calcFit(const List<point>&, const label faci);
+
+
+public:
+
+ TypeName("quadraticFitData");
+
+
+ // Constructors
+
+ explicit quadraticFitData
+ (
+ const fvMesh& mesh,
+ scalar cWeightDim
+ );
+
+
+ // Destructor
+
+ virtual ~quadraticFitData()
+ {}
+
+
+ // Member functions
+
+
+ //- Return reference to the stencil
+ const extendedStencil& stencil() const
+ {
+ return stencil_;
+ }
+
+ //- Return reference to fit coefficients
+ const List<scalarList>& fit() const
+ {
+ return fit_;
+ }
+
+ //- Delete the data when the mesh moves not implemented
+ virtual bool movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--
GitLab