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

Merge branch 'master' of ssh://hunt/~OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2008-09-29 10:14:50 +01:00
parent 5a9cfc954b 318b71e206
commit 2b163418bc
47 changed files with 3479 additions and 618 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
View File

@ -139,7 +139,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
IOobject LESPropertiesHeader
(
"RASProperties",
"LESProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,

4
src/ODE/ODE/ODE.H
View File

@ -34,7 +34,7 @@ Description
#define ODE_H
#include "scalarField.H"
#include "Matrix.H"
#include "scalarMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,7 +80,7 @@ public:
const scalar x,
const scalarField& y,
scalarField& dfdx,
Matrix<scalar>& dfdy
scalarSquareMatrix& dfdy
) const = 0;
};

19
src/ODE/ODESolvers/KRR4/KRR4.C
View File

@ -25,7 +25,6 @@ License
\*---------------------------------------------------------------------------*/
#include "KRR4.H"
#include "simpleMatrix.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -36,11 +35,11 @@ namespace Foam
{
addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
const scalar
const scalar
KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
const scalar
const scalar
KRR4::gamma = 1.0/2.0,
KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
@ -81,8 +80,8 @@ void Foam::KRR4::solve
const ODE& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
@ -109,14 +108,14 @@ void Foam::KRR4::solve
a_[i][i] += 1.0/(gamma*h);
}
simpleMatrix<scalar>::LUDecompose(a_, pivotIndices_);
LUDecompose(a_, pivotIndices_);
for (register label i=0; i<n_; i++)
{
g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
}
simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g1_);
LUBacksubstitute(a_, pivotIndices_, g1_);
for (register label i=0; i<n_; i++)
{
@ -131,7 +130,7 @@ void Foam::KRR4::solve
g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
}
simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g2_);
LUBacksubstitute(a_, pivotIndices_, g2_);
for (register label i=0; i<n_; i++)
{
@ -146,7 +145,7 @@ void Foam::KRR4::solve
g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
}
simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g3_);
LUBacksubstitute(a_, pivotIndices_, g3_);
for (register label i=0; i<n_; i++)
{
@ -154,7 +153,7 @@ void Foam::KRR4::solve
+ (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
}
simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g4_);
LUBacksubstitute(a_, pivotIndices_, g4_);
for (register label i=0; i<n_; i++)
{

6
src/ODE/ODESolvers/KRR4/KRR4.H
View File

@ -62,15 +62,15 @@ class KRR4
mutable scalarField g4_;
mutable scalarField yErr_;
mutable scalarField dfdx_;
mutable Matrix<scalar> dfdy_;
mutable Matrix<scalar> a_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
static const int maxtry = 40;
static const scalar safety, grow, pgrow, shrink, pshrink, errcon;
static const scalar
static const scalar
gamma,
a21, a31, a32,
c21, c31, c32, c41, c42, c43,

10
src/ODE/ODESolvers/SIBS/SIBS.H
View File

@ -60,8 +60,8 @@ class SIBS
static const scalar safe1, safe2, redMax, redMin, scaleMX;
mutable scalarField a_;
mutable Matrix<scalar> alpha_;
mutable Matrix<scalar> d_p_;
mutable scalarSquareMatrix alpha_;
mutable scalarSquareMatrix d_p_;
mutable scalarField x_p_;
mutable scalarField err_;
@ -69,7 +69,7 @@ class SIBS
mutable scalarField ySeq_;
mutable scalarField yErr_;
mutable scalarField dfdx_;
mutable Matrix<scalar> dfdy_;
mutable scalarSquareMatrix dfdy_;
mutable label first_, kMax_, kOpt_;
mutable scalar epsOld_, xNew_;
@ -82,7 +82,7 @@ void SIMPR
const scalarField& y,
const scalarField& dydx,
const scalarField& dfdx,
const Matrix<scalar>& dfdy,
const scalarSquareMatrix& dfdy,
const scalar deltaX,
const label nSteps,
scalarField& yEnd
@ -97,7 +97,7 @@ void polyExtrapolate
scalarField& yz,
scalarField& dy,
scalarField& x_p,
Matrix<scalar>& d_p
scalarSquareMatrix& d_p
) const;

13
src/ODE/ODESolvers/SIBS/SIMPR.C
View File

