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