src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgPolicy/BlockAamgPolicy.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 Class
  25     BlockAamgPolicy
  26
  27 Description
  28     Agglomerative AMG policy, adjusted for BlockLduMatrix
  29
  30 Author
  31     Klas Jareteg, 2012-12-13
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "BlockAamgPolicy.H"
  36 #include "coeffFields.H"
  37 #include "addToRunTimeSelectionTable.H"
  38 #include "BlockGAMGInterfaceField.H"
  39 #include "processorLduInterfaceField.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 template<class Type>
  44 const Foam::scalar Foam::BlockAamgPolicy<Type>::weightFactor_ = 0.65;
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 template<class Type>
  50 void Foam::BlockAamgPolicy<Type>::calcChild()
  51 {
  52     // Algorithm:
  53     // 1) Create temporary equation addressing using a double-pass algorithm.
  54     //    to create the offset table.
  55     // 2) Loop through all equations and for each equation find the best fit
  56     //    neighbour.  If all neighbours are grouped, add equation to best group
  57
  58     // Create row-based addressing
  59
  60     const label nRows = matrix_.lduAddr().size();
  61
  62     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
  63     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
  64
  65     // For each equation get number of coefficients in a row
  66     labelList cols(upperAddr.size() + lowerAddr.size());
  67     labelList cIndex(upperAddr.size() + lowerAddr.size());
  68     labelList rowOffset(nRows + 1, 0);
  69
  70     // Count the number of coefficients
  71     forAll (upperAddr, coeffI)
  72     {
  73         rowOffset[upperAddr[coeffI]]++;
  74     }
  75
  76     forAll (lowerAddr, coeffI)
  77     {
  78         rowOffset[lowerAddr[coeffI]]++;
  79     }
  80
  81     label nCoeffs = 0;
  82
  83     forAll (rowOffset, eqnI)
  84     {
  85         nCoeffs += rowOffset[eqnI];
  86     }
  87
  88     rowOffset[nRows] = nCoeffs;
  89
  90     for (label eqnI = nRows - 1; eqnI >= 0; --eqnI)
  91     {
  92         rowOffset[eqnI] = rowOffset[eqnI  + 1] - rowOffset[eqnI];
  93     }
  94
  95     rowOffset[0] = 0;
  96
  97     // Create column and coefficient index array
  98     {
  99         // Use child array to count number of entries per row.
 100         // Reset the list for counting
 101         labelList& nPerRow = child_;
 102         nPerRow = 0;
 103
 104         forAll (upperAddr, coeffI)
 105         {
 106             cols[rowOffset[upperAddr[coeffI]] + nPerRow[upperAddr[coeffI]]] =
 107                 lowerAddr[coeffI];
 108
 109             cIndex[rowOffset[upperAddr[coeffI]] + nPerRow[upperAddr[coeffI]]] =
 110                 coeffI;
 111
 112             nPerRow[upperAddr[coeffI]]++;
 113         }
 114
 115         forAll (lowerAddr, coeffI)
 116         {
 117             cols[rowOffset[lowerAddr[coeffI]] + nPerRow[lowerAddr[coeffI]]] =
 118                 upperAddr[coeffI];
 119
 120             cIndex[rowOffset[lowerAddr[coeffI]] + nPerRow[lowerAddr[coeffI]]] =
 121                 coeffI;
 122
 123             nPerRow[lowerAddr[coeffI]]++;
 124         }
 125
 126         // Reset child array
 127         child_ = -1;
 128     }
 129
 130
 131     // Calculate agglomeration
 132
 133     // Get matrix coefficients
 134     const CoeffField<Type>& diag = matrix_.diag();
 135     const CoeffField<Type>& upper = matrix_.upper();
 136
 137     labelList sizeOfGroups(nRows, 0);
 138
 139     nCoarseEqns_ = 0;
 140
 141     label indexU, indexG, colI, curEqn, nextEqn, groupPassI;
 142
 143     scalar magRowDiag, magColDiag, weight, weightU, weightG;
 144
 145     // Magnitudes pre-calculated
 146     Field<scalar> magDiag(diag.size());
 147     Field<scalar> magUpper(upper.size());
 148
 149     normPtr_->coeffMag(diag,magDiag);
 150     normPtr_->coeffMag(upper,magUpper);
 151
 152     for (label eqnI = 0; eqnI < nRows; eqnI++)
 153     {
 154         if (child_[eqnI] == -1)
 155         {
 156             curEqn = eqnI;
 157
 158             indexU = -1;
 159             indexG = -1;
 160
 161             child_[curEqn] = nCoarseEqns_;
 162
 163             magRowDiag = magDiag[curEqn];
 164
 165             for (groupPassI = 1; groupPassI < groupSize_; groupPassI++)
 166             {
 167                 weightU = 0;
 168                 weightG = 0;
 169
 170                 indexU = -1;
 171                 indexG = -1;
 172
 173                 for
 174                 (
 175                     label rowCoeffI = rowOffset[curEqn];
 176                     rowCoeffI < rowOffset[curEqn + 1];
 177                     rowCoeffI++
 178                 )
 179                 {
 180                     colI = cols[rowCoeffI];
 181
 182                     magColDiag = magDiag[colI];
 183
 184                     weight = (magUpper[cIndex[rowCoeffI]]
 185                         / max(magRowDiag, magColDiag));
 186
 187                     if (child_[colI] == -1)
 188                     {
 189                         if (indexU == -1 || weight > weightU)
 190                         {
 191                             indexU = rowCoeffI;
 192                             weightU = weight;
 193                         }
 194                     }
 195                     else if (child_[curEqn] != child_[colI])
 196                     {
 197                         if (indexG == -1 || weight > weightG)
 198                         {
 199                             indexG = rowCoeffI;
 200                             weightG = weight;
 201                         }
 202                     }
 203                 }
 204
 205                 if
 206                 (
 207                     indexU != -1
 208                  && (indexG == -1 || weightU >= weightFactor_*weightG)
 209                 )
 210                 {
 211                     // Found new element of group.  Add it and use as
 212                     // start of next search
 213
 214                     nextEqn = cols[indexU];
 215
 216                     child_[nextEqn] = child_[curEqn];
 217                     sizeOfGroups[child_[curEqn]];
 218
 219                     curEqn = nextEqn;
 220                 }
 221                 else
 222                 {
 223                     // Group full or cannot be extended
 224                     break;
 225                 }
 226             }
 227
 228             if
 229             (
 230                 groupPassI > 1
 231              || indexG == -1
 232              || sizeOfGroups[child_[cols[indexG]]] > (groupSize_ + 2)
 233             )
 234             {
 235                 sizeOfGroups[child_[eqnI]]++;
 236                 nCoarseEqns_++;
 237             }
 238             else
 239             {
 240                 child_[eqnI] = child_[cols[indexG]];
 241                 sizeOfGroups[child_[cols[indexG]]]++;
 242             }
 243         }
 244     }
 245
 246     // The decision on parallel agglomeration needs to be made for the
 247     // whole gang of processes; otherwise I may end up with a different
 248     // number of agglomeration levels on different processors.
 249
 250     if
 251     (
 252         nCoarseEqns_ > BlockAmgPolicy<Type>::minCoarseEqns()
 253      && 3*nCoarseEqns_ <= 2*nRows
 254     )
 255     {
 256         coarsen_ = true;
 257     }
 258
 259     reduce(coarsen_, andOp<bool>());
 260
 261 }
 262
 263
 264 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 265
 266 template<class Type>
 267 Foam::BlockAamgPolicy<Type>::BlockAamgPolicy
 268 (
 269     const BlockLduMatrix<Type>& matrix,
 270     const dictionary& dict,
 271     const label groupSize,
 272     const label minCoarseEqns
 273 )
 274 :
 275     BlockAmgPolicy<Type>(matrix, dict, groupSize, minCoarseEqns),
 276     matrix_(matrix),
 277     normPtr_
 278     (
 279         BlockCoeffNorm<Type>::New
 280         (
 281             dict
 282         )
 283     ),
 284     child_(matrix_.lduAddr().size()),
 285     groupSize_(groupSize),
 286     nCoarseEqns_(0),
 287     coarsen_(false)
 288 {
 289     calcChild();
 290 }
 291
 292 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 293
 294 template<class Type>
 295 Foam::BlockAamgPolicy<Type>::~BlockAamgPolicy()
 296 {}
 297
 298
 299 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 300
 301 template<class Type>
 302 Foam::autoPtr<Foam::BlockLduMatrix<Type> >
 303 Foam::BlockAamgPolicy<Type>::restrictMatrix() const
 304 {
 305     if (!coarsen_)
 306     {
 307         FatalErrorIn
 308         (
 309             "autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const"
 310         )   << "Requesting coarse matrix when it cannot be created"
 311             << abort(FatalError);
 312     }
 313
 314     // Construct the coarse matrix and ldu addressing for the next level
 315     // Algorithm:
 316     // 1) Loop through all fine coeffs. If the child labels on two sides are
 317     //    different, this creates a coarse coeff. Define owner and neighbour
 318     //    for this coeff based on cluster IDs.
 319     // 2) Check if the coeff has been seen before. If yes, add the coefficient
 320     //    to the appropriate field (stored with the equation). If no, create
 321     //    a new coeff with neighbour ID and add the coefficient
 322     // 3) Once all the coeffs have been created, loop through all clusters and
 323     //    insert the coeffs in the upper order. At the same time, collect the
 324     //    owner and neighbour addressing.
 325     // 4) Agglomerate the diagonal by summing up the fine diagonal
 326
 327     // Get addressing
 328     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
 329     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 330
 331     const label nFineCoeffs = upperAddr.size();
 332
 333 #   ifdef FULLDEBUG
 334     if (child_.size() != matrix_.lduAddr().size())
 335     {
 336         FatalErrorIn
 337         (
 338             "autoPtr<BlockLduMatrix<Type> > BlockAamgPolicy<Type>::"
 339             "restrictMatrix() const"
 340         )   << "Child array does not correspond to fine level. " << endl
 341             << " Child size: " << child_.size()
 342             << " number of equations: " << matrix_.lduAddr().size()
 343             << abort(FatalError);
 344     }
 345 #   endif
 346
 347
 348     // Storage for block neighbours and coefficients
 349
 350     // Guess initial maximum number of neighbours in block
 351     label maxNnbrs = 10;
 352
 353     // Number of neighbours per block
 354     labelList blockNnbrs(nCoarseEqns_, 0);
 355
 356     // Setup initial packed storage for neighbours and coefficients
 357     labelList blockNbrsData(maxNnbrs*nCoarseEqns_);
 358
 359     // Create face-restriction addressing
 360     labelList coeffRestrictAddr(nFineCoeffs);
 361
 362     // Initial neighbour array (not in upper-triangle order)
 363     labelList initCoarseNeighb(nFineCoeffs);
 364
 365     // Counter for coarse coeffs
 366     label nCoarseCoeffs = 0;
 367
 368     // Loop through all fine coeffs
 369     forAll (upperAddr, fineCoeffi)
 370     {
 371         label rmUpperAddr = child_[upperAddr[fineCoeffi]];
 372         label rmLowerAddr = child_[lowerAddr[fineCoeffi]];
 373
 374         if (rmUpperAddr == rmLowerAddr)
 375         {
 376             // For each fine coeff inside of a coarse cluster keep the address
 377             // of the cluster corresponding to the coeff in the
 378             // coeffRestrictAddr as a negative index
 379             coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1);
 380         }
 381         else
 382         {
 383             // This coeff is a part of a coarse coeff
 384
 385             label cOwn = rmUpperAddr;
 386             label cNei = rmLowerAddr;
 387
 388             // Get coarse owner and neighbour
 389             if (rmUpperAddr > rmLowerAddr)
 390             {
 391                 cOwn = rmLowerAddr;
 392                 cNei = rmUpperAddr;
 393             }
 394
 395             // Check the neighbour to see if this coeff has already been found
 396             bool nbrFound = false;
 397             label& ccnCoeffs = blockNnbrs[cOwn];
 398
 399             for (int i = 0; i < ccnCoeffs; i++)
 400             {
 401                 if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei)
 402                 {
 403                     nbrFound = true;
 404                     coeffRestrictAddr[fineCoeffi] =
 405                         blockNbrsData[maxNnbrs*cOwn + i];
 406                     break;
 407                 }
 408             }
 409
 410             if (!nbrFound)
 411             {
 412                 if (ccnCoeffs >= maxNnbrs)
 413                 {
 414                     // Double the size of list and copy data
 415                     label oldMaxNnbrs = maxNnbrs;
 416                     maxNnbrs *= 2;
 417
 418                     // Resize and copy list
 419                     const labelList oldBlockNbrsData = blockNbrsData;
 420                     blockNbrsData.setSize(maxNnbrs*nCoarseEqns_);
 421
 422                     forAll (blockNnbrs, i)
 423                     {
 424                         for (int j = 0; j < blockNnbrs[i]; j++)
 425                         {
 426                             blockNbrsData[maxNnbrs*i + j] =
 427                                 oldBlockNbrsData[oldMaxNnbrs*i + j];
 428                         }
 429                     }
 430                 }
 431
 432                 blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs;
 433                 initCoarseNeighb[nCoarseCoeffs] = cNei;
 434                 coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs;
 435                 ccnCoeffs++;
 436
 437                 // New coarse coeff created
 438                 nCoarseCoeffs++;
 439             }
 440         }
 441     } // End for all fine coeffs
 442
 443
 444     // Renumber into upper-triangular order
 445
 446     // All coarse owner-neighbour storage
 447     labelList coarseOwner(nCoarseCoeffs);
 448     labelList coarseNeighbour(nCoarseCoeffs);
 449     labelList coarseCoeffMap(nCoarseCoeffs);
 450
 451     label coarseCoeffi = 0;
 452
 453     forAll (blockNnbrs, cci)
 454     {
 455         label* cCoeffs = &blockNbrsData[maxNnbrs*cci];
 456         label ccnCoeffs = blockNnbrs[cci];
 457
 458         for (int i = 0; i < ccnCoeffs; i++)
 459         {
 460             coarseOwner[coarseCoeffi] = cci;
 461             coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]];
 462             coarseCoeffMap[cCoeffs[i]] = coarseCoeffi;
 463             coarseCoeffi++;
 464         }
 465     }
 466
 467     forAll(coeffRestrictAddr, fineCoeffi)
 468     {
 469         if (coeffRestrictAddr[fineCoeffi] >= 0)
 470         {
 471             coeffRestrictAddr[fineCoeffi] =
 472                 coarseCoeffMap[coeffRestrictAddr[fineCoeffi]];
 473         }
 474     }
 475
 476     // Clear the temporary storage for the coarse matrix data
 477     blockNnbrs.setSize(0);
 478     blockNbrsData.setSize(0);
 479     initCoarseNeighb.setSize(0);
 480     coarseCoeffMap.setSize(0);
 481
 482     // Create coarse-level coupled interfaces
 483
 484     // Create coarse interfaces, addressing and coefficients
 485     const label interfaceSize =
 486         const_cast<BlockLduMatrix<Type>& >(matrix_).interfaces().size();
 487
 488     const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
 489         interfaceFields =
 490         const_cast<BlockLduMatrix<Type>&>(matrix_).interfaces();
 491
 492     // Set the coarse interfaces and coefficients
 493     lduInterfacePtrsList* coarseInterfacesPtr =
 494         new lduInterfacePtrsList(interfaceSize);
 495     lduInterfacePtrsList& coarseInterfaces = *coarseInterfacesPtr;
 496
 497     labelListList coarseInterfaceAddr(interfaceSize);
 498
 499     // Add the coarse level
 500
 501     // Set the coarse ldu addressing onto the list
 502     lduPrimitiveMesh* coarseAddrPtr =
 503         new lduPrimitiveMesh
 504         (
 505             nCoarseEqns_,
 506             coarseOwner,
 507             coarseNeighbour,
 508             true
 509         );
 510
 511     // Initialise transfer of restrict addressing on the interface
 512     forAll (interfaceFields, intI)
 513     {
 514         if (interfaceFields.set(intI))
 515         {
 516             interfaceFields[intI].interface().initInternalFieldTransfer
 517             (
 518                 Pstream::blocking,
 519                 child_
 520             );
 521         }
 522     }
 523
 524     // Store coefficients to avoid tangled communications
 525     // HJ, 1/Apr/2009
 526     FieldField<Field, label> fineInterfaceAddr(interfaceFields.size());
 527
 528     forAll (interfaceFields, intI)
 529     {
 530         if (interfaceFields.set(intI))
 531         {
 532             const lduInterface& fineInterface =
 533                 interfaceFields[intI].interface();
 534
 535             fineInterfaceAddr.set
 536             (
 537                 intI,
 538                 new labelField
 539                 (
 540                     fineInterface.internalFieldTransfer
 541                     (
 542                         Pstream::blocking,
 543                         child_
 544                     )
 545                 )
 546             );
 547         }
 548     }
 549
 550     // Create GAMG interfaces
 551     forAll (interfaceFields, intI)
 552     {
 553         if (interfaceFields.set(intI))
 554         {
 555             const lduInterface& fineInterface =
 556                 interfaceFields[intI].interface();
 557
 558             coarseInterfaces.set
 559             (
 560                 intI,
 561                 GAMGInterface::New
 562                 (
 563                     *coarseAddrPtr,
 564                     fineInterface,
 565                     fineInterface.interfaceInternalField(child_),
 566                     fineInterfaceAddr[intI]
 567                 ).ptr()
 568             );
 569         }
 570     }
 571
 572     forAll (interfaceFields, intI)
 573     {
 574         if (interfaceFields.set(intI))
 575         {
 576             const GAMGInterface& coarseInterface =
 577                 refCast<const GAMGInterface>(coarseInterfaces[intI]);
 578
 579             coarseInterfaceAddr[intI] = coarseInterface.faceCells();
 580         }
 581     }
 582
 583     // Add interfaces
 584     coarseAddrPtr->addInterfaces
 585     (
 586         *coarseInterfacesPtr,
 587         coarseInterfaceAddr,
 588         matrix_.patchSchedule()
 589     );
 590
 591     // Set the coarse level matrix
 592     BlockLduMatrix<Type>* coarseMatrixPtr =
 593         new BlockLduMatrix<Type>(*coarseAddrPtr);
 594     BlockLduMatrix<Type>& coarseMatrix = *coarseMatrixPtr;
 595
 596     typename BlockLduInterfaceFieldPtrsList<Type>::Type&
 597         coarseInterfaceFieldsTransfer =
 598           coarseMatrix.interfaces();
 599
 600     // Aggolmerate the upper and lower coupled coefficients
 601     forAll (interfaceFields, intI)
 602     {
 603         if (interfaceFields.set(intI))
 604         {
 605             const GAMGInterface& coarseInterface =
 606                 refCast<const GAMGInterface>(coarseInterfaces[intI]);
 607
 608             coarseInterfaceFieldsTransfer.set
 609             (
 610                 intI,
 611                 BlockGAMGInterfaceField<Type>::New
 612                 (
 613                     coarseInterface,
 614                     (interfaceFields[intI])
 615                 ).ptr()
 616             );
 617
 618             coarseMatrix.coupleUpper().set
 619             (
 620                 intI,
 621                 coarseInterface.agglomerateBlockCoeffs
 622                 (
 623                     matrix_.coupleUpper()[intI]
 624                 )
 625             );
 626
 627             coarseMatrix.coupleLower().set
 628             (
 629                 intI,
 630                 coarseInterface.agglomerateBlockCoeffs
 631                 (
 632                     matrix_.coupleLower()[intI]
 633                 )
 634             );
 635         }
 636     }
 637
 638     // Matrix restriction done!
 639
 640     typedef CoeffField<Type> TypeCoeffField;
 641
 642     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 643
 644     TypeCoeffField& coarseUpper = coarseMatrix.upper();
 645     TypeCoeffField& coarseDiag = coarseMatrix.diag();
 646     const TypeCoeffField& fineUpper = matrix_.upper();
 647     const TypeCoeffField& fineDiag = matrix_.diag();
 648
 649     // KRJ: 2013-01-31: Many cases needed as there are different combinations
 650     if (matrix_.asymmetric())
 651     {
 652         TypeCoeffField& coarseLower = coarseMatrix.lower();
 653         const TypeCoeffField& fineLower = matrix_.lower();
 654
 655         if (fineDiag.activeType() == blockCoeffBase::SQUARE)
 656         {
 657             if (fineUpper.activeType() == blockCoeffBase::SQUARE)
 658             {
 659                 squareTypeField& activeCoarseUpper = coarseUpper.asSquare();
 660                 squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
 661                 const squareTypeField& activeFineUpper = fineUpper.asSquare();
 662
 663                 squareTypeField& activeCoarseLower = coarseLower.asSquare();
 664                 const squareTypeField& activeFineLower = fineLower.asSquare();
 665
 666                 restrictResidual(fineDiag, coarseDiag);
 667
 668                 forAll(coeffRestrictAddr, fineCoeffI)
 669                 {
 670                     label cCoeff = coeffRestrictAddr[fineCoeffI];
 671
 672                     if (cCoeff >= 0)
 673                     {
 674                         activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI];
 675                         activeCoarseLower[cCoeff] += activeFineLower[fineCoeffI];
 676                     }
 677                     else
 678                     {
 679                         // Add the fine face coefficients into the diagonal.
 680                         activeCoarseDiag[-1 - cCoeff] +=
 681                             activeFineUpper[fineCoeffI] + activeFineLower[fineCoeffI];
 682                     }
 683                 }
 684             }
 685             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
 686             {
 687                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 688                     << "Matrix diagonal of square type and upper of linear type is not implemented"
 689                     << abort(FatalError);
 690             }
 691             else
 692             {
 693                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 694                     << "Matrix diagonal of square type and upper of scalar type is not implemented"
 695                     << abort(FatalError);
 696             }
 697         }
 698         else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
 699         {
 700                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 701                     << "Matrix diagonal of linear type not implemented"
 702                     << abort(FatalError);
 703         }
 704         else
 705         {
 706                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 707                     << "Matrix diagonal of scalar type not implemented"
 708                     << abort(FatalError);
 709         }
 710     }
 711     else
 712     {
 713         if (fineDiag.activeType() == blockCoeffBase::SQUARE)
 714         {
 715             if (fineUpper.activeType() == blockCoeffBase::SQUARE)
 716             {
 717                 squareTypeField& activeCoarseUpper = coarseUpper.asSquare();
 718                 squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
 719                 const squareTypeField& activeFineUpper = fineUpper.asSquare();
 720
 721                 restrictResidual(fineDiag, coarseDiag);
 722
 723                 forAll(coeffRestrictAddr, fineCoeffI)
 724                 {
 725                     label cCoeff = coeffRestrictAddr[fineCoeffI];
 726
 727                     if (cCoeff >= 0)
 728                     {
 729                         activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI];
 730                     }
 731                     else
 732                     {
 733                         // Add the fine face coefficient into the diagonal.
 734                         activeCoarseDiag[-1 - cCoeff] += 2*activeFineUpper[fineCoeffI];
 735                     }
 736                 }
 737             }
 738             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
 739             {
 740                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 741                     << "Matrix diagonal of square type and upper of linear type is not implemented"
 742                     << abort(FatalError);
 743             }
 744             else
 745             {
 746                 FatalErrorIn("autoPtr<amgMatrix> BlockAamgPolicy<Type>::restrictMatrix() const")
 747                     << "Matrix diagonal of square type and upper of scalar type is not implemented"
 748                     << abort(FatalError);
 749             }
 750
 751         }
 752         else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
 753         {
 754                 FatalErrorIn
 755                 (
 756                     "autoPtr<amgMatrix> BlockAamgPolicy<Type>::"
 757                     "restrictMatrix() const"
 758                 )   << "Matrix diagonal of linear type not implemented"
 759                     << abort(FatalError);
 760         }
 761         else
 762         {
 763                 FatalErrorIn
 764                 (
 765                     "autoPtr<amgMatrix> BlockAamgPolicy<Type>::"
 766                     "restrictMatrix() const"
 767                 )   << "Matrix diagonal of scalar type not implemented"
 768                     << abort(FatalError);
 769         }
 770     }
 771
 772     return autoPtr<BlockLduMatrix<Type> >
 773     (
 774         new BlockLduMatrix<Type>
 775         (
 776             coarseMatrix
 777         )
 778     );
 779 }
 780
 781
 782 template<class Type>
 783 void Foam::BlockAamgPolicy<Type>::restrictResidual
 784 (
 785     const CoeffField<Type>& res,
 786     CoeffField<Type>& coarseRes
 787 ) const
 788 {
 789     typedef CoeffField<Type> TypeCoeffField;
 790
 791     if (res.activeType() == blockCoeffBase::SQUARE &&
 792             coarseRes.activeType() == blockCoeffBase::SQUARE)
 793     {
 794         typedef typename TypeCoeffField::squareTypeField squareTypeField;
 795
 796         squareTypeField& activeCoarseRes = coarseRes.asSquare();
 797         const squareTypeField& activeRes = res.asSquare();
 798
 799         // KRJ: 2013-02-05: To be done by correct pTraits?
 800         forAll (coarseRes, i)
 801         {
 802             activeCoarseRes[i] *= 0.0;
 803         }
 804
 805         forAll (res, i)
 806         {
 807             activeCoarseRes[child_[i]] += activeRes[i];
 808         }
 809     }
 810     else
 811     {
 812         FatalErrorIn("void  BlockAamgPolicy<Type>::restrictResidual() const")
 813             << "Only present for square type coeff type"
 814             << abort(FatalError);
 815     }
 816 }
 817
 818
 819 template<class Type>
 820 void Foam::BlockAamgPolicy<Type>::restrictResidual
 821 (
 822     const Field<Type>& res,
 823     Field<Type>& coarseRes
 824 ) const
 825 {
 826     coarseRes = pTraits<Type>::zero;
 827
 828     forAll (res, i)
 829     {
 830         coarseRes[child_[i]] += res[i];
 831     }
 832 }
 833
 834
 835 template<class Type>
 836 void Foam::BlockAamgPolicy<Type>::prolongateCorrection
 837 (
 838     Field<Type>& x,
 839     const Field<Type>& coarseX
 840 ) const
 841 {
 842     forAll (x, i)
 843     {
 844         x[i] += coarseX[child_[i]];
 845     }
 846 }
 847
 848
 849 // ************************************************************************* //