@ -25,7 +25,6 @@ License
\*---------------------------------------------------------------------------*/
#include "SIBS.H"
#include "simpleMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -36,7 +35,7 @@ void Foam::SIBS::SIMPR
const scalarField& y,
const scalarField& dydx,
const scalarField& dfdx,
const Matrix<scalar>& dfdy,
const scalarSquareMatrix& dfdy,
const scalar deltaX,
const label nSteps,
scalarField& yEnd
@ -44,7 +43,7 @@ void Foam::SIBS::SIMPR
{
scalar h = deltaX/nSteps;
Matrix<scalar> a(n_, n_);
scalarSquareMatrix a(n_);
for (register label i=0; i<n_; i++)
{
for (register label j=0; j<n_; j++)
@ -55,14 +54,14 @@ void Foam::SIBS::SIMPR
}
labelList pivotIndices(n_);
simpleMatrix<scalar>::LUDecompose(a, pivotIndices);
LUDecompose(a, pivotIndices);
for (register label i=0; i<n_; i++)
{
yEnd[i] = h*(dydx[i] + h*dfdx[i]);
}
simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
LUBacksubstitute(a, pivotIndices, yEnd);
scalarField del(yEnd);
scalarField ytemp(n_);
@ -83,7 +82,7 @@ void Foam::SIBS::SIMPR
yEnd[i] = h*yEnd[i] - del[i];
}
simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
LUBacksubstitute(a, pivotIndices, yEnd);
for (register label i=0; i<n_; i++)
{
@ -99,7 +98,7 @@ void Foam::SIBS::SIMPR
yEnd[i] = h*yEnd[i] - del[i];
}
simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
LUBacksubstitute(a, pivotIndices, yEnd);
for (register label i=0; i<n_; i++)
{

2
src/ODE/ODESolvers/SIBS/polyExtrapolate.C
View File

@ -36,7 +36,7 @@ void Foam::SIBS::polyExtrapolate
scalarField& yz,
scalarField& dy,
scalarField& x,
Matrix<scalar>& d
scalarSquareMatrix& d
) const
{
label n = yz.size();

121
src/OSspecific/Unix/regularExpression.H Normal file
View File

@ -0,0 +1,121 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::regularExpression
Description
Wrapper around regular expressions.
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef regularExpression_H
#define regularExpression_H
#include <sys/types.h>
#include <regex.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class regularExpression Declaration
\*---------------------------------------------------------------------------*/
class regularExpression
{
// Private data
//- Precompiled regular expression
regex_t* preg_;
// Private member functions
//- Disallow default bitwise copy construct
regularExpression(const regularExpression&);
//- Disallow default bitwise assignment
void operator=(const regularExpression&);
public:
// Constructors
//- Construct from string
inline regularExpression(const string& s)
{
preg_ = new regex_t;
if (regcomp(preg_, s.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
{
FatalErrorIn
(
"regularExpression::regularExpression(const char*)"
) << "Failed to compile regular expression " << s
<< exit(FatalError);
}
}
// Destructor
//- Construct from string
inline ~regularExpression()
{
if (preg_)
{
regfree(preg_);
delete preg_;
}
}
// Member functions
//- Matches?
inline bool matches(const string& s)
{
size_t nmatch = 0;
regmatch_t *pmatch = NULL;
return regexec(preg_, s.c_str(), nmatch, pmatch, 0) == 0;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

5
src/OpenFOAM/Make/files
View File

@ -177,8 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C
matrices/solution/solution.C
scalarMatrix = matrices/scalarMatrix
$(scalarMatrix)/scalarMatrix.C
scalarMatrices = matrices/scalarMatrices
$(scalarMatrices)/scalarMatrices.C
$(scalarMatrices)/SVD/SVD.C
LUscalarMatrix = matrices/LUscalarMatrix
$(LUscalarMatrix)/LUscalarMatrix.C

2
src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C
View File

@ -22,8 +22,6 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/
#include "Pstream.H"

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

@ -0,0 +1,100 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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>
template<class Form>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, 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;
}
// ************************************************************************* //

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

@ -0,0 +1,103 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Form, class Type> class Matrix;
/*---------------------------------------------------------------------------*\
Class DiagonalMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class DiagonalMatrix
:
public List<Type>
{
public:
// Constructors
//- Construct from diagonal component of a Matrix
template<class Form>
DiagonalMatrix<Type>(const Matrix<Form, 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
// ************************************************************************* //

8
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
View File

@ -31,9 +31,9 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::LUscalarMatrix::LUscalarMatrix(const Matrix<scalar>& matrix)
Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
:
scalarMatrix(matrix),
scalarSquareMatrix(matrix),
pivotIndices_(n())
{
LUDecompose(*this, pivotIndices_);
@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
nCells += lduMatrices[i].size();
}
Matrix<scalar> m(nCells, nCells, 0.0);
scalarSquareMatrix m(nCells, 0.0);
transfer(m);
convert(lduMatrices);
}
@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
else
{
label nCells = ldum.lduAddr().size();
Matrix<scalar> m(nCells, nCells, 0.0);
scalarSquareMatrix m(nCells, 0.0);
transfer(m);
convert(ldum, interfaceCoeffs, interfaces);
}

10
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
View File

@ -36,7 +36,7 @@ SourceFiles
#ifndef LUscalarMatrix_H
#define LUscalarMatrix_H
#include "scalarMatrix.H"
#include "scalarMatrices.H"
#include "labelList.H"
#include "FieldField.H"
#include "lduInterfaceFieldPtrsList.H"
@ -55,7 +55,7 @@ class procLduMatrix;
class LUscalarMatrix
:
public scalarMatrix
public scalarSquareMatrix
{
// Private data
@ -78,7 +78,7 @@ class LUscalarMatrix
void convert(const PtrList<procLduMatrix>& lduMatrices);
//- Print the ratio of the mag-sum of the off-diagonal coefficients
//- Print the ratio of the mag-sum of the off-diagonal coefficients
// to the mag-diagonal
void printDiagonalDominance() const;
@ -87,8 +87,8 @@ public:
// Constructors
//- Construct from Matrix<scalar> and perform LU decomposition
LUscalarMatrix(const Matrix<scalar>&);
//- Construct from scalarSquareMatrix and perform LU decomposition
LUscalarMatrix(const scalarSquareMatrix&);
//- Construct from lduMatrix and perform LU decomposition
LUscalarMatrix

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

176
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 Form, class Type>
void Foam::Matrix<Form, 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 Form, class Type>
Foam::Matrix<Form, 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 Form, class Type>
const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
{
Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
Matrix<Form, Type>* nullPtr = reinterpret_cast<Matrix<Form, Type>*>(NULL);
return *nullPtr;
}
template<class T>
Foam::Matrix<T>::Matrix(const label n, const label m)
template<class Form, class Type>
Foam::Matrix<Form, 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<Form, 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 Form, class Type>
Foam::Matrix<Form, 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<Form, 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 Form, class Type>
Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, 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 Form, class Type>
void Foam::Matrix<Form, 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 Form, class Type>
void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
{
clear();
@ -169,14 +169,32 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
}
template<class Form, class Type>
Form Foam::Matrix<Form, Type>::T() const
{
const Matrix<Form, Type>& A = *this;
Form At(m(), n());
for (register label i=0; i<n(); i++)
{
for (register label j=0; j<m(); j++)
{
At[j][i] = A[i][j];
}
}
return At;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class T>
void Foam::Matrix<T>::operator=(const T& t)
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator=(const Type& t)
{
if (v_)
{
T* v = v_[0];
Type* v = v_[0];
label nm = n_*m_;
for (register label i=0; i<nm; i++)
@ -188,12 +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 Form, class Type>
void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
{
if (this == &a)
{
FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
FatalErrorIn("Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)")
<< "attempted assignment to self"
<< abort(FatalError);
}
@ -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 Form, class Type>
const Type& Foam::max(const Matrix<Form, 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<Form, 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 Form, class Type>
const Type& Foam::min(const Matrix<Form, 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<Form, 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 Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a)
{
Matrix<T> na(a.n_, a.m_);
Form na(a.n_, a.m_);
if (a.n_ && a.m_)
{
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 Form, class Type>
Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, 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<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of columns: "
FatalErrorIn
(
"Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_
<< abort(FatalError);
}
Matrix<T> ab(a.n_, a.m_);
Form ab(a.n_, a.m_);
T* abv = ab.v_[0];
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 Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, 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<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of columns: "
FatalErrorIn
(
"Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_
<< abort(FatalError);
}
Matrix<T> ab(a.n_, a.m_);
Form ab(a.n_, a.m_);
T* abv = ab.v_[0];
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,17 +405,17 @@ 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 Form, class Type>
Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
{
Matrix<T> sa(a.n_, a.m_);
Form sa(a.n_, a.m_);
if (a.n_ && a.m_)
{
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_;;
label nm = a.n_*a.m_;
for (register label i=0; i<nm; i++)
{
sav[i] = s*av[i];

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

@ -51,25 +51,26 @@ namespace Foam
// Forward declaration of friend functions and operators
template<class T> class Matrix;
template<class Form, class Type> class Matrix;
template<class T> const T& max(const Matrix<T>&);
template<class T> const T& min(const Matrix<T>&);
template<class Form, class Type> Istream& operator>>
(
Istream&,
Matrix<Form, 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 T> Istream& operator>>(Istream&, Matrix<T>&);
template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
template<class Form, class Type> Ostream& operator<<
(
Ostream&,
const Matrix<Form, Type>&
);
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
\*---------------------------------------------------------------------------*/
template<class T>
template<class Form, class Type>
class Matrix
{
// Private data
@ -78,7 +79,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 +98,16 @@ public:
//- Construct with given number of rows and columns
// and value for all elements.
Matrix(const label n, const label m, const T&);
Matrix(const label n, const label m, const Type&);
//- Copy constructor.
Matrix(const Matrix<T>&);
Matrix(const Matrix<Form, Type>&);
//- Construct from Istream.
Matrix(Istream&);
//- Clone
inline autoPtr<Matrix<T> > clone() const;
inline autoPtr<Matrix<Form, Type> > clone() const;
// Destructor
@ -117,7 +118,7 @@ public:
// Member functions
//- Return a null Matrix
static const Matrix<T>& null();
static const Matrix<Form, Type>& null();
// Access
@ -148,48 +149,64 @@ public:
//- Transfer the contents of the argument Matrix into this Matrix
// and annull the argument Matrix.
void transfer(Matrix<T>&);
void transfer(Matrix<Form, Type>&);
//- Return the transpose of the matrix
Form T() const;
// Member operators
//- Return subscript-checked element of Matrix.
inline T* operator[](const label);
inline Type* operator[](const label);
//- Return subscript-checked element of constant Matrix.
inline const T* operator[](const label) const;
inline const Type* operator[](const label) const;
//- Assignment operator. Takes linear time.
void operator=(const Matrix<T>&);
void operator=(const Matrix<Form, Type>&);
//- Assignment of all entries to the given value
void operator=(const T&);
// Friend functions
friend const T& max<T>(const Matrix<T>&);
friend const T& min<T>(const Matrix<T>&);
// Friend operators
friend Matrix<T> operator-<T>(const Matrix<T>&);
friend Matrix<T> operator+<T>(const Matrix<T>&, const Matrix<T>&);
friend Matrix<T> operator-<T>(const Matrix<T>&, const Matrix<T>&);
friend Matrix<T> operator*<T>(const scalar, const Matrix<T>&);
void operator=(const Type&);
// IOstream operators
//- Read Matrix from Istream, discarding contents of existing Matrix.
friend Istream& operator>> <T>(Istream&, Matrix<T>&);
friend Istream& operator>> <Form, Type>(Istream&, Matrix<Form, Type>&);
// Write Matrix to Ostream.
friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
friend Ostream& operator<< <Form, Type>(Ostream&, const Matrix<Form, Type>&);
};
// Global functions and operators
template<class Form, class Type> const Type& max(const Matrix<Form, Type>&);
template<class Form, class Type> const Type& min(const Matrix<Form, Type>&);
template<class Form, class Type> Form operator-(const Matrix<Form, Type>&);
template<class Form, class Type> Form operator+
(
const Matrix<Form, Type>&,
const Matrix<Form, Type>&
);
template<class Form, class Type> Form operator-
(
const Matrix<Form, Type>&,
const Matrix<Form, Type>&
);
template<class Form, class Type> Form operator*
(
const scalar,
const Matrix<Form, Type>&
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

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 Form, class Type>
inline Foam::Matrix<Form, 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 Form, class Type>
inline Foam::autoPtr<Foam::Matrix<Form, Type> > Foam::Matrix<Form, Type>::clone() const
{
return autoPtr<Matrix<T> >(new Matrix<T>(*this));
return autoPtr<Matrix<Form, Type> >(new Matrix<Form, Type>(*this));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return the number of rows
template<class T>
inline Foam::label Foam::Matrix<T>::n() const
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::n() const
{
return n_;
}
//- Return the number of columns
template<class T>
inline Foam::label Foam::Matrix<T>::m() const
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::m() const
{
return m_;
}
//- Return the number of columns
template<class T>
inline Foam::label Foam::Matrix<T>::size() const
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, 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 Form, class Type>
inline void Foam::Matrix<Form, Type>::checki(const label i) const
{
if (!n_)
{
FatalErrorIn("Matrix<T>::checki(const label)")
FatalErrorIn("Matrix<Form, Type>::checki(const label)")
<< "attempt to access element from zero sized row"
<< abort(FatalError);
}
else if (i<0 || i>=n_)
{
FatalErrorIn("Matrix<T>::checki(const label)")
FatalErrorIn("Matrix<Form, Type>::checki(const label)")
<< "index " << i << " out of range 0 ... " << n_-1
<< abort(FatalError);
}
}
// Check index j is within valid range (0 ... n-1).
template<class T>
inline void Foam::Matrix<T>::checkj(const label j) const
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checkj(const label j) const
{
if (!m_)
{
FatalErrorIn("Matrix<T>::checkj(const label)")
FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
<< "attempt to access element from zero sized column"
<< abort(FatalError);
}
else if (j<0 || j>=m_)
{
FatalErrorIn("Matrix<T>::checkj(const label)")
FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
<< "index " << j << " out of range 0 ... " << m_-1
<< abort(FatalError);
}
@ -108,9 +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 Form, class Type>
inline Type* Foam::Matrix<Form, 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 Form, class Type>
inline const Type* Foam::Matrix<Form, 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 Form, class Type>
Foam::Matrix<Form, 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 Form, class Type>
Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
{
// Anull matrix
M.clear();
is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&)");
token firstToken(is);
is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&) : reading first token");
if (firstToken.isLabel())
{
@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
label nm = M.n_*M.m_;
// Read list contents depending on data format
if (is.format() == IOstream::ASCII || !contiguous<T>())
if (is.format() == IOstream::ASCII || !contiguous<Type>())
{
// Read beginning of contents
char listDelimiter = is.readBeginList("Matrix");
@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
if (nm)
{
M.allocate();
T* v = M.v_[0];
Type* v = M.v_[0];
if (listDelimiter == token::BEGIN_LIST)
{
@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
is.fatalCheck
(
"operator>>(Istream&, Matrix<T>&) : "
"operator>>(Istream&, Matrix<Form, Type>&) : "
"reading entry"
);
}
@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
}
else
{
T element;
Type element;
is >> element;
is.fatalCheck
(
"operator>>(Istream&, Matrix<T>&) : "
"operator>>(Istream&, Matrix<Form, Type>&) : "
"reading the single entry"
);
@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
if (nm)
{
M.allocate();
T* v = M.v_[0];
Type* v = M.v_[0];
is.read(reinterpret_cast<char*>(v), nm*sizeof(T));
is.read(reinterpret_cast<char*>(v), nm*sizeof(Type));
is.fatalCheck
(
"operator>>(Istream&, Matrix<T>&) : "
"operator>>(Istream&, Matrix<Form, Type>&) : "
"reading the binary block"
);
}
@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
}
else
{
FatalIOErrorIn("operator>>(Istream&, Matrix<T>&)", is)
FatalIOErrorIn("operator>>(Istream&, Matrix<Form, Type>&)", is)
<< "incorrect first token, expected <int>, found "
<< firstToken.info()
<< exit(FatalIOError);
@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
}
template<class T>
Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
template<class Form, class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
{
label nm = M.n_*M.m_;
os << M.n() << token::SPACE << M.m();
// Write list contents depending on data format
if (os.format() == IOstream::ASCII || !contiguous<T>())
if (os.format() == IOstream::ASCII || !contiguous<Type>())
{
if (nm)
{
bool uniform = false;
const T* v = M.v_[0];
const Type* v = M.v_[0];
if (nm > 1 && contiguous<T>())
if (nm > 1 && contiguous<Type>())
{
uniform = true;
@ -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));
}
}

92
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H Normal file
View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::RectangularMatrix
Description
A templated 2D rectangular matrix of objects of \<T\>, where the n x n matrix
dimension is known and used for subscript bounds checking, etc.
SourceFiles
RectangularMatrixI.H
RectangularMatrix.C
\*---------------------------------------------------------------------------*/
#ifndef RectangularMatrix_H
#define RectangularMatrix_H
#include "Matrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class RectangularMatrix
:
public Matrix<RectangularMatrix<Type>, Type>
{
public:
// Constructors
//- Null constructor.
inline RectangularMatrix();
//- Construct given number of rows and columns,
inline RectangularMatrix(const label m, const label n);
//- Construct with given number of rows and columns
// and value for all elements.
inline RectangularMatrix(const label m, const label n, const Type&);
//- Construct from Istream.
inline RectangularMatrix(Istream&);
//- Clone
inline autoPtr<RectangularMatrix<Type> > clone() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
# include "RectangularMatrixI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

70
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H Normal file
View File

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix()
:
Matrix<RectangularMatrix<Type>, Type>()
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const label m,
const label n
)
:
Matrix<RectangularMatrix<Type>, Type>(m, n)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const label m,
const label n,
const Type& t
)
:
Matrix<RectangularMatrix<Type>, Type>(m, n, t)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix(Istream& is)
:
Matrix<RectangularMatrix<Type>, Type>(is)
{}
template<class Type>
inline Foam::autoPtr<Foam::RectangularMatrix<Type> >
Foam::RectangularMatrix<Type>::clone() const
{
return autoPtr<RectangularMatrix<Type> >(new RectangularMatrix<Type>(*this));
}
// ************************************************************************* //

81
src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H → src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
View File

@ -23,22 +23,22 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::scalarMatrix
Foam::SquareMatrix
Description
Foam::scalarMatrix
A templated 2D square matrix of objects of \<T\>, where the n x n matrix
dimension is known and used for subscript bounds checking, etc.
SourceFiles
scalarMatrix.C
SquareMatrixI.H
SquareMatrix.C
\*---------------------------------------------------------------------------*/
#ifndef scalarMatrix_H
#define scalarMatrix_H
#ifndef SquareMatrix_H
#define SquareMatrix_H
#include "Matrix.H"
#include "scalarField.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,65 +46,38 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class scalarMatrix Declaration
Class Matrix Declaration
\*---------------------------------------------------------------------------*/
class scalarMatrix
template<class Type>
class SquareMatrix
:
public Matrix<scalar>
public Matrix<SquareMatrix<Type>, Type>
{
public:
// Constructors
//- Construct null
scalarMatrix();
//- Null constructor.
inline SquareMatrix();
//- Construct given size
scalarMatrix(const label);
//- Construct given number of rows/columns.
inline SquareMatrix(const label n);
//- Construct from Matrix<scalar>
scalarMatrix(const Matrix<scalar>&);
//- Construct given number of rows and columns,
// It checks that m == n.
inline SquareMatrix(const label m, const label n);
//- Construct from Istream
scalarMatrix(Istream&);
//- Construct with given number of rows/columns
// and value for all elements.
inline SquareMatrix(const label n, const Type&);
//- Construct from Istream.
inline SquareMatrix(Istream&);
// Member Functions
//- 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);
//- 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;
//- LU decompose the matrix with pivoting
static void LUDecompose
(
Matrix<scalar>& matrix,
labelList& pivotIndices
);
//- LU back-substitution with given source, returning the solution
// in the source
template<class T>
static void LUBacksubstitute
(
const Matrix<scalar>& luMmatrix,
const labelList& pivotIndices,
Field<T>& 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);
//- Clone
inline autoPtr<SquareMatrix<Type> > clone() const;
};
@ -114,9 +87,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "scalarMatrixTemplates.C"
#endif
# include "SquareMatrixI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

75
src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H Normal file
View File

@ -0,0 +1,75 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix()
:
Matrix<SquareMatrix<Type>, Type>()
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
:
Matrix<SquareMatrix<Type>, Type>(n, n)
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(const label m, const label n)
:
Matrix<SquareMatrix<Type>, Type>(m, n)
{
if (m != n)
{
FatalErrorIn
(
"SquareMatrix<Type>::SquareMatrix(const label m, const label n)"
) << "m != n for constructing a square matrix" << exit(FatalError);
}
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(const label n, const Type& t)
:
Matrix<SquareMatrix<Type>, Type>(n, t)
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
:
Matrix<SquareMatrix<Type>, Type>(is)
{}
template<class Type>
inline Foam::autoPtr<Foam::SquareMatrix<Type> >
Foam::SquareMatrix<Type>::clone() const
{
return autoPtr<SquareMatrix<Type> >(new SquareMatrix<Type>(*this));
}
// ************************************************************************* //

403
src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C Normal file
View File

@ -0,0 +1,403 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "scalarMatrices.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * //
Foam::SVD::SVD(const scalarRectangularMatrix& 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"
"(scalarRectangularMatrix& A, const scalar minCondition)"
) << "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
/*scalarRectangularMatrix 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/scalarMatrices/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
Foam::SVD
Description
Singular value decomposition of a rectangular matrix.
SourceFiles
SVDI.H
SVD.C
\*---------------------------------------------------------------------------*/
#ifndef SVD_H
#define SVD_H
#include "scalarMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/*---------------------------------------------------------------------------*\
Class SVD Declaration
\*---------------------------------------------------------------------------*/
class SVD
{
// Private data
//- Rectangular matrix with the same dimensions as the input
scalarRectangularMatrix U_;
//- square matrix V
scalarRectangularMatrix V_;
//- The singular values
DiagonalMatrix<scalar> S_;
//- The matrix product V S^(-1) U^T
scalarRectangularMatrix 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 scalarRectangularMatrix& A, const scalar minCondition = 0);
// Access functions
//- Return U
inline const scalarRectangularMatrix& U() const;
//- Return the square matrix V
inline const scalarRectangularMatrix& V() const;
//- Return the singular values
inline const scalarDiagonalMatrix& S() const;
//- Return VSinvUt (the pseudo inverse)
inline const scalarRectangularMatrix& 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/scalarMatrices/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::scalarRectangularMatrix& Foam::SVD::U() const
{
return U_;
}
inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const
{
return V_;
}
inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const
{
return S_;
}
inline const Foam::scalarRectangularMatrix& 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;
}
// ************************************************************************* //

293
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C Normal file
View File

@ -0,0 +1,293 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "scalarMatrices.H"
#include "SVD.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::LUDecompose
(
scalarSquareMatrix& matrix,
labelList& pivotIndices
)
{
label n = matrix.n();
scalar vv[n];
for (register label i=0; i<n; i++)
{
scalar largestCoeff = 0.0;
scalar temp;
const scalar* __restrict__ matrixi = matrix[i];
for (register label j=0; j<n; j++)
{
if ((temp = mag(matrixi[j])) > largestCoeff)
{
largestCoeff = temp;
}
}
if (largestCoeff == 0.0)
{
FatalErrorIn
(
"LUdecompose"
"(scalarSquareMatrix& matrix, labelList& rowIndices)"
) << "Singular matrix" << exit(FatalError);
}
vv[i] = 1.0/largestCoeff;
}
for (register label j=0; j<n; j++)
{
scalar* __restrict__ matrixj = matrix[j];
for (register label i=0; i<j; i++)
{
scalar* __restrict__ matrixi = matrix[i];
scalar sum = matrixi[j];
for (register label k=0; k<i; k++)
{
sum -= matrixi[k]*matrix[k][j];
}
matrixi[j] = sum;
}
label iMax = 0;
scalar largestCoeff = 0.0;
for (register label i=j; i<n; i++)
{
scalar* __restrict__ matrixi = matrix[i];
scalar sum = matrixi[j];
for (register label k=0; k<j; k++)
{
sum -= matrixi[k]*matrix[k][j];
}
matrixi[j] = sum;
scalar temp;
if ((temp = vv[i]*mag(sum)) >= largestCoeff)
{
largestCoeff = temp;
iMax = i;
}
}
pivotIndices[j] = iMax;
if (j != iMax)
{
scalar* __restrict__ matrixiMax = matrix[iMax];
for (register label k=0; k<n; k++)
{
Swap(matrixj[k], matrixiMax[k]);
}
vv[iMax] = vv[j];
}
if (matrixj[j] == 0.0)
{
matrixj[j] = SMALL;
}
if (j != n-1)
{
scalar rDiag = 1.0/matrixj[j];
for (register label i=j+1; i<n; i++)
{
matrix[i][j] *= rDiag;
}
}
}
}
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::multiply
(
scalarRectangularMatrix& ans, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B
)
{
if (A.m() != B.n())
{
FatalErrorIn
(
"multiply("
"scalarRectangularMatrix& answer "
"const scalarRectangularMatrix& A, "
"const scalarRectangularMatrix& B)"
) << "A and B must have identical inner dimensions but A.m = "
<< A.m() << " and B.n = " << B.n()
<< abort(FatalError);
}
ans = scalarRectangularMatrix(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
(
scalarRectangularMatrix& ans, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B,
const scalarRectangularMatrix& C
)
{
if (A.m() != B.n())
{
FatalErrorIn
(
"multiply("
"const scalarRectangularMatrix& A, "
"const scalarRectangularMatrix& B, "
"const scalarRectangularMatrix& C, "
"scalarRectangularMatrix& 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 scalarRectangularMatrix& A, "
"const scalarRectangularMatrix& B, "
"const scalarRectangularMatrix& C, "
"scalarRectangularMatrix& answer)"
) << "B and C must have identical inner dimensions but B.m = "
<< B.m() << " and C.n = " << C.n()
<< abort(FatalError);
}
ans = scalarRectangularMatrix(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
(
scalarRectangularMatrix& ans, // value changed in return
const scalarRectangularMatrix& A,
const DiagonalMatrix<scalar>& B,
const scalarRectangularMatrix& C
)
{
if (A.m() != B.size())
{
FatalErrorIn
(
"multiply("
"const scalarRectangularMatrix& A, "
"const DiagonalMatrix<scalar>& B, "
"const scalarRectangularMatrix& C, "
"scalarRectangularMatrix& 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 scalarRectangularMatrix& A, "
"const DiagonalMatrix<scalar>& B, "
"const scalarRectangularMatrix& C, "
"scalarRectangularMatrix& answer)"
) << "B and C must have identical inner dimensions but B.m = "
<< B.size() << " and C.n = " << C.n()
<< abort(FatalError);
}
ans = scalarRectangularMatrix(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::RectangularMatrix<Foam::scalar> Foam::SVDinv
(
const scalarRectangularMatrix& A,
scalar minCondition
)
{
SVD svd(A, minCondition);
return svd.VSinvUt();
}
// ************************************************************************* //

137
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H Normal file
View File

@ -0,0 +1,137 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
scalarMatrices
Description
Scalar matrices
SourceFiles
scalarMatrices.C
scalarMatricesTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef scalarMatrices_H
#define scalarMatrices_H
#include "RectangularMatrix.H"
#include "SquareMatrix.H"
#include "DiagonalMatrix.H"
#include "scalarField.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef RectangularMatrix<scalar> scalarRectangularMatrix;
typedef SquareMatrix<scalar> scalarSquareMatrix;
typedef DiagonalMatrix<scalar> scalarDiagonalMatrix;
//- Solve the matrix using Gaussian elimination with pivoting,
// returning the solution in the source
template<class Type>
void solve(scalarSquareMatrix& matrix, Field<Type>& source);
//- Solve the matrix using Gaussian elimination with pivoting
// and return the solution
template<class Type>
void solve
(
Field<Type>& psi,
const scalarSquareMatrix& matrix,
const Field<Type>& source
);
//- LU decompose the matrix with pivoting
void LUDecompose
(
scalarSquareMatrix& matrix,
labelList& pivotIndices
);
//- LU back-substitution with given source, returning the solution
// in the source
template<class Type>
void LUBacksubstitute
(
const scalarSquareMatrix& luMmatrix,
const labelList& pivotIndices,
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 Type>
void LUsolve(scalarSquareMatrix& matrix, Field<Type>& source);
void multiply
(
scalarRectangularMatrix& answer, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B
);
void multiply
(
scalarRectangularMatrix& answer, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B,
const scalarRectangularMatrix& C
);
void multiply
(
scalarRectangularMatrix& answer, // value changed in return
const scalarRectangularMatrix& A,
const DiagonalMatrix<scalar>& B,
const scalarRectangularMatrix& C
);
//- Return the inverse of matrix A using SVD
scalarRectangularMatrix SVDinv
(
const scalarRectangularMatrix& A,
scalar minCondition = 0
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "scalarMatricesTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

49
src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C → src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
View File

@ -24,16 +24,16 @@ License
\*---------------------------------------------------------------------------*/
#include "scalarMatrix.H"
#include "scalarMatrices.H"
#include "Swap.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
void Foam::scalarMatrix::solve
template<class Type>
void Foam::solve
(
Matrix<scalar>& tmpMatrix,
Field<T>& sourceSol
scalarSquareMatrix& tmpMatrix,
Field<Type>& sourceSol
)
{
label n = tmpMatrix.n();
@ -68,7 +68,7 @@ void Foam::scalarMatrix::solve
// Check that the system of equations isn't singular
if (mag(tmpMatrix[i][i]) < 1e-20)
{
FatalErrorIn("scalarMatrix::solve()")
FatalErrorIn("solve(scalarSquareMatrix&, Field<Type>& sourceSol)")
<< "Singular Matrix"
<< exit(FatalError);
}
@ -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,21 +101,26 @@ void Foam::scalarMatrix::solve
}
template<class T>
void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
template<class Type>
void Foam::solve
(
Field<Type>& psi,
const scalarSquareMatrix& matrix,
const Field<Type>& source
)
{
Matrix<scalar> tmpMatrix = *this;
scalarSquareMatrix tmpMatrix = matrix;
psi = source;
solve(tmpMatrix, psi);
}
template<class T>
void Foam::scalarMatrix::LUBacksubstitute
template<class Type>
void Foam::LUBacksubstitute
(
const Matrix<scalar>& luMatrix,
const scalarSquareMatrix& luMatrix,
const labelList& pivotIndices,
Field<T>& sourceSol
Field<Type>& sourceSol
)
{
label n = luMatrix.n();
@ -125,7 +130,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 +141,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 +151,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 +164,11 @@ void Foam::scalarMatrix::LUBacksubstitute
}
template<class T>
void Foam::scalarMatrix::LUsolve
template<class Type>
void Foam::LUsolve
(
Matrix<scalar>& matrix,
Field<T>& sourceSol
scalarSquareMatrix& matrix,
Field<Type>& sourceSol
)
{
labelList pivotIndices(matrix.n());

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

@ -1,161 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "scalarMatrix.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::scalarMatrix::scalarMatrix()
{}
Foam::scalarMatrix::scalarMatrix(const label mSize)
:
Matrix<scalar>(mSize, mSize, 0.0)
{}
Foam::scalarMatrix::scalarMatrix(const Matrix<scalar>& matrix)
:
Matrix<scalar>(matrix)
{}
Foam::scalarMatrix::scalarMatrix(Istream& is)
:
Matrix<scalar>(is)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::scalarMatrix::LUDecompose
(
Matrix<scalar>& matrix,
labelList& pivotIndices
)
{
label n = matrix.n();
scalar vv[n];
for (register label i=0; i<n; i++)
{
scalar largestCoeff = 0.0;
scalar temp;
const scalar* __restrict__ matrixi = matrix[i];
for (register label j=0; j<n; j++)
{
if ((temp = mag(matrixi[j])) > largestCoeff)
{
largestCoeff = temp;
}
}
if (largestCoeff == 0.0)
{
FatalErrorIn
(
"scalarMatrix::LUdecompose"
"(Matrix<scalar>& matrix, labelList& rowIndices)"
) << "Singular matrix" << exit(FatalError);
}
vv[i] = 1.0/largestCoeff;
}
for (register label j=0; j<n; j++)
{
scalar* __restrict__ matrixj = matrix[j];
for (register label i=0; i<j; i++)
{
scalar* __restrict__ matrixi = matrix[i];
scalar sum = matrixi[j];
for (register label k=0; k<i; k++)
{
sum -= matrixi[k]*matrix[k][j];
}
matrixi[j] = sum;
}
label iMax = 0;
scalar largestCoeff = 0.0;
for (register label i=j; i<n; i++)
{
scalar* __restrict__ matrixi = matrix[i];
scalar sum = matrixi[j];
for (register label k=0; k<j; k++)
{
sum -= matrixi[k]*matrix[k][j];
}
matrixi[j] = sum;
scalar temp;
if ((temp = vv[i]*mag(sum)) >= largestCoeff)
{
largestCoeff = temp;
iMax = i;
}
}
pivotIndices[j] = iMax;
if (j != iMax)
{
scalar* __restrict__ matrixiMax = matrix[iMax];
for (register label k=0; k<n; k++)
{
Swap(matrixj[k], matrixiMax[k]);
}
vv[iMax] = vv[j];
}
if (matrixj[j] == 0.0)
{
matrixj[j] = SMALL;
}
if (j != n-1)
{
scalar rDiag = 1.0/matrixj[j];
for (register label i=j+1; i<n; i++)
{
matrix[i][j] *= rDiag;
}
}
}
}
// ************************************************************************* //

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

@ -28,55 +28,55 @@ 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)
scalarSquareMatrix(mSize),
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 scalarSquareMatrix& matrix,
const Field<Type>& source
)
:
scalarMatrix(matrix),
scalarSquareMatrix(matrix),
source_(source)
{}
template<class T>
Foam::simpleMatrix<T>::simpleMatrix(Istream& is)
template<class Type>
Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
:
scalarMatrix(is),
scalarSquareMatrix(is),
source_(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_;
scalarSquareMatrix tmpMatrix = *this;
Field<Type> sourceSol = source_;
scalarMatrix::solve(tmpMatrix, sourceSol);
Foam::solve(tmpMatrix, sourceSol);
return sourceSol;
}
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_;
scalarSquareMatrix luMatrix = *this;
Field<Type> sourceSol = source_;
scalarMatrix::LUsolve(luMatrix, sourceSol);
Foam::LUsolve(luMatrix, sourceSol);
return sourceSol;
}
@ -84,82 +84,82 @@ 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);
}
scalarMatrix::operator=(m);
scalarSquareMatrix::operator=(m);
source_ = m.source_;
}
// * * * * * * * * * * * * * * * 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),
static_cast<const scalarSquareMatrix&>(m1)
+ static_cast<const scalarSquareMatrix&>(m2),
m1.source_ + m2.source_
);
}
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),
static_cast<const scalarSquareMatrix&>(m1)
- static_cast<const scalarSquareMatrix&>(m2),
m1.source_ - m2.source_
);
}
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_;
os << static_cast<const scalarSquareMatrix&>(m) << nl << m.source_;
return os;
}

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

@ -36,7 +36,7 @@ SourceFiles
#ifndef simpleMatrix_H
#define simpleMatrix_H
#include "scalarMatrix.H"
#include "scalarMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,35 +45,14 @@ namespace Foam
// Forward declaration of friend functions and operators
template<class T>
template<class Type>
class simpleMatrix;
template<class T>
simpleMatrix<T> operator+
(
const simpleMatrix<T>&,
const simpleMatrix<T>&
);
template<class T>
simpleMatrix<T> operator-
(
const simpleMatrix<T>&,
const simpleMatrix<T>&
);
template<class T>
simpleMatrix<T> operator*
(
const scalar,
const simpleMatrix<T>&
);
template<class T>
template<class Type>
Ostream& operator<<
(
Ostream&,
const simpleMatrix<T>&
const simpleMatrix<Type>&
);
@ -81,14 +60,14 @@ Ostream& operator<<
Class simpleMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class T>
template<class Type>
class simpleMatrix
:
public scalarMatrix
public scalarSquareMatrix
{
// Private data
Field<T> source_;
Field<Type> source_;
public:
@ -99,25 +78,25 @@ public:
simpleMatrix(const label);
//- Construct from components
simpleMatrix(const scalarMatrix&, const Field<T>&);
simpleMatrix(const scalarSquareMatrix&, 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,49 +104,52 @@ 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>&);
// Friend Operators
friend simpleMatrix<T> operator+ <T>
(
const simpleMatrix<T>&,
const simpleMatrix<T>&
);
friend simpleMatrix<T> operator- <T>
(
const simpleMatrix<T>&,
const simpleMatrix<T>&
);
friend simpleMatrix<T> operator* <T>
(
const scalar,
const simpleMatrix<T>&
);
void operator=(const simpleMatrix<Type>&);
// Ostream Operator
friend Ostream& operator<< <T>
friend Ostream& operator<< <Type>
(
Ostream&,
const simpleMatrix<T>&
const simpleMatrix<Type>&
);
};
// Global operators
template<class Type>
simpleMatrix<Type> operator+
(
const simpleMatrix<Type>&,
const simpleMatrix<Type>&
);
template<class Type>
simpleMatrix<Type> operator-
(
const simpleMatrix<Type>&,
const simpleMatrix<Type>&
);
template<class Type>
simpleMatrix<Type> operator*
(
const scalar,
const simpleMatrix<Type>&
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

31
src/OpenFOAM/meshes/meshShapes/face/face.C
View File

@ -526,7 +526,7 @@ Foam::point Foam::face::centre(const pointField& meshPoints) const
}
Foam::vector Foam::face::normal(const pointField& meshPoints) const
Foam::vector Foam::face::normal(const pointField& p) const
{
// Calculate the normal by summing the face triangle normals.
// Changed to deal with small concavity by using a central decomposition
@ -539,38 +539,43 @@ Foam::vector Foam::face::normal(const pointField& meshPoints) const
{
return triPointRef
(
meshPoints[operator[](0)],
meshPoints[operator[](1)],
meshPoints[operator[](2)]
p[operator[](0)],
p[operator[](1)],
p[operator[](2)]
).normal();
}
vector n = vector::zero;
point centrePoint = Foam::average(points(meshPoints));
label nPoints = size();
point nextPoint = centrePoint;
register label pI;
point centrePoint = vector::zero;
for (pI = 0; pI < nPoints; pI++)
{
centrePoint += p[operator[](pI)];
}
centrePoint /= nPoints;
vector n = vector::zero;
point nextPoint = centrePoint;
for (pI = 0; pI < nPoints; pI++)
{
if (pI < nPoints - 1)
{
nextPoint = meshPoints[operator[](pI + 1)];
nextPoint = p[operator[](pI + 1)];
}
else
{
nextPoint = meshPoints[operator[](0)];
nextPoint = p[operator[](0)];
}
// Note: for best accuracy, centre point always comes last
//
n += triPointRef
(
meshPoints[operator[](pI)],
p[operator[](pI)],
nextPoint,
centrePoint
).normal();

23
src/OpenFOAM/primitives/Lists/stringListTemplates.C
View File

@ -25,9 +25,7 @@ License
\*---------------------------------------------------------------------------*/
#include "labelList.H"
#include <sys/types.h>
#include <regex.h>
#include "regularExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -41,24 +39,12 @@ labelList findStrings(const string& regexp, const StringList& sl)
{
labelList matches(sl.size());
regex_t *preg = new regex_t;
if (regcomp(preg, regexp.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
{
WarningIn("findStrings(const string& regexp, const stringList& sl)")
<< "Failed to compile regular expression " << regexp
<< endl;
return matches;
}
size_t nmatch = 0;
regmatch_t *pmatch = NULL;
regularExpression re(regexp);
label matchi = 0;
forAll(sl, i)
{
if (regexec(preg, sl[i].c_str(), nmatch, pmatch, 0) == 0)
if (re.matches(sl[i]))
{
matches[matchi++] = i;
}
@ -66,9 +52,6 @@ labelList findStrings(const string& regexp, const StringList& sl)
matches.setSize(matchi);
regfree(preg);
delete preg;
return matches;
}

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

11
src/meshTools/sets/topoSets/faceSet.C
View File

@ -205,11 +205,12 @@ void faceSet::sync(const polyMesh& mesh)
reduce(nAdded, sumOp<label>());
if (nAdded > 0)
{
Info<< "Added an additional " << nAdded << " faces on coupled patches. "
<< "(processorPolyPatch, cyclicPolyPatch)" << endl;
}
//if (nAdded > 0)
//{
// Info<< "Added an additional " << nAdded
// << " faces on coupled patches. "
// << "(processorPolyPatch, cyclicPolyPatch)" << endl;
//}
}

1
src/postProcessing/Allwmake
View File

@ -6,5 +6,6 @@ wmake libo postCalc
wmake libso forces
wmake libso fieldAverage
wmake libso foamCalcFunctions
wmake libso minMaxFields
# ----------------------------------------------------------------- end-of-file

38
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -114,7 +114,7 @@ Foam::scalarField Foam::chemistryModel::omega
forAll(reactions_, i)
{
const reaction& R = reactions_[i];
scalar omegai = omega
(
R, c, T, p, pf, cf, lRef, pr, cr, rRef
@ -164,13 +164,13 @@ Foam::scalar Foam::chemistryModel::omega
pf = 1.0;
pr = 1.0;
label Nl = R.lhs().size();
label Nr = R.rhs().size();
label slRef = 0;
lRef = R.lhs()[slRef].index;
pf = kf;
for(label s=1; s<Nl; s++)
{
@ -212,7 +212,7 @@ Foam::scalar Foam::chemistryModel::omega
label srRef = 0;
rRef = R.rhs()[srRef].index;
// find the matrix element and element position for the rhs
pr = kr;
for(label s=1; s<Nr; s++)
@ -250,7 +250,7 @@ Foam::scalar Foam::chemistryModel::omega
{
pr *= pow(cr, exp-1.0);
}
}
return pf*cf - pr*cr;
@ -313,12 +313,12 @@ void Foam::chemistryModel::jacobian
const scalar t,
const scalarField& c,
scalarField& dcdt,
Matrix<scalar>& dfdc
scalarSquareMatrix& dfdc
) const
{
scalar T = c[Ns_];
scalar p = c[Ns_ + 1];
scalarField c2(Ns(), 0.0);
for(label i=0; i<Ns(); i++)
{
@ -470,23 +470,23 @@ Foam::tmp<Foam::volScalarField> Foam::chemistryModel::tc() const
scalar pi = thermo_.p()[celli];
scalarField c(Ns_);
scalar cSum = 0.0;
for(label i=0; i<Ns_; i++)
{
scalar Yi = Y_[i][celli];
c[i] = rhoi*Yi/specieThermo_[i].W();
cSum += c[i];
}
forAll(reactions_, i)
{
const reaction& R = reactions_[i];
omega
(
R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef
);
forAll(R.rhs(), s)
{
scalar sr = R.rhs()[s].stoichCoeff;
@ -544,22 +544,22 @@ void Foam::chemistryModel::calculate()
{
RR_[i][celli] = 0.0;
}
scalar rhoi = rho_[celli];
scalar Ti = thermo_.T()[celli];
scalar pi = thermo_.p()[celli];
scalarField c(Ns_);
scalarField dcdt(nEqns(), 0.0);
for(label i=0; i<Ns_; i++)
{
scalar Yi = Y_[i][celli];
c[i] = rhoi*Yi/specieThermo_[i].W();
}
dcdt = omega(c, Ti, pi);
for(label i=0; i<Ns_; i++)
{
RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
@ -624,7 +624,7 @@ Foam::scalar Foam::chemistryModel::solve(const scalar t0, const scalar deltaT)
for(label i=0; i<Ns_; i++)
{
mixture += (c[i]/cTot)*specieThermo_[i];
}
}
Ti = mixture.TH(hi, Ti);
timeLeft -= dt;
@ -639,7 +639,7 @@ Foam::scalar Foam::chemistryModel::solve(const scalar t0, const scalar deltaT)
for(label i=0; i<Ns_; i++)
{
WTot += c[i]*specieThermo_[i].W();
}
}
WTot /= cTot;
for(label i=0; i<Ns_; i++)

3
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
View File

@ -39,7 +39,6 @@ SourceFiles
#include "hCombustionThermo.H"
#include "reactingMixture.H"
#include "Matrix.H"
#include "ODE.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -251,7 +250,7 @@ public:
const scalar t,
const scalarField& c,
scalarField& dcdt,
Matrix<scalar>& dfdc
scalarSquareMatrix& dfdc
) const;
//- Calculates the reaction rates

16
wmake/wmakeScheduler
View File

@ -137,6 +137,16 @@ done
# is returned and not of colouring pipe.
set -o pipefail
# Define function to colour output by argument 1
colourPipe(){
if [ "$1" ]; then
(while read line; do setterm -foreground $1; echo "$line" ; done; setterm -foreground default)
else
cat
fi
}
colourIndex=0
while :
@ -156,14 +166,14 @@ do
if lockfile -r0 "$lockFile" 2>/dev/null; then
if [ "$WM_COLOURS" ]; then
# Set colour
colourString=`setterm -foreground ${colours[$colourIndex]}`
colour="${colours[$colourIndex]}"
if [ "$host" = "$HOST" ]; then
eval $* 2>&1 | sed -e "s/^/$colourString/"
eval $* 2>&1 | colourPipe "$colour"
elif [ -n "$JOB_ID" ]; then
qrsh -inherit -v PWD $host "$rcmd"
else
ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1 | sed -e "s/^/$colourString/"
ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1 | colourPipe "$colour"
fi
retval=$?
else