src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.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     BlockMatrixAgglomeration
  26
  27 Description
  28     Agglomerative block matrix AMG corsening, adjusted for BlockLduMatrix
  29
  30 Author
  31     Klas Jareteg, 2012-12-13
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "BlockMatrixAgglomeration.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::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65;
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 template<class Type>
  50 void Foam::BlockMatrixAgglomeration<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_ > BlockMatrixCoarsening<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::BlockMatrixAgglomeration<Type>::BlockMatrixAgglomeration
 268 (
 269     const BlockLduMatrix<Type>& matrix,
 270     const dictionary& dict,
 271     const label groupSize,
 272     const label minCoarseEqns
 273 )
 274 :
 275     BlockMatrixCoarsening<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::BlockMatrixAgglomeration<Type>::~BlockMatrixAgglomeration()
 296 {}
 297
 298
 299 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 300
 301 template<class Type>
 302 Foam::autoPtr<Foam::BlockLduMatrix<Type> >
 303 Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
 304 {
 305     if (!coarsen_)
 306     {
 307         FatalErrorIn
 308         (
 309             "autoPtr<amgMatrix> "
 310             "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
 311         )   << "Requesting coarse matrix when it cannot be created"
 312             << abort(FatalError);
 313     }
 314
 315     // Construct the coarse matrix and ldu addressing for the next level
 316     // Algorithm:
 317     // 1) Loop through all fine coeffs. If the child labels on two sides are
 318     //    different, this creates a coarse coeff. Define owner and neighbour
 319     //    for this coeff based on cluster IDs.
 320     // 2) Check if the coeff has been seen before. If yes, add the coefficient
 321     //    to the appropriate field (stored with the equation). If no, create
 322     //    a new coeff with neighbour ID and add the coefficient
 323     // 3) Once all the coeffs have been created, loop through all clusters and
 324     //    insert the coeffs in the upper order. At the same time, collect the
 325     //    owner and neighbour addressing.
 326     // 4) Agglomerate the diagonal by summing up the fine diagonal
 327
 328     // Get addressing
 329     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
 330     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 331
 332     const label nFineCoeffs = upperAddr.size();
 333
 334 #   ifdef FULLDEBUG
 335     if (child_.size() != matrix_.lduAddr().size())
 336     {
 337         FatalErrorIn
 338         (
 339             "autoPtr<BlockLduMatrix<Type> > BlockMatrixAgglomeration<Type>::"
 340             "restrictMatrix() const"
 341         )   << "Child array does not correspond to fine level. " << endl
 342             << " Child size: " << child_.size()
 343             << " number of equations: " << matrix_.lduAddr().size()
 344             << abort(FatalError);
 345     }
 346 #   endif
 347
 348
 349     // Storage for block neighbours and coefficients
 350
 351     // Guess initial maximum number of neighbours in block
 352     label maxNnbrs = 10;
 353
 354     // Number of neighbours per block
 355     labelList blockNnbrs(nCoarseEqns_, 0);
 356
 357     // Setup initial packed storage for neighbours and coefficients
 358     labelList blockNbrsData(maxNnbrs*nCoarseEqns_);
 359
 360     // Create face-restriction addressing
 361     labelList coeffRestrictAddr(nFineCoeffs);
 362
 363     // Initial neighbour array (not in upper-triangle order)
 364     labelList initCoarseNeighb(nFineCoeffs);
 365
 366     // Counter for coarse coeffs
 367     label nCoarseCoeffs = 0;
 368
 369     // Loop through all fine coeffs
 370     forAll (upperAddr, fineCoeffi)
 371     {
 372         label rmUpperAddr = child_[upperAddr[fineCoeffi]];
 373         label rmLowerAddr = child_[lowerAddr[fineCoeffi]];
 374
 375         if (rmUpperAddr == rmLowerAddr)
 376         {
 377             // For each fine coeff inside of a coarse cluster keep the address
 378             // of the cluster corresponding to the coeff in the
 379             // coeffRestrictAddr as a negative index
 380             coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1);
 381         }
 382         else
 383         {
 384             // This coeff is a part of a coarse coeff
 385
 386             label cOwn = rmUpperAddr;
 387             label cNei = rmLowerAddr;
 388
 389             // Get coarse owner and neighbour
 390             if (rmUpperAddr > rmLowerAddr)
 391             {
 392                 cOwn = rmLowerAddr;
 393                 cNei = rmUpperAddr;
 394             }
 395
 396             // Check the neighbour to see if this coeff has already been found
 397             bool nbrFound = false;
 398             label& ccnCoeffs = blockNnbrs[cOwn];
 399
 400             for (int i = 0; i < ccnCoeffs; i++)
 401             {
 402                 if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei)
 403                 {
 404                     nbrFound = true;
 405                     coeffRestrictAddr[fineCoeffi] =
 406                         blockNbrsData[maxNnbrs*cOwn + i];
 407                     break;
 408                 }
 409             }
 410
 411             if (!nbrFound)
 412             {
 413                 if (ccnCoeffs >= maxNnbrs)
 414                 {
 415                     // Double the size of list and copy data
 416                     label oldMaxNnbrs = maxNnbrs;
 417                     maxNnbrs *= 2;
 418
 419                     // Resize and copy list
 420                     const labelList oldBlockNbrsData = blockNbrsData;
 421                     blockNbrsData.setSize(maxNnbrs*nCoarseEqns_);
 422
 423                     forAll (blockNnbrs, i)
 424                     {
 425                         for (int j = 0; j < blockNnbrs[i]; j++)
 426                         {
 427                             blockNbrsData[maxNnbrs*i + j] =
 428                                 oldBlockNbrsData[oldMaxNnbrs*i + j];
 429                         }
 430                     }
 431                 }
 432
 433                 blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs;
 434                 initCoarseNeighb[nCoarseCoeffs] = cNei;
 435                 coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs;
 436                 ccnCoeffs++;
 437
 438                 // New coarse coeff created
 439                 nCoarseCoeffs++;
 440             }
 441         }
 442     } // End for all fine coeffs
 443
 444
 445     // Renumber into upper-triangular order
 446
 447     // All coarse owner-neighbour storage
 448     labelList coarseOwner(nCoarseCoeffs);
 449     labelList coarseNeighbour(nCoarseCoeffs);
 450     labelList coarseCoeffMap(nCoarseCoeffs);
 451
 452     label coarseCoeffi = 0;
 453
 454     forAll (blockNnbrs, cci)
 455     {
 456         label* cCoeffs = &blockNbrsData[maxNnbrs*cci];
 457         label ccnCoeffs = blockNnbrs[cci];
 458
 459         for (int i = 0; i < ccnCoeffs; i++)
 460         {
 461             coarseOwner[coarseCoeffi] = cci;
 462             coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]];
 463             coarseCoeffMap[cCoeffs[i]] = coarseCoeffi;
 464             coarseCoeffi++;
 465         }
 466     }
 467
 468     forAll(coeffRestrictAddr, fineCoeffi)
 469     {
 470         if (coeffRestrictAddr[fineCoeffi] >= 0)
 471         {
 472             coeffRestrictAddr[fineCoeffi] =
 473                 coarseCoeffMap[coeffRestrictAddr[fineCoeffi]];
 474         }
 475     }
 476
 477     // Clear the temporary storage for the coarse matrix data
 478     blockNnbrs.setSize(0);
 479     blockNbrsData.setSize(0);
 480     initCoarseNeighb.setSize(0);
 481     coarseCoeffMap.setSize(0);
 482
 483     // Create coarse-level coupled interfaces
 484
 485     // Create coarse interfaces, addressing and coefficients
 486     const label interfaceSize =
 487         const_cast<BlockLduMatrix<Type>& >(matrix_).interfaces().size();
 488
 489     const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
 490         interfaceFields =
 491         const_cast<BlockLduMatrix<Type>&>(matrix_).interfaces();
 492
 493     // Set the coarse interfaces and coefficients
 494     lduInterfacePtrsList* coarseInterfacesPtr =
 495         new lduInterfacePtrsList(interfaceSize);
 496     lduInterfacePtrsList& coarseInterfaces = *coarseInterfacesPtr;
 497
 498     labelListList coarseInterfaceAddr(interfaceSize);
 499
 500     // Add the coarse level
 501
 502     // Set the coarse ldu addressing onto the list
 503     lduPrimitiveMesh* coarseAddrPtr =
 504         new lduPrimitiveMesh
 505         (
 506             nCoarseEqns_,
 507             coarseOwner,
 508             coarseNeighbour,
 509             true
 510         );
 511
 512     // Initialise transfer of restrict addressing on the interface
 513     forAll (interfaceFields, intI)
 514     {
 515         if (interfaceFields.set(intI))
 516         {
 517             interfaceFields[intI].interface().initInternalFieldTransfer
 518             (
 519                 Pstream::blocking,
 520                 child_
 521             );
 522         }
 523     }
 524
 525     // Store coefficients to avoid tangled communications
 526     // HJ, 1/Apr/2009
 527     FieldField<Field, label> fineInterfaceAddr(interfaceFields.size());
 528
 529     forAll (interfaceFields, intI)
 530     {
 531         if (interfaceFields.set(intI))
 532         {
 533             const lduInterface& fineInterface =
 534                 interfaceFields[intI].interface();
 535
 536             fineInterfaceAddr.set
 537             (
 538                 intI,
 539                 new labelField
 540                 (
 541                     fineInterface.internalFieldTransfer
 542                     (
 543                         Pstream::blocking,
 544                         child_
 545                     )
 546                 )
 547             );
 548         }
 549     }
 550
 551     // Create GAMG interfaces
 552     forAll (interfaceFields, intI)
 553     {
 554         if (interfaceFields.set(intI))
 555         {
 556             const lduInterface& fineInterface =
 557                 interfaceFields[intI].interface();
 558
 559             coarseInterfaces.set
 560             (
 561                 intI,
 562                 GAMGInterface::New
 563                 (
 564                     *coarseAddrPtr,
 565                     fineInterface,
 566                     fineInterface.interfaceInternalField(child_),
 567                     fineInterfaceAddr[intI]
 568                 ).ptr()
 569             );
 570         }
 571     }
 572
 573     forAll (interfaceFields, intI)
 574     {
 575         if (interfaceFields.set(intI))
 576         {
 577             const GAMGInterface& coarseInterface =
 578                 refCast<const GAMGInterface>(coarseInterfaces[intI]);
 579
 580             coarseInterfaceAddr[intI] = coarseInterface.faceCells();
 581         }
 582     }
 583
 584     // Add interfaces
 585     coarseAddrPtr->addInterfaces
 586     (
 587         *coarseInterfacesPtr,
 588         coarseInterfaceAddr,
 589         matrix_.patchSchedule()
 590     );
 591
 592     // Set the coarse level matrix
 593     BlockLduMatrix<Type>* coarseMatrixPtr =
 594         new BlockLduMatrix<Type>(*coarseAddrPtr);
 595     BlockLduMatrix<Type>& coarseMatrix = *coarseMatrixPtr;
 596
 597     typename BlockLduInterfaceFieldPtrsList<Type>::Type&
 598         coarseInterfaceFieldsTransfer =
 599           coarseMatrix.interfaces();
 600
 601     // Aggolmerate the upper and lower coupled coefficients
 602     forAll (interfaceFields, intI)
 603     {
 604         if (interfaceFields.set(intI))
 605         {
 606             const GAMGInterface& coarseInterface =
 607                 refCast<const GAMGInterface>(coarseInterfaces[intI]);
 608
 609             coarseInterfaceFieldsTransfer.set
 610             (
 611                 intI,
 612                 BlockGAMGInterfaceField<Type>::New
 613                 (
 614                     coarseInterface,
 615                     (interfaceFields[intI])
 616                 ).ptr()
 617             );
 618
 619             coarseMatrix.coupleUpper().set
 620             (
 621                 intI,
 622                 coarseInterface.agglomerateBlockCoeffs
 623                 (
 624                     matrix_.coupleUpper()[intI]
 625                 )
 626             );
 627
 628             coarseMatrix.coupleLower().set
 629             (
 630                 intI,
 631                 coarseInterface.agglomerateBlockCoeffs
 632                 (
 633                     matrix_.coupleLower()[intI]
 634                 )
 635             );
 636         }
 637     }
 638
 639     // Matrix restriction done!
 640
 641     typedef CoeffField<Type> TypeCoeffField;
 642
 643     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 644
 645     TypeCoeffField& coarseUpper = coarseMatrix.upper();
 646     TypeCoeffField& coarseDiag = coarseMatrix.diag();
 647     const TypeCoeffField& fineUpper = matrix_.upper();
 648     const TypeCoeffField& fineDiag = matrix_.diag();
 649
 650     // KRJ: 2013-01-31: Many cases needed as there are different combinations
 651     if (matrix_.asymmetric())
 652     {
 653         TypeCoeffField& coarseLower = coarseMatrix.lower();
 654         const TypeCoeffField& fineLower = matrix_.lower();
 655
 656         if (fineDiag.activeType() == blockCoeffBase::SQUARE)
 657         {
 658             if (fineUpper.activeType() == blockCoeffBase::SQUARE)
 659             {
 660                 squareTypeField& activeCoarseUpper = coarseUpper.asSquare();
 661                 squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
 662                 const squareTypeField& activeFineUpper = fineUpper.asSquare();
 663
 664                 squareTypeField& activeCoarseLower = coarseLower.asSquare();
 665                 const squareTypeField& activeFineLower = fineLower.asSquare();
 666
 667                 restrictResidual(fineDiag, coarseDiag);
 668
 669                 forAll(coeffRestrictAddr, fineCoeffI)
 670                 {
 671                     label cCoeff = coeffRestrictAddr[fineCoeffI];
 672
 673                     if (cCoeff >= 0)
 674                     {
 675                         activeCoarseUpper[cCoeff]
 676                             += activeFineUpper[fineCoeffI];
 677                         activeCoarseLower[cCoeff]
 678                             += activeFineLower[fineCoeffI];
 679                     }
 680                     else
 681                     {
 682                         // Add the fine face coefficients into the diagonal.
 683                         activeCoarseDiag[-1 - cCoeff] +=
 684                             activeFineUpper[fineCoeffI]
 685                           + activeFineLower[fineCoeffI];
 686                     }
 687                 }
 688             }
 689             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
 690             {
 691                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 692                     << "Matrix diagonal of square type and upper of linear type is not implemented"
 693                     << abort(FatalError);
 694             }
 695             else
 696             {
 697                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 698                     << "Matrix diagonal of square type and upper of scalar type is not implemented"
 699                     << abort(FatalError);
 700             }
 701         }
 702         else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
 703         {
 704                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 705                     << "Matrix diagonal of linear type not implemented"
 706                     << abort(FatalError);
 707         }
 708         else
 709         {
 710                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 711                     << "Matrix diagonal of scalar type not implemented"
 712                     << abort(FatalError);
 713         }
 714     }
 715     else
 716     {
 717         if (fineDiag.activeType() == blockCoeffBase::SQUARE)
 718         {
 719             if (fineUpper.activeType() == blockCoeffBase::SQUARE)
 720             {
 721                 squareTypeField& activeCoarseUpper = coarseUpper.asSquare();
 722                 squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
 723                 const squareTypeField& activeFineUpper = fineUpper.asSquare();
 724
 725                 restrictResidual(fineDiag, coarseDiag);
 726
 727                 forAll(coeffRestrictAddr, fineCoeffI)
 728                 {
 729                     label cCoeff = coeffRestrictAddr[fineCoeffI];
 730
 731                     if (cCoeff >= 0)
 732                     {
 733                         activeCoarseUpper[cCoeff]
 734                             += activeFineUpper[fineCoeffI];
 735                     }
 736                     else
 737                     {
 738                         // Add the fine face coefficient into the diagonal.
 739                         activeCoarseDiag[-1 - cCoeff]
 740                             += 2*activeFineUpper[fineCoeffI];
 741                     }
 742                 }
 743             }
 744             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
 745             {
 746                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 747                     << "Matrix diagonal of square type and upper of linear type is not implemented"
 748                     << abort(FatalError);
 749             }
 750             else
 751             {
 752                 FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
 753                     << "Matrix diagonal of square type and upper of scalar type is not implemented"
 754                     << abort(FatalError);
 755             }
 756
 757         }
 758         else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
 759         {
 760                 FatalErrorIn
 761                 (
 762                     "autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::"
 763                     "restrictMatrix() const"
 764                 )   << "Matrix diagonal of linear type not implemented"
 765                     << abort(FatalError);
 766         }
 767         else
 768         {
 769                 FatalErrorIn
 770                 (
 771                     "autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::"
 772                     "restrictMatrix() const"
 773                 )   << "Matrix diagonal of scalar type not implemented"
 774                     << abort(FatalError);
 775         }
 776     }
 777
 778     return autoPtr<BlockLduMatrix<Type> >
 779     (
 780         new BlockLduMatrix<Type>
 781         (
 782             coarseMatrix
 783         )
 784     );
 785 }
 786
 787
 788 template<class Type>
 789 void Foam::BlockMatrixAgglomeration<Type>::restrictResidual
 790 (
 791     const CoeffField<Type>& res,
 792     CoeffField<Type>& coarseRes
 793 ) const
 794 {
 795     typedef CoeffField<Type> TypeCoeffField;
 796
 797     if (res.activeType() == blockCoeffBase::SQUARE &&
 798             coarseRes.activeType() == blockCoeffBase::SQUARE)
 799     {
 800         typedef typename TypeCoeffField::squareTypeField squareTypeField;
 801
 802         squareTypeField& activeCoarseRes = coarseRes.asSquare();
 803         const squareTypeField& activeRes = res.asSquare();
 804
 805         // KRJ: 2013-02-05: To be done by correct pTraits?
 806         forAll (coarseRes, i)
 807         {
 808             activeCoarseRes[i] *= 0.0;
 809         }
 810
 811         forAll (res, i)
 812         {
 813             activeCoarseRes[child_[i]] += activeRes[i];
 814         }
 815     }
 816     else
 817     {
 818         FatalErrorIn("void  BlockMatrixAgglomeration<Type>::restrictResidual() const")
 819             << "Only present for square type coeff type"
 820             << abort(FatalError);
 821     }
 822 }
 823
 824
 825 template<class Type>
 826 void Foam::BlockMatrixAgglomeration<Type>::restrictResidual
 827 (
 828     const Field<Type>& res,
 829     Field<Type>& coarseRes
 830 ) const
 831 {
 832     coarseRes = pTraits<Type>::zero;
 833
 834     forAll (res, i)
 835     {
 836         coarseRes[child_[i]] += res[i];
 837     }
 838 }
 839
 840
 841 template<class Type>
 842 void Foam::BlockMatrixAgglomeration<Type>::prolongateCorrection
 843 (
 844     Field<Type>& x,
 845     const Field<Type>& coarseX
 846 ) const
 847 {
 848     forAll (x, i)
 849     {
 850         x[i] += coarseX[child_[i]];
 851     }
 852 }
 853
 854
 855 // ************************************************************************* //