Skip to content
GitLab
Commit 7d56dcf4 authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/hunt2/OpenFOAM/OpenFOAM-dev

parents 0dc639e9 df283b54
Hide whitespace changes
Inline Side-by-side
src/OpenFOAM/Make/files
View file @ 7d56dcf4
......@@ -179,6 +179,7 @@ matrices/solution/solution.C
scalarMatrix = matrices/scalarMatrix
$(scalarMatrix)/scalarMatrix.C
$(scalarMatrix)/SVD/SVD.C
LUscalarMatrix = matrices/LUscalarMatrix
$(LUscalarMatrix)/LUscalarMatrix.C
......
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C 0 → 100644
View file @ 7d56dcf4
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "DiagonalMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Type>& a)
:
List<Type>(min(a.n(), a.m()))
{
forAll(*this, i)
{
this->operator[](i) = a[i][i];
}
}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size)
:
List<Type>(size)
{}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size, const Type& val)
:
List<Type>(size, val)
{}
template<class Type>
Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::invert()
{
forAll(*this, i)
{
Type x = this->operator[](i);
if (mag(x) < VSMALL)
{
this->operator[](i) = Type(0);
}
else
{
this->operator[](i) = Type(1)/x;
}
}
return this;
}
template<class Type>
Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& A)
{
DiagonalMatrix<Type> Ainv = A;
forAll(A, i)
{
Type x = A[i];
if (mag(x) < VSMALL)
{
Ainv[i] = Type(0);
}
else
{
Ainv[i] = Type(1)/x;
}
}
return Ainv;
}
// ************************************************************************* //
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H 0 → 100644
View file @ 7d56dcf4
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::DiagonalMatrix<Type>
Description
DiagonalMatrix<Type> is a 2D diagonal matrix of objects
of type Type, size nxn
SourceFiles
DiagonalMatrix.C
\*---------------------------------------------------------------------------*/
#ifndef DiagonalMatrix_H
#define DiagonalMatrix_H
#include "List.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * //
template<class Type> class Matrix;
/*---------------------------------------------------------------------------*\
Class DiagonalMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class DiagonalMatrix
:
public List<Type>
{
public:
// Constructors
//- Construct from diagonal component of a Matrix
DiagonalMatrix<Type>(const Matrix<Type>&);
//- Construct empty from size
DiagonalMatrix<Type>(const label size);
//- Construct from size and a value
DiagonalMatrix<Type>(const label, const Type&);
// Member functions
//- Invert the diaganol matrix and return itself
DiagonalMatrix<Type>& invert();
};
// Global functions
//- Return the diagonal Matrix inverse
template<class Type>
DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>&);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "DiagonalMatrix.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
View file @ 7d56dcf4
......@@ -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)
);
}
}
......
src/OpenFOAM/containers/Matrix/Matrix.C → src/OpenFOAM/matrices/Matrix/Matrix.C
View file @ 7d56dcf4
......@@ -28,13 +28,13 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class T>
void Foam::Matrix<T>::allocate()
template<class Type>
void Foam::Matrix<Type>::allocate()
{
if (n_ && m_)
{
v_ = new T*[n_];
v_[0] = new T[n_*m_];
v_ = new Type*[n_];
v_[0] = new Type[n_*m_];
for (register label i=1; i<n_; i++)
{
......@@ -46,12 +46,12 @@ void Foam::Matrix<T>::allocate()
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
template<class T>
Foam::Matrix<T>::~Matrix()
template<class Type>
Foam::Matrix<Type>::~Matrix()
{
if (v_)
{
delete[] (v_[0]);
delete[] (v_[0]);
delete[] v_;
}
}
......@@ -59,16 +59,16 @@ Foam::Matrix<T>::~Matrix()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
const Foam::Matrix<T>& Foam::Matrix<T>::null()
template<class Type>
const Foam::Matrix<Type>& Foam::Matrix<Type>::null()
{
Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
Matrix<Type>* nullPtr = reinterpret_cast<Matrix<Type>*>(NULL);
return *nullPtr;
}
template<class T>
Foam::Matrix<T>::Matrix(const label n, const label m)
template<class Type>
Foam::Matrix<Type>::Matrix(const label n, const label m)
:
n_(n),
m_(m),
......@@ -76,7 +76,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
{
if (n_ < 0 || m_ < 0)
{
FatalErrorIn("Matrix<T>::Matrix(const label n, const label m)")
FatalErrorIn("Matrix<Type>::Matrix(const label n, const label m)")
<< "bad n, m " << n_ << ", " << m_
<< abort(FatalError);
}
......@@ -85,8 +85,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
}
template<class T>
Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
template<class Type>
Foam::Matrix<Type>::Matrix(const label n, const label m, const Type& a)
:
n_(n),
m_(m),
......@@ -96,7 +96,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
{
FatalErrorIn
(
"Matrix<T>::Matrix(const label n, const label m, const T&)"
"Matrix<Type>::Matrix(const label n, const label m, const T&)"
) << "bad n, m " << n_ << ", " << m_
<< abort(FatalError);
}
......@@ -105,7 +105,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
if (v_)
{
T* v = v_[0];
Type* v = v_[0];
label nm = n_*m_;
......@@ -117,8 +117,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
}
template<class T>
Foam::Matrix<T>::Matrix(const Matrix<T>& a)
template<class Type>
Foam::Matrix<Type>::Matrix(const Matrix<Type>& a)
:
n_(a.n_),
m_(a.m_),
......@@ -127,8 +127,8 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
if (a.v_)
{
allocate();
T* v = v_[0];
const T* av = a.v_[0];
Type* v = v_[0];
const Type* av = a.v_[0];
label nm = n_*m_;
for (register label i=0; i<nm; i++)
......@@ -138,13 +138,13 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
}
}
template<class T>
void Foam::Matrix<T>::clear()
template<class Type>
void Foam::Matrix<Type>::clear()
{
if (v_)
{
delete[] (v_[0]);
delete[] (v_[0]);
delete[] v_;
}
n_ = 0;
......@@ -153,8 +153,8 @@ void Foam::Matrix<T>::clear()
}
template<class T>
void Foam::Matrix<T>::transfer(Matrix<T>& a)
template<class Type>
void Foam::Matrix<Type>::transfer(Matrix<Type>& a)
{
clear();
......@@ -171,12 +171,30 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class T>
void Foam::Matrix<T>::operator=(const T& t)
template<class Type>
Foam::Matrix<Type> Foam::Matrix<Type>::T() const
{
const Matrix<Type>& A = *this;
Matrix<Type> At(m(), n());
for (register label i=0; i<n(); i++)
{
for (register label j=0; j<m(); j++)
{
At[j][i] = A[i][j];
}
}
return At;
}
template<class Type>
void Foam::Matrix<Type>::operator=(const Type& t)
{
if (v_)
{
T* v = v_[0];
Type* v = v_[0];
label nm = n_*m_;
for (register label i=0; i<nm; i++)
......@@ -188,12 +206,12 @@ void Foam::Matrix<T>::operator=(const T& t)
// Assignment operator. Takes linear time.
template<class T>
void Foam::Matrix<T>::operator=(const Matrix<T>& a)
template<class Type>
void Foam::Matrix<Type>::operator=(const Matrix<Type>& a)
{
if (this == &a)
{
FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
FatalErrorIn("Matrix<Type>::operator=(const Matrix<Type>&)")
<< "attempted assignment to self"
<< abort(FatalError);
}
......@@ -204,12 +222,12 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
n_ = a.n_;
m_ = a.m_;
allocate();
}
}
if (v_)
{
T* v = v_[0];
const T* av = a.v_[0];
Type* v = v_[0];
const Type* av = a.v_[0];
label nm = n_*m_;
for (register label i=0; i<nm; i++)
......@@ -222,15 +240,15 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
template<class T>
const T& Foam::max(const Matrix<T>& a)
template<class Type>
const Type& Foam::max(const Matrix<Type>& a)
{
label nm = a.n_*a.m_;
if (nm)
{
label curMaxI = 0;
const T* v = a.v_[0];
const Type* v = a.v_[0];
for (register label i=1; i<nm; i++)
{
......@@ -244,7 +262,7 @@ const T& Foam::max(const Matrix<T>& a)
}
else
{
FatalErrorIn("max(const Matrix<T>&)")
FatalErrorIn("max(const Matrix<Type>&)")
<< "matrix is empty"
<< abort(FatalError);
......@@ -254,15 +272,15 @@ const T& Foam::max(const Matrix<T>& a)
}
template<class T>
const T& Foam::min(const Matrix<T>& a)
template<class Type>
const Type& Foam::min(const Matrix<Type>& a)
{
label nm = a.n_*a.m_;
if (nm)
{
label curMinI = 0;
const T* v = a.v_[0];
const Type* v = a.v_[0];
for (register label i=1; i<nm; i++)
{
......@@ -276,7 +294,7 @@ const T& Foam::min(const Matrix<T>& a)
}
else
{
FatalErrorIn("min(const Matrix<T>&)")
FatalErrorIn("min(const Matrix<Type>&)")
<< "matrix is empty"
<< abort(FatalError);
......@@ -288,15 +306,15 @@ const T& Foam::min(const Matrix<T>& a)
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class T>
Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
template<class Type>
Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a)
{
Matrix<T> na(a.n_, a.m_);
Matrix<Type> na(a.n_, a.m_);
if (a.n_ && a.m_)
{
T* nav = na.v_[0];
const T* av = a.v_[0];
Type* nav = na.v_[0];
const Type* av = a.v_[0];
label nm = a.n_*a.m_;
for (register label i=0; i<nm; i++)
......@@ -309,30 +327,34 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
}
template<class T>
Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
template<class Type>
Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<Type>& b)
{
if (a.n_ != b.n_)
{
FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of rows: "
FatalErrorIn
(
"Matrix<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of columns: "
FatalErrorIn
(
"Matrix<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_
<< abort(FatalError);
}
Matrix<T> ab(a.n_, a.m_);
Matrix<Type> ab(a.n_, a.m_);
T* abv = ab.v_[0];
const T* av = a.v_[0];
const T* bv = b.v_[0];
Type* abv = ab.v_[0];
const Type* av = a.v_[0];
const Type* bv = b.v_[0];
label nm = a.n_*a.m_;;
for (register label i=0; i<nm; i++)
......@@ -344,30 +366,34 @@ Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
}
template<class T>
Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
template<class Type>
Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<Type>& b)
{
if (a.n_ != b.n_)
{
FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of rows: "
FatalErrorIn
(
"Matrix<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_
<< abort(FatalError);
}
if (a.m_ != b.m_)
{
FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
<< "attempted add matrices with different number of columns: "
FatalErrorIn
(
"Matrix<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
) << "attempted add matrices with different number of columns: "