FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

lduMatrixTemplates.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025     lduMatrix member H operations.
00026 
00027 \*---------------------------------------------------------------------------*/
00028 
00029 #include <OpenFOAM/lduMatrix.H>
00030 
00031 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00032 
00033 template<class Type>
00034 Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const
00035 {
00036     tmp<Field<Type> > tHpsi
00037     (
00038         new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
00039     );
00040 
00041     if (lowerPtr_ || upperPtr_)
00042     {
00043         Field<Type> & Hpsi = tHpsi();
00044 
00045         Type* __restrict__ HpsiPtr = Hpsi.begin();
00046 
00047         const Type* __restrict__ psiPtr = psi.begin();
00048 
00049         const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
00050         const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
00051 
00052         const scalar* __restrict__ lowerPtr = lower().begin();
00053         const scalar* __restrict__ upperPtr = upper().begin();
00054 
00055         register const label nFaces = upper().size();
00056 
00057         for (register label face=0; face<nFaces; face++)
00058         {
00059             HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
00060             HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
00061         }
00062     }
00063 
00064     return tHpsi;
00065 }
00066 
00067 template<class Type>
00068 Foam::tmp<Foam::Field<Type> >
00069 Foam::lduMatrix::H(const tmp<Field<Type> >& tpsi) const
00070 {
00071     tmp<Field<Type> > tHpsi(H(tpsi()));
00072     tpsi.clear();
00073     return tHpsi;
00074 }
00075 
00076 
00077 template<class Type>
00078 Foam::tmp<Foam::Field<Type> >
00079 Foam::lduMatrix::faceH(const Field<Type>& psi) const
00080 {
00081     if (lowerPtr_ || upperPtr_)
00082     {
00083         const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
00084         const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
00085 
00086         const unallocLabelList& l = lduAddr().lowerAddr();
00087         const unallocLabelList& u = lduAddr().upperAddr();
00088 
00089         tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
00090         Field<Type> & faceHpsi = tfaceHpsi();
00091 
00092         for (register label face=0; face<l.size(); face++)
00093         {
00094             faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
00095         }
00096 
00097         return tfaceHpsi;
00098     }
00099     else
00100     {
00101         FatalErrorIn("lduMatrix::faceH(const Field<Type>& psi) const")
00102             << "Cannot calculate faceH"
00103                " the matrix does not have any off-diagonal coefficients."
00104             << exit(FatalError);
00105 
00106         return tmp<Field<Type> >(NULL);
00107     }
00108 }
00109 
00110 
00111 template<class Type>
00112 Foam::tmp<Foam::Field<Type> >
00113 Foam::lduMatrix::faceH(const tmp<Field<Type> >& tpsi) const
00114 {
00115     tmp<Field<Type> > tfaceHpsi(faceH(tpsi()));
00116     tpsi.clear();
00117     return tfaceHpsi;
00118 }
00119 
00120 
00121 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines