src/OpenFOAM/matrices/constraint/constraintTools.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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     A storage mechanism which allows setting of the fixed value and
  27     consequently recovering the equation for a single row of the matrix as
  28     well as the source. The equation is taken out of the matrix using a
  29     variant of compact matrix storage mechanism.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "constraint.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  41
  42 template<class Type>
  43 template<template<class> class Matrix>
  44 void constraint<Type>::setMatrix
  45 (
  46     const Matrix<Type>& matrix
  47 )
  48 {
  49     if (matrixCoeffsSet_)
  50     {
  51         FatalErrorIn
  52         (
  53             "const scalarField& constraint<Type>::setMatrix"
  54         )   << "(const Matrix<Type>& matrix)"
  55             << "matrix coefficients already set"
  56             << abort(FatalError);
  57     }
  58
  59     matrixCoeffsSet_ = true;
  60
  61     if (matrix.hasDiag())
  62     {
  63         diagCoeff_ = matrix.diag()[rowID_];
  64     }
  65
  66     source_ = matrix.source()[rowID_];
  67
  68     const label startFaceOwn =
  69         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
  70     const label endFaceOwn =
  71         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_ + 1];
  72
  73     const label startFaceNbr =
  74         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
  75     const label endFaceNbr =
  76         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_ + 1];
  77
  78     const unallocLabelList& losort = matrix.psi().mesh().lduAddr().losortAddr();
  79
  80     if (matrix.hasUpper())
  81     {
  82         // get the upper coefficients
  83
  84         const scalarField& matrixUpper = matrix.upper();
  85
  86         // owner side
  87         upperCoeffsOwnerPtr_ = new scalarField(endFaceOwn - startFaceOwn);
  88
  89         scalarField& uOwn = *upperCoeffsOwnerPtr_;
  90
  91         label faceIndex = startFaceOwn;
  92
  93         forAll(uOwn, uOwnI)
  94         {
  95             uOwn[uOwnI] = matrixUpper[faceIndex];
  96
  97             faceIndex++;
  98         }
  99
 100         // neighbour side
 101         upperCoeffsNeighbourPtr_ = new scalarField(endFaceNbr - startFaceNbr);
 102
 103         scalarField& uNbr = *upperCoeffsNeighbourPtr_;
 104
 105         faceIndex = startFaceNbr;
 106
 107         forAll (uNbr, uNbrI)
 108         {
 109             uNbr[uNbrI] = matrixUpper[losort[faceIndex]];
 110
 111             faceIndex++;
 112         }
 113     }
 114
 115     if (matrix.hasLower())
 116     {
 117         // get the lower coefficients
 118
 119         const scalarField& matrixLower = matrix.lower();
 120
 121         // owner side
 122         lowerCoeffsOwnerPtr_ = new scalarField(endFaceOwn - startFaceOwn);
 123
 124         scalarField& lOwn = *lowerCoeffsOwnerPtr_;
 125
 126         label faceIndex = startFaceOwn;
 127
 128         forAll (lOwn, lOwnI)
 129         {
 130             lOwn[lOwnI] = matrixLower[faceIndex];
 131
 132             faceIndex++;
 133         }
 134
 135         // neighbour side
 136         lowerCoeffsNeighbourPtr_ = new scalarField(endFaceNbr - startFaceNbr);
 137
 138         scalarField& lNbr = *lowerCoeffsNeighbourPtr_;
 139
 140         faceIndex = startFaceNbr;
 141
 142         forAll (lNbr, lNbrI)
 143         {
 144             lNbr[lNbrI] = matrixLower[losort[faceIndex]];
 145
 146             faceIndex++;
 147         }
 148     }
 149 }
 150
 151
 152 template<class Type>
 153 template<template<class> class Matrix>
 154 void constraint<Type>::eliminateEquation
 155 (
 156     Matrix<Type>& matrix,
 157     const label rowID,
 158     const Type& value
 159 )
 160 {
 161     Field<Type>& source = matrix.source();
 162
 163     const label startFaceOwn =
 164         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID];
 165     const label endFaceOwn =
 166         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID + 1];
 167
 168     const label startFaceNbr =
 169         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID];
 170     const label endFaceNbr =
 171         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID + 1];
 172
 173     const unallocLabelList& owner = matrix.psi().mesh().lduAddr().lowerAddr();
 174     const unallocLabelList& neighbour =
 175         matrix.psi().mesh().lduAddr().upperAddr();
 176     const unallocLabelList& losort = matrix.psi().mesh().lduAddr().losortAddr();
 177
 178     // My index =  rowID
 179     if (matrix.symmetric())
 180     {
 181         // get the coefficients
 182         scalarField* coeffs = NULL;
 183
 184         if (matrix.hasUpper())
 185         {
 186             coeffs = &matrix.upper();
 187         }
 188         else if (matrix.hasLower())
 189         {
 190             coeffs = &matrix.lower();
 191         }
 192
 193         scalarField& matrixCoeffs = *coeffs;
 194
 195         // owner side
 196         label faceIndex = startFaceOwn;
 197
 198         while (faceIndex < endFaceOwn)
 199         {
 200             // add contribution to source of the neighbour (I am the owner)
 201             source[neighbour[faceIndex]] -= matrixCoeffs[faceIndex]*value;
 202
 203             // I am the owner =  owner[faceIndex]
 204             // set off-diagonal coefficient to zero
 205             matrixCoeffs[faceIndex] = 0.0;
 206
 207             faceIndex++;
 208         }
 209
 210         // neighbour side
 211         faceIndex = startFaceNbr;
 212
 213         while (faceIndex < endFaceNbr)
 214         {
 215             // add contribution to source of owner (I am the neighbour)
 216             source[owner[losort[faceIndex]]] -=
 217                 matrixCoeffs[losort[faceIndex]]*value;
 218
 219             // I am the neighbour = neighbour[losort[faceIndex]]
 220             // set off-diagonal coefficient to zero
 221             matrixCoeffs[losort[faceIndex]] = 0.0;
 222
 223             faceIndex++;
 224         }
 225     }
 226     else if (matrix.asymmetric())
 227     {
 228         scalarField& matrixUpper = matrix.upper();
 229
 230         scalarField& matrixLower = matrix.lower();
 231
 232         // owner side
 233         label faceIndex = startFaceOwn;
 234
 235         while (faceIndex < endFaceOwn)
 236         {
 237             // add contribution to source of the neighbour (I am the owner)
 238             source[neighbour[faceIndex]] -= matrixLower[faceIndex]*value;
 239
 240             // set off-diagonal coefficient to zero
 241             matrixLower[faceIndex] = 0.0;
 242
 243             faceIndex++;
 244         }
 245
 246         // neighbour side
 247         faceIndex = startFaceNbr;
 248
 249         while (faceIndex < endFaceNbr)
 250         {
 251             // add contribution to source of owner (I am the neighbour)
 252             source[owner[losort[faceIndex]]] -=
 253                 matrixUpper[losort[faceIndex]]*value;
 254
 255             // set off-diagonal coefficient to zero
 256             matrixUpper[losort[faceIndex]] = 0.0;
 257
 258             faceIndex++;
 259         }
 260     }
 261 }
 262
 263
 264 template<class Type>
 265 template<template<class> class Matrix>
 266 void constraint<Type>::eliminateEquation(Matrix<Type>& matrix) const
 267 {
 268     eliminateEquation(matrix, rowID_, value_);
 269 }
 270
 271
 272 template<class Type>
 273 template<template<class> class Matrix>
 274 void constraint<Type>::eliminateEquation
 275 (
 276     Matrix<Type>& matrix,
 277     const direction d,
 278     scalarField& sourceCmpt
 279 ) const
 280 {
 281     const Type& fc = fixedComponents();
 282
 283     const scalar fcOfD = componentOfValue(fc, d);
 284
 285     if (fcOfD > SMALL)
 286     {
 287         const label startFaceOwn =
 288             matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
 289         const label endFaceOwn =
 290             matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_ + 1];
 291
 292         const label startFaceNbr =
 293             matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
 294         const label endFaceNbr =
 295             matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_ + 1];
 296
 297         const unallocLabelList& owner =
 298             matrix.psi().mesh().lduAddr().lowerAddr();
 299         const unallocLabelList& neighbour =
 300             matrix.psi().mesh().lduAddr().upperAddr();
 301         const unallocLabelList& losort =
 302             matrix.psi().mesh().lduAddr().losortAddr();
 303
 304         // My index =  rowID_
 305         if (matrix.symmetric())
 306         {
 307             // get the coefficients
 308             scalarField* coeffs = NULL;
 309
 310             if (matrix.hasUpper())
 311             {
 312                 coeffs = &matrix.upper();
 313             }
 314             else if (matrix.hasLower())
 315             {
 316                 coeffs = &matrix.lower();
 317             }
 318
 319             scalarField& matrixCoeffs = *coeffs;
 320
 321             // owner side
 322             label faceIndex = startFaceOwn;
 323
 324             while (faceIndex < endFaceOwn)
 325             {
 326                 // add contribution to source of the neighbour (I am the owner)
 327                 sourceCmpt[neighbour[faceIndex]] -=
 328                     fcOfD*matrixCoeffs[faceIndex]*componentOfValue(value(), d);
 329
 330                 // I am the owner =  owner[faceIndex]
 331                 // set off-diagonal coefficient to zero
 332                 matrixCoeffs[faceIndex] *= 1.0 - fcOfD;
 333
 334                 faceIndex++;
 335             }
 336
 337             // neighbour side
 338             faceIndex = startFaceNbr;
 339
 340             while (faceIndex < endFaceNbr)
 341             {
 342                 // add contribution to source of owner (I am the neighbour)
 343                 sourceCmpt[owner[losort[faceIndex]]] -=
 344                     fcOfD*matrixCoeffs[losort[faceIndex]]
 345                     *componentOfValue(value(), d);
 346
 347                 // I am the neighbour = neighbour[losort[faceIndex]]
 348                 // set off-diagonal coefficient to zero
 349                 matrixCoeffs[losort[faceIndex]] *= 1.0 - fcOfD;
 350
 351                 faceIndex++;
 352             }
 353         }
 354         else if (matrix.asymmetric())
 355         {
 356             scalarField& matrixUpper = matrix.upper();
 357
 358             scalarField& matrixLower = matrix.lower();
 359
 360             // owner side
 361             label faceIndex = startFaceOwn;
 362
 363             while (faceIndex < endFaceOwn)
 364             {
 365                 // add contribution to source of the neighbour (I am the owner)
 366                 sourceCmpt[neighbour[faceIndex]] -=
 367                     fcOfD*matrixLower[faceIndex]*componentOfValue(value(), d);
 368
 369                 // set off-diagonal coefficient to zero
 370                 matrixLower[faceIndex] *= 1.0 - fcOfD;
 371
 372                 faceIndex++;
 373             }
 374
 375             // neighbour side
 376             faceIndex = startFaceNbr;
 377
 378             while (faceIndex < endFaceNbr)
 379             {
 380                 // add contribution to source of owner (I am the neighbour)
 381                 sourceCmpt[owner[losort[faceIndex]]] -=
 382                     fcOfD*matrixUpper[losort[faceIndex]]
 383                     *componentOfValue(value(), d);
 384
 385                 // set off-diagonal coefficient to zero
 386                 matrixUpper[losort[faceIndex]] *= 1.0 - fcOfD;
 387
 388                 faceIndex++;
 389             }
 390         }
 391     }
 392 }
 393
 394
 395 template<class Type>
 396 template<template<class> class Matrix>
 397 void constraint<Type>::setSource
 398 (
 399     Matrix<Type>& matrix,
 400     const label rowID,
 401     const Type& value
 402 )
 403 {
 404     matrix.source()[rowID] = matrix.diag()[rowID]*value;
 405     matrix.psi()[rowID] = value;
 406 }
 407
 408
 409 template<class Type>
 410 template<template<class> class Matrix>
 411 void constraint<Type>::setSource(Matrix<Type>& matrix) const
 412 {
 413     setSource(matrix, rowID_, value_);
 414 }
 415
 416
 417 template<class Type>
 418 template<template<class> class Matrix>
 419 void constraint<Type>::setSourceDiag
 420 (
 421     Matrix<Type>& matrix,
 422     const direction d,
 423     scalarField& psiCmpt,
 424     scalarField& sourceCmpt
 425 ) const
 426 {
 427     const Type& fc = fixedComponents();
 428
 429     if (componentOfValue(fc, d) > SMALL)
 430     {
 431         sourceCmpt[rowID()] = componentOfValue(fc, d)*matrix.diag()[rowID()]
 432             *componentOfValue(value(), d);
 433
 434         // set the solution to the right value as well
 435         psiCmpt[rowID()] =
 436             componentOfValue(fc, d)*componentOfValue(value(), d);
 437     }
 438 }
 439
 440
 441 template<class Type>
 442 template<template<class> class Matrix>
 443 void constraint<Type>::reconstructMatrix
 444 (
 445     Matrix<Type>& matrix
 446 ) const
 447 {
 448     if (!matrixCoeffsSet_)
 449     {
 450         FatalErrorIn
 451         (
 452             "void constraint<Type>::reconstructMatrix("
 453             "Matrix<Type>& matrix)"
 454         )   << "matrix coefficients not set"
 455             << abort(FatalError);
 456     }
 457
 458     if (matrix.hasDiag())
 459     {
 460         matrix.diag()[rowID_] = diagCoeff_;
 461     }
 462
 463     const label startFaceOwn =
 464         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
 465
 466     const label startFaceNbr =
 467         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
 468
 469     const unallocLabelList& losort = matrix.psi().mesh().lduAddr().losortAddr();
 470
 471     if (matrix.hasUpper())
 472     {
 473         // get the upper coefficients
 474
 475         scalarField& matrixUpper = matrix.upper();
 476
 477         // owner side
 478         const scalarField& uOwn = upperCoeffsOwner();
 479
 480         label faceIndex = startFaceOwn;
 481
 482         forAll(uOwn, uOwnI)
 483         {
 484             matrixUpper[faceIndex] = uOwn[uOwnI];
 485
 486             faceIndex++;
 487         }
 488
 489         // neighbour side
 490         const scalarField& uNbr = upperCoeffsNeighbour();
 491
 492         faceIndex = startFaceNbr;
 493
 494         forAll (uNbr, uNbrI)
 495         {
 496             matrixUpper[losort[faceIndex]] = uNbr[uNbrI];
 497
 498             faceIndex++;
 499         }
 500     }
 501
 502     if (matrix.hasLower())
 503     {
 504         // get the lower coefficients
 505
 506         scalarField& matrixLower = matrix.lower();
 507
 508         // owner side
 509         const scalarField& lOwn = lowerCoeffsOwner();
 510
 511         label faceIndex = startFaceOwn;
 512
 513         forAll (lOwn, lOwnI)
 514         {
 515             matrixLower[faceIndex] = lOwn[lOwnI];
 516
 517             faceIndex++;
 518         }
 519
 520         // neighbour side
 521         const scalarField& lNbr = lowerCoeffsNeighbour();
 522
 523         faceIndex = startFaceNbr;
 524
 525         forAll (lNbr, lNbrI)
 526         {
 527             matrixLower[losort[faceIndex]] = lNbr[lNbrI];
 528
 529             faceIndex++;
 530         }
 531     }
 532 }
 533
 534
 535 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 536
 537 } // End namespace Foam
 538
 539 // ************************************************************************* //