src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockConstraint/BlockConstraintTools.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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     A storage mechanism which allows setting of the fixed value and
  26     consequently recovering the equation for a single row of the matrix as
  27     well as b. The equation is taken out of the matrix using a
  28     variant of compact matrix storage mechanism.
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "error.H"
  36
  37 #include "BlockConstraint.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  45
  46 template<class Type>
  47 void BlockConstraint<Type>::setMatrix
  48 (
  49     const BlockLduMatrix<Type>& matrix,
  50     const TypeField& x,
  51     const TypeField& b
  52 )
  53 {
  54     if (matrixCoeffsSet_)
  55     {
  56         FatalErrorIn
  57         (
  58             "void BlockConstraint<Type>::setMatrix\n"
  59             "(\n"
  60             "    const BlockLduMatrix<Type>& matrix,\n"
  61             "    const TypeField& x,\n"
  62             "    const TypeField& b\n"
  63             ")"
  64         )   << "matrix coefficients already set"
  65             << abort(FatalError);
  66     }
  67
  68     matrixCoeffsSet_ = true;
  69
  70     if (matrix.thereIsDiag())
  71     {
  72         diagCoeff_ = matrix.diag().getCoeff(rowID_);
  73     }
  74
  75     b_ = b[rowID_];
  76
  77     const label startFaceOwn =
  78         matrix.lduAddr().ownerStartAddr()[rowID_];
  79     const label endFaceOwn =
  80         matrix.lduAddr().ownerStartAddr()[rowID_ + 1];
  81     const label ownSize = endFaceOwn - startFaceOwn;
  82
  83     const label startFaceNbr =
  84         matrix.lduAddr().losortStartAddr()[rowID_];
  85     const label endFaceNbr =
  86         matrix.lduAddr().losortStartAddr()[rowID_ + 1];
  87     const label nbrSize = endFaceNbr - startFaceNbr;
  88
  89     const unallocLabelList& losort = matrix.lduAddr().losortAddr();
  90
  91     // Create losort addressing
  92     labelList losortAddr(nbrSize);
  93
  94     forAll (losortAddr, laI)
  95     {
  96         losortAddr[laI] = losort[startFaceNbr + laI];
  97     }
  98
  99     if (matrix.thereIsUpper())
 100     {
 101         // Get the upper coefficients
 102
 103         const TypeCoeffField& matrixUpper = matrix.upper();
 104
 105         // owner side
 106         upperCoeffsOwnerPtr_ = new TypeCoeffField(ownSize);
 107         TypeCoeffField& uOwn = *upperCoeffsOwnerPtr_;
 108
 109         matrixUpper.getSubset(uOwn, startFaceOwn, ownSize);
 110
 111         // neighbour side
 112         upperCoeffsNeighbourPtr_ = new TypeCoeffField(nbrSize);
 113         TypeCoeffField& uNbr = *upperCoeffsNeighbourPtr_;
 114
 115         matrixUpper.getSubset(uNbr, losortAddr);
 116     }
 117
 118     if (matrix.thereIsLower())
 119     {
 120         // Get the lower coefficients
 121
 122         const TypeCoeffField& matrixLower = matrix.lower();
 123
 124         // owner side
 125         lowerCoeffsOwnerPtr_ = new TypeCoeffField(ownSize);
 126         TypeCoeffField& lOwn = *lowerCoeffsOwnerPtr_;
 127
 128         matrixLower.getSubset(lOwn, startFaceOwn, ownSize);
 129
 130         // neighbour side
 131         lowerCoeffsNeighbourPtr_ = new TypeCoeffField(nbrSize);
 132         TypeCoeffField& lNbr = *lowerCoeffsNeighbourPtr_;
 133
 134         matrixLower.getSubset(lNbr, losortAddr);
 135     }
 136 }
 137
 138
 139 template<class Type>
 140 void BlockConstraint<Type>::eliminateEquation
 141 (
 142     BlockLduMatrix<Type>& matrix,
 143     TypeField& b
 144 ) const
 145 {
 146     const label startFaceOwn =
 147         matrix.lduAddr().ownerStartAddr()[rowID_];
 148     const label endFaceOwn =
 149         matrix.lduAddr().ownerStartAddr()[rowID_ + 1];
 150     const label ownSize = endFaceOwn - startFaceOwn;
 151
 152     const label startFaceNbr =
 153         matrix.lduAddr().losortStartAddr()[rowID_];
 154     const label endFaceNbr =
 155         matrix.lduAddr().losortStartAddr()[rowID_ + 1];
 156     const label nbrSize = endFaceNbr - startFaceNbr;
 157
 158     const unallocLabelList& owner = matrix.lduAddr().lowerAddr();
 159     const unallocLabelList& neighbour = matrix.lduAddr().upperAddr();
 160     const unallocLabelList& losort = matrix.lduAddr().losortAddr();
 161
 162     // Create losort addressing
 163     labelList losortAddr(nbrSize);
 164
 165     forAll (losortAddr, laI)
 166     {
 167         losortAddr[laI] = losort[startFaceNbr + laI];
 168     }
 169
 170     typename BlockCoeff<Type>::multiply mult;
 171
 172     // My index = rowID_
 173     if (matrix.symmetric())
 174     {
 175         TypeField bOwn;
 176         TypeField bNbr;
 177
 178         TypeCoeffField& upperLower = matrix.upper();
 179         upperLower.zeroOutSubset(startFaceOwn, ownSize);
 180         upperLower.zeroOutSubset(losortAddr);
 181
 182         bOwn = upperCoeffsOwner()*value_;
 183         bNbr = upperCoeffsNeighbour()*value_;
 184
 185         // owner side
 186
 187         forAll (bOwn, soI)
 188         {
 189             // add contribution to b of the neighbour (I am the owner)
 190             b[neighbour[startFaceOwn + soI]] -= bOwn[soI];
 191         }
 192
 193         // neighbour side
 194
 195         forAll (bNbr, snI)
 196         {
 197             // add contribution to b of owner (I am the neighbour)
 198             b[owner[losort[startFaceNbr + snI]]] -= bNbr[snI];
 199         }
 200     }
 201     else if (matrix.asymmetric())
 202     {
 203         // Do upper
 204         TypeCoeffField& matrixUpper = matrix.upper();
 205
 206         matrixUpper.zeroOutSubset(startFaceOwn, ownSize);
 207         TypeField bOwn = lowerCoeffsOwner()*value_;
 208
 209         // Do lower
 210         TypeCoeffField& matrixLower = matrix.lower();
 211
 212         matrixLower.zeroOutSubset(losortAddr);
 213         TypeField bNbr = upperCoeffsNeighbour()*value_;
 214
 215         // owner side
 216
 217         forAll (bOwn, soI)
 218         {
 219             // add contribution to b of the neighbour (I am the owner)
 220             b[neighbour[startFaceOwn + soI]] -= bOwn[soI];
 221         }
 222
 223         // neighbour side
 224
 225         forAll (bNbr, snI)
 226         {
 227             // add contribution to b of owner (I am the neighbour)
 228             b[owner[losort[startFaceNbr + snI]]] -= bNbr[snI];
 229         }
 230     }
 231 }
 232
 233
 234 template<class Type>
 235 void BlockConstraint<Type>::setSourceDiag
 236 (
 237     BlockLduMatrix<Type>& matrix,
 238     Field<Type>& x,
 239     Field<Type>& b
 240 ) const
 241 {
 242     const Type& fc = fixedComponents();
 243
 244     typename BlockCoeff<Type>::multiply mult;
 245
 246     if (mag(fc) > SMALL)
 247     {
 248         b[rowID()] =
 249             cmptMultiply
 250             (
 251                 fc,
 252                 mult(matrix.diag().getCoeff(rowID()), value())
 253             );
 254
 255         // set the solution to the right value as well
 256         x[rowID()] = cmptMultiply(fc, value());
 257     }
 258 }
 259
 260
 261 template<class Type>
 262 void BlockConstraint<Type>::reconstructMatrix
 263 (
 264     BlockLduMatrix<Type>& matrix
 265 ) const
 266 {
 267     if (!matrixCoeffsSet_)
 268     {
 269         FatalErrorIn
 270         (
 271             "void BlockConstraint<Type>::reconstructMatrix("
 272             "BlockLduMatrix<Type>& matrix)"
 273         )   << "matrix coefficients not set"
 274             << abort(FatalError);
 275     }
 276
 277     if (matrix.thereIsDiag())
 278     {
 279         matrix.diag().setCoeff(rowID_, diagCoeff_);
 280     }
 281
 282     const label startFaceOwn =
 283         matrix.lduAddr().ownerStartAddr()[rowID_];
 284     const label endFaceOwn =
 285         matrix.lduAddr().ownerStartAddr()[rowID_ + 1];
 286     const label ownSize = endFaceOwn - startFaceOwn;
 287
 288     const label startFaceNbr =
 289         matrix.lduAddr().losortStartAddr()[rowID_];
 290     const label endFaceNbr =
 291         matrix.lduAddr().losortStartAddr()[rowID_ + 1];
 292     const label nbrSize = endFaceNbr - startFaceNbr;
 293
 294     const unallocLabelList& losort = matrix.lduAddr().losortAddr();
 295
 296     // Create losort addressing
 297     labelList losortAddr(nbrSize);
 298
 299     forAll (losortAddr, laI)
 300     {
 301         losortAddr[laI] = losort[startFaceNbr + laI];
 302     }
 303
 304     if (matrix.thereIsUpper())
 305     {
 306         // Get the upper coefficients
 307
 308         TypeCoeffField& matrixUpper = matrix.upper();
 309
 310         // owner side
 311         const TypeCoeffField& uOwn = upperCoeffsOwner();
 312
 313         matrixUpper.setSubset(uOwn, startFaceOwn, ownSize);
 314
 315         // neighbour side
 316         const TypeCoeffField& uNbr = upperCoeffsNeighbour();
 317
 318         matrixUpper.setSubset(uNbr, losortAddr);
 319     }
 320
 321     if (matrix.thereIsLower())
 322     {
 323         // Get the lower coefficients
 324
 325         TypeCoeffField& matrixLower = matrix.lower();
 326
 327         // owner side
 328         const TypeCoeffField& lOwn = lowerCoeffsOwner();
 329
 330         matrixLower.setSubset(lOwn, startFaceOwn, ownSize);
 331
 332         // neighbour side
 333         const TypeCoeffField& lNbr = lowerCoeffsNeighbour();
 334
 335         matrixLower.setSubset(lNbr, losortAddr);
 336     }
 337 }
 338
 339
 340 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 341
 342 } // End namespace Foam
 343
 344 // ************************************************************************* //