zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Added support for quadraticFit interpolation scheme.

Browse Source 
Added the quadraticFit interpolation scheme.
This commit is contained in:
henry
2008-09-27 09:25:30 +01:00
parent 63d5360166
commit 26f2c24d84
24 changed files with 2722 additions and 249 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/OpenFOAM/Make/files
View File

@ -179,6 +179,7 @@ matrices/solution/solution.C
scalarMatrix = matrices/scalarMatrix
$(scalarMatrix)/scalarMatrix.C
$(scalarMatrix)/SVD/SVD.C
LUscalarMatrix = matrices/LUscalarMatrix
$(LUscalarMatrix)/LUscalarMatrix.C

99
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C Normal file
View File

@ -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;
}
// ************************************************************************* //

102
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H Normal file
View File

@ -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
// ************************************************************************* //

14
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
View File

@ -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)
);
}
}

174
src/OpenFOAM/containers/Matrix/Matrix.C → src/OpenFOAM/matrices/Matrix/Matrix.C
View File

@ -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++)

83
src/OpenFOAM/containers/Matrix/Matrix.H → src/OpenFOAM/matrices/Matrix/Matrix.H
View File

@ -51,25 +51,37 @@ namespace Foam
// Forward declaration of friend functions and operators
template<class T> class Matrix;
template<class Type> class Matrix;
template<class T> const T& max(const Matrix<T>&);
template<class T> const T& min(const Matrix<T>&);
template<class Type> const Type& max(const Matrix<Type>&);
template<class Type> const Type& min(const Matrix<Type>&);
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 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 T> Istream& operator>>(Istream&, Matrix<T>&);
template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
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>&);
};

52
src/OpenFOAM/containers/Matrix/MatrixI.H → src/OpenFOAM/matrices/Matrix/MatrixI.H
View File

@ -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);

44
src/OpenFOAM/containers/Matrix/MatrixIO.C → src/OpenFOAM/matrices/Matrix/MatrixIO.C
View File

@ -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));
}
}

400
src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C Normal file
View File

@ -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;
}
*/
}
// ************************************************************************* //

128
src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H Normal file
View File

@ -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
// ************************************************************************* //

75
src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H Normal file
View File

@ -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;
}
// ************************************************************************* //

159
src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
View File

@ -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();
}
// ************************************************************************* //

46
src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
View File

@ -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

26
src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
View File

@ -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());

68
src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
View File

@ -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;

66
src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
View File

@ -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>&
);
};

5
src/finiteVolume/Make/files
View File

@ -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

525
src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C Normal file
View File

@ -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);
}
// ************************************************************************* //

151
src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H Normal file
View File

@ -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
// ************************************************************************* //

129
src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C Normal file
View File

@ -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;
}
// ************************************************************************* //

36
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C Normal file
View File

@ -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);
}
// ************************************************************************* //

138
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H Normal file
View File

@ -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
// ************************************************************************* //

310
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C Normal file
View File

@ -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;
}
// ************************************************************************* //

140
src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H Normal file
View File

@ -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
// ************************************************************************* //