src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     lduMatrix member H operations.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "lduMatrix.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 template<class Type>
  35 Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const
  36 {
  37     tmp<Field<Type> > tHpsi
  38     (
  39         new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
  40     );
  41
  42     if (lowerPtr_ || upperPtr_)
  43     {
  44         Field<Type> & Hpsi = tHpsi();
  45
  46         Type* __restrict__ HpsiPtr = Hpsi.begin();
  47
  48         const Type* __restrict__ psiPtr = psi.begin();
  49
  50         const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
  51         const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
  52
  53         const scalar* __restrict__ lowerPtr = lower().begin();
  54         const scalar* __restrict__ upperPtr = upper().begin();
  55
  56         register const label nFaces = upper().size();
  57
  58         for (register label face=0; face<nFaces; face++)
  59         {
  60             HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
  61             HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
  62         }
  63     }
  64
  65     return tHpsi;
  66 }
  67
  68 template<class Type>
  69 Foam::tmp<Foam::Field<Type> >
  70 Foam::lduMatrix::H(const tmp<Field<Type> >& tpsi) const
  71 {
  72     tmp<Field<Type> > tHpsi(H(tpsi()));
  73     tpsi.clear();
  74     return tHpsi;
  75 }
  76
  77
  78 template<class Type>
  79 Foam::tmp<Foam::Field<Type> >
  80 Foam::lduMatrix::faceH(const Field<Type>& psi) const
  81 {
  82     if (lowerPtr_ || upperPtr_)
  83     {
  84         const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
  85         const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
  86
  87         const unallocLabelList& l = lduAddr().lowerAddr();
  88         const unallocLabelList& u = lduAddr().upperAddr();
  89
  90         tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
  91         Field<Type> & faceHpsi = tfaceHpsi();
  92
  93         for (register label face=0; face<l.size(); face++)
  94         {
  95             faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
  96         }
  97
  98         return tfaceHpsi;
  99     }
 100     else
 101     {
 102         FatalErrorIn("lduMatrix::faceH(const Field<Type>& psi) const")
 103             << "Cannot calculate faceH"
 104                " the matrix does not have any off-diagonal coefficients."
 105             << exit(FatalError);
 106
 107         return tmp<Field<Type> >(NULL);
 108     }
 109 }
 110
 111
 112 template<class Type>
 113 Foam::tmp<Foam::Field<Type> >
 114 Foam::lduMatrix::faceH(const tmp<Field<Type> >& tpsi) const
 115 {
 116     tmp<Field<Type> > tfaceHpsi(faceH(tpsi()));
 117     tpsi.clear();
 118     return tfaceHpsi;
 119 }
 120
 121
 122 // ************************************************************************* //