src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixDecoupledHOps.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Vector-matrix multiplication operations for a block matrix
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "BlockLduMatrix.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 template<class Type>
  34 Foam::tmp<Foam::Field<Type> >
  35 Foam::BlockLduMatrix<Type>::decoupledH(const Field<Type>& x) const
  36 {
  37     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
  38     typedef typename TypeCoeffField::linearTypeField linearTypeField;
  39
  40     // Create result
  41     tmp<Field<Type> > tresult
  42     (
  43         new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
  44     );
  45     Field<Type>& result = tresult();
  46
  47     const unallocLabelList& u = lduAddr().upperAddr();
  48     const unallocLabelList& l = lduAddr().lowerAddr();
  49
  50     const TypeCoeffField& Upper = this->upper();
  51
  52     // Create multiplication function object
  53     typename BlockCoeff<Type>::multiply mult;
  54
  55     // Lower multiplication
  56
  57     if (symmetric())
  58     {
  59         if (Upper.activeType() == blockCoeffBase::SCALAR)
  60         {
  61             const scalarTypeField& activeUpper = Upper.asScalar();
  62
  63             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  64             {
  65                 result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
  66             }
  67         }
  68         else if (Upper.activeType() == blockCoeffBase::LINEAR)
  69         {
  70             const linearTypeField& activeUpper = Upper.asLinear();
  71
  72             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  73             {
  74                 result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
  75             }
  76         }
  77     }
  78     else // Asymmetric matrix
  79     {
  80         const TypeCoeffField& Lower = this->lower();
  81
  82         if (Lower.activeType() == blockCoeffBase::SCALAR)
  83         {
  84             const scalarTypeField& activeLower = Lower.asScalar();
  85
  86             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  87             {
  88                 result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
  89             }
  90         }
  91         else if (Lower.activeType() == blockCoeffBase::LINEAR)
  92         {
  93             const linearTypeField& activeLower = Lower.asLinear();
  94
  95             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  96             {
  97                 result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
  98             }
  99         }
 100     }
 101
 102
 103     // Upper multiplication
 104
 105     if (Upper.activeType() == blockCoeffBase::SCALAR)
 106     {
 107         const scalarTypeField& activeUpper = Upper.asScalar();
 108
 109         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 110         {
 111             result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
 112         }
 113     }
 114     else if (Upper.activeType() == blockCoeffBase::LINEAR)
 115     {
 116         const linearTypeField& activeUpper = Upper.asLinear();
 117
 118         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 119         {
 120             result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
 121         }
 122     }
 123
 124     return tresult;
 125 }
 126
 127
 128 template<class Type>
 129 Foam::tmp<Foam::Field<Type> >
 130 Foam::BlockLduMatrix<Type>::decoupledFaceH(const Field<Type>& x) const
 131 {
 132     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 133     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 134
 135     const unallocLabelList& u = lduAddr().upperAddr();
 136     const unallocLabelList& l = lduAddr().lowerAddr();
 137
 138     // Create result
 139     tmp<Field<Type> > tresult(new Field<Type>(u.size(), pTraits<Type>::zero));
 140     Field<Type>& result = tresult();
 141
 142     const TypeCoeffField& Upper = this->upper();
 143
 144     // Create multiplication function object
 145     typename BlockCoeff<Type>::multiply mult;
 146
 147     // Lower multiplication
 148
 149     if (symmetric())
 150     {
 151         if (Upper.activeType() == blockCoeffBase::SCALAR)
 152         {
 153             const scalarTypeField& activeUpper = Upper.asScalar();
 154
 155             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 156             {
 157                 // This can be optimised with a subtraction.
 158                 // Currently not done for clarity.  HJ, 31/Oct/2007
 159                 result[coeffI] =
 160                     mult(activeUpper[coeffI], x[u[coeffI]])
 161                   - mult(activeUpper[coeffI], x[l[coeffI]]);
 162             }
 163         }
 164         else if (Upper.activeType() == blockCoeffBase::LINEAR)
 165         {
 166             const linearTypeField& activeUpper = Upper.asLinear();
 167
 168             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 169             {
 170                 // This can be optimised with a subtraction.
 171                 // Currently not done for clarity.  HJ, 31/Oct/2007
 172                 result[coeffI] =
 173                     mult(activeUpper[coeffI], x[u[coeffI]])
 174                   - mult(activeUpper[coeffI], x[l[coeffI]]);
 175             }
 176         }
 177     }
 178     else // Asymmetric matrix
 179     {
 180         const TypeCoeffField& Lower = this->lower();
 181
 182         if (Lower.activeType() == blockCoeffBase::SCALAR)
 183         {
 184             const scalarTypeField& activeUpper = Upper.asScalar();
 185             const scalarTypeField& activeLower = Lower.asScalar();
 186
 187             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 188             {
 189                 result[coeffI] =
 190                     mult(activeUpper[coeffI], x[u[coeffI]])
 191                   - mult(activeLower[coeffI], x[l[coeffI]]);
 192             }
 193         }
 194         else if (Lower.activeType() == blockCoeffBase::LINEAR)
 195         {
 196             const linearTypeField& activeUpper = Upper.asLinear();
 197             const linearTypeField& activeLower = Lower.asLinear();
 198
 199             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 200             {
 201                 result[coeffI] =
 202                     mult(activeUpper[coeffI], x[u[coeffI]])
 203                   - mult(activeLower[coeffI], x[l[coeffI]]);
 204             }
 205         }
 206     }
 207
 208     return tresult;
 209 }
 210
 211
 212 // ************************************************************************* //