src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
   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 \*----------------------------------------------------------------------------*/
  26
  27 #include "referredCellList.H"
  28 #include "interactionLists.H"
  29 #include "polyBoundaryMeshEntries.H"
  30 #include "PstreamCombineReduceOps.H"
  31 #include "Time.H"
  32 #include "globalMeshData.H"
  33 #include "processorPolyPatch.H"
  34 #include "cyclicPolyPatch.H"
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 void Foam::referredCellList::buildReferredCellList
  39 (
  40     bool pointPointListBuild
  41 )
  42 {
  43     Info << nl << "Building list of referred interaction neighbours" << endl;
  44
  45     const polyMesh& mesh(il_.mesh());
  46
  47     DynamicList<referredCell> referredInteractionList;
  48
  49     // realCellsWithinRCutMaxOfAnyReferringPatch
  50     DynamicList<label> rCellsWRRP;
  51
  52     // realFacesWithinRCutMaxOfAnyReferringPatch
  53     DynamicList<label> rFacesWRRP;
  54
  55     // realEdgesWithinRCutMaxOfAnyReferringPatch
  56     DynamicList<label> rEdgesWRRP;
  57
  58     // realPointsWithinRCutMaxOfAnyReferringPatch
  59     DynamicList<label> rPointsWRRP;
  60
  61     labelListList processorPatchSegmentMapping
  62     (
  63         mesh.globalData().processorPatches().size()
  64     );
  65
  66     List<vectorList> allNeighbourFaceCentres
  67     (
  68         mesh.globalData().processorPatches().size()
  69     );
  70
  71     List<vectorList> allNeighbourFaceAreas
  72     (
  73         mesh.globalData().processorPatches().size()
  74     );
  75
  76     label nUndecomposedPatches = 0;
  77
  78     if (Pstream::parRun())
  79     {
  80         dictionary patchDictionary;
  81
  82         DynamicList<word> patchNames;
  83
  84         Time undecomposedTime
  85         (
  86             Time::controlDictName,
  87             mesh.time().rootPath(),
  88             mesh.time().caseName().path()
  89         );
  90
  91         IOobject undecomposedBoundaryHeader
  92         (
  93             "boundary",
  94             undecomposedTime.constant(),
  95             polyMesh::meshSubDir,
  96             undecomposedTime,
  97             IOobject::MUST_READ,
  98             IOobject::NO_WRITE,
  99             false
 100         );
 101
 102         if (undecomposedBoundaryHeader.headerOk())
 103         {
 104             polyBoundaryMeshEntries undecomposedPatchEntries
 105             (
 106                 undecomposedBoundaryHeader
 107             );
 108
 109             forAll(undecomposedPatchEntries, patchi)
 110             {
 111                 patchNames.append
 112                 (
 113                     undecomposedPatchEntries[patchi].keyword()
 114                 );
 115
 116                 patchDictionary.add
 117                 (
 118                     undecomposedPatchEntries[patchi]
 119                 );
 120             }
 121         }
 122         else
 123         {
 124             FatalErrorIn ("referredCellList.C")
 125                 << nl << "unable to read undecomposed boundary file from "
 126                 << "constant/polyMesh" << nl
 127                 << abort(FatalError);
 128         }
 129
 130         labelIOList faceProcAddressing
 131         (
 132             IOobject
 133             (
 134                 "faceProcAddressing",
 135                 mesh.time().constant(),
 136                 polyMesh::meshSubDir,
 137                 mesh,
 138                 IOobject::MUST_READ,
 139                 IOobject::NO_WRITE,
 140                 false
 141             )
 142         );
 143
 144         labelList procPatches(mesh.globalData().processorPatches());
 145
 146         nUndecomposedPatches = patchNames.size();
 147
 148         // processorPatchSegmentMapping works by mapping the original patch and
 149         // half that a face on a processor patch was on before decomposition.
 150         // This creates a patch segment for each half of each original (cyclic)
 151         // patch which can be assessed separately.  There are n =
 152         // patchNames.size() original patches, k = 0 to n-1.  The mapping is:
 153         // value = 0: originally an internal face value = k, was originally on
 154         // the on the patch k-1, 1st half value = -k, was originally on the on
 155         // the patch k-1, 2nd half
 156
 157         forAll(procPatches,pP)
 158         {
 159             const processorPolyPatch& patch = refCast<const processorPolyPatch>
 160             (
 161                 mesh.boundaryMesh()[procPatches[pP]]
 162             );
 163
 164             labelList& procPatchSegMap = processorPatchSegmentMapping[pP];
 165
 166             procPatchSegMap.setSize(patch.size());
 167
 168             forAll (patch, pI)
 169             {
 170                 label decomposedMeshFace = patch.start() + pI;
 171
 172                 label faceProcAdd = faceProcAddressing[decomposedMeshFace];
 173
 174                 label globalFace = abs(faceProcAdd)-1;
 175
 176                 label minStart = -1;
 177
 178                 forAll(patchNames, pN)
 179                 {
 180                     if (patchDictionary.found(patchNames[pN]))
 181                     {
 182                         const dictionary& patchDict =
 183                         patchDictionary.subDict(patchNames[pN]);
 184
 185                         label startFace
 186                         (
 187                             readLabel
 188                             (
 189                                 patchDict.lookup("startFace")
 190                             )
 191                         );
 192
 193                         label nFaces(readLabel(patchDict.lookup("nFaces")));
 194
 195                         if (minStart < 0 || startFace < minStart)
 196                         {
 197                             minStart = startFace;
 198                         }
 199
 200                         if
 201                         (
 202                             globalFace >= startFace
 203                          && globalFace < startFace + nFaces/2
 204                         )
 205                         {
 206                             procPatchSegMap[pI] = pN + 1;
 207                         }
 208                         else if
 209                         (
 210                             globalFace >= startFace + nFaces/2
 211                          && globalFace < startFace + nFaces
 212                         )
 213                         {
 214                             procPatchSegMap[pI] = -(pN + 1);
 215                         }
 216                     }
 217                 }
 218
 219                 if (globalFace < minStart)
 220                 {
 221                     procPatchSegMap[pI] = 0;
 222                 }
 223             }
 224         }
 225
 226         forAll(procPatches,pP)
 227         {
 228             const processorPolyPatch& patch = refCast<const processorPolyPatch>
 229             (
 230                 mesh.boundaryMesh()[procPatches[pP]]
 231             );
 232
 233             {
 234                 OPstream toNeighbProc
 235                 (
 236                     Pstream::blocking,
 237                     patch.neighbProcNo()
 238                 );
 239
 240                 toNeighbProc << patch.faceCentres() << patch.faceAreas();
 241             }
 242         }
 243
 244         forAll(procPatches,pP)
 245         {
 246             const processorPolyPatch& patch = refCast<const processorPolyPatch>
 247             (
 248                 mesh.boundaryMesh()[procPatches[pP]]
 249             );
 250
 251             vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
 252
 253             neighbFaceCentres.setSize(patch.size());
 254
 255             vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
 256
 257             neighbFaceAreas.setSize(patch.size());
 258
 259             {
 260                 IPstream fromNeighbProc
 261                 (
 262                     Pstream::blocking,
 263                     patch.neighbProcNo()
 264                 );
 265
 266                 fromNeighbProc >> neighbFaceCentres >> neighbFaceAreas;
 267             }
 268         }
 269
 270         // *************************************************************
 271         // Tests that all processor patch segments from different
 272         // original patches prior to decomposition produce the same
 273         // transform. Check before 1st iteration.
 274         // *************************************************************
 275
 276         forAll(procPatches,pP)
 277         {
 278             const processorPolyPatch& patch = refCast<const processorPolyPatch>
 279             (
 280                 mesh.boundaryMesh()[procPatches[pP]]
 281             );
 282
 283             const vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
 284
 285             const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
 286
 287             label nUP;
 288
 289             for
 290             (
 291                 nUP = -nUndecomposedPatches;
 292                 nUP <= nUndecomposedPatches;
 293                 nUP++
 294             )
 295             {
 296                 DynamicList<vector> refOff;
 297
 298                 DynamicList<tensor> refTrans;
 299
 300                 forAll (patch, faceI)
 301                 {
 302                     if (processorPatchSegmentMapping[pP][faceI] == nUP)
 303                     {
 304                         referredCell testRefCell
 305                         (
 306                             mesh,
 307                             -1,
 308                             -1,
 309                             patch.faceCentres()[faceI],
 310                             neighbFaceCentres[faceI],
 311                             patch.faceNormals()[faceI],
 312                             neighbFaceAreas[faceI]
 313                             /(mag(neighbFaceAreas[faceI]) + VSMALL)
 314                         );
 315
 316                         refOff.append(testRefCell.offset());
 317
 318                         refTrans.append(testRefCell.rotation());
 319                     }
 320                 }
 321
 322                 refOff.shrink();
 323
 324                 refTrans.shrink();
 325
 326                 if (refOff.size())
 327                 {
 328                     if
 329                     (
 330                         sum(mag(refOff-refOff[0]))/refOff.size()
 331                             > interactionLists::transTol
 332                      || sum(mag(refTrans-refTrans[0]))/refTrans.size()
 333                             > interactionLists::transTol
 334                     )
 335                     {
 336                         FatalErrorIn ("referredCellList.C")
 337                             << nl << "Face pairs on patch "
 338                             << patch.name()
 339                             << ", segment " << patchNames[nUP]
 340                             << " do not give the same referring "
 341                             << " transformations to within tolerance of "
 342                             << interactionLists::transTol << nl
 343                             << " Referring offsets:" << refOff << nl
 344                             << " Average sum of mag difference: "
 345                             << sum(mag(refOff-refOff[0]))/refOff.size() << nl
 346                             << " Referring transforms:" << refTrans << nl
 347                             << " Average sum of mag difference: "
 348                             << sum(mag(refTrans-refTrans[0]))/refTrans.size()
 349                             << nl << abort(FatalError);
 350                     }
 351                 }
 352             }
 353         }
 354     }
 355
 356     label cellsReferredThisIteration = 1;
 357
 358     label iterationNo = 0;
 359
 360     while (cellsReferredThisIteration)
 361     {
 362         label refIntListStartSize = referredInteractionList.size();
 363
 364         forAll (mesh.boundaryMesh(), patchI)
 365         {
 366             // Treat local cyclics on each processor before processor
 367             // boundaries.  Separate treatment allows the serial version to run
 368             // transparently.
 369
 370             if (mesh.boundaryMesh()[patchI].type() == "cyclic")
 371             {
 372                 const cyclicPolyPatch& patch = refCast<const cyclicPolyPatch>
 373                 (
 374                     mesh.boundaryMesh()[patchI]
 375                 );
 376
 377                 if (patch.size())
 378                 {
 379                     if (iterationNo == 0)
 380                     {
 381                         // Tests that all combinations of face pairs produce the
 382                         // same transform.  Only check on 1st iteration.
 383
 384                         label faceL;
 385                         // A face in the 1st half of the patch
 386
 387                         label faceM;
 388                         // Face corresponding to faceL in the 2nd half of the
 389                         // patch. Assumes correct face ordering.
 390
 391                         vectorList refOff(patch.size()/2);
 392
 393                         List<tensor> refTrans(patch.size()/2);
 394
 395                         for
 396                         (
 397                             faceL = 0, faceM = patch.size()/2;
 398                             faceL < patch.size()/2;
 399                             faceL++, faceM++
 400                         )
 401                         {
 402                             referredCell testRefCell
 403                             (
 404                                 mesh,
 405                                 -1,
 406                                 -1,
 407                                 patch.faceCentres()[faceL],
 408                                 patch.faceCentres()[faceM],
 409                                 patch.faceNormals()[faceL],
 410                                 patch.faceNormals()[faceM]
 411                             );
 412
 413                             refOff[faceL] = testRefCell.offset();
 414
 415                             refTrans[faceL] = testRefCell.rotation();
 416                         }
 417
 418                         if
 419                         (
 420                             sum(mag(refOff - refOff[0]))/(patch.size()/2)
 421                                 > interactionLists::transTol
 422                          || sum(mag(refTrans - refTrans[0]))/(patch.size()/2)
 423                                 > interactionLists::transTol
 424                         )
 425                         {
 426                             FatalErrorIn ("referredCellList.C")
 427                                 << nl << "Face pairs on patch "
 428                                 << patch.name()
 429                                 << " do not give the same referring "
 430                                 << " transformations to within tolerance of "
 431                                 << interactionLists::transTol << nl
 432                                 << " Referring offsets:" << refOff << nl
 433                                 << " Average sum of mag difference: "
 434                                 << sum(mag(refOff-refOff[0]))/refOff.size()
 435                                 << nl
 436                                 << " Referring transforms:" << refTrans << nl
 437                                 << " Average sum of mag difference: "
 438                                 << sum(mag(refTrans-refTrans[0]))
 439                                    /refTrans.size()
 440                                 << nl << abort(FatalError);
 441                         }
 442                     }
 443
 444                     // *********************************************************
 445                     // 1st half of face list - 1st side of boundary
 446                     // *********************************************************
 447
 448                     label faceI;
 449
 450                     DynamicList<label> meshFacesOnThisSegment;
 451
 452                     for (faceI = 0; faceI < patch.size()/2; faceI++)
 453                     {
 454                         // unable to use the normal accessors for the polyPatch
 455                         // because points on separate halves need used
 456                         // separately
 457
 458                         meshFacesOnThisSegment.append(faceI + patch.start());
 459                     }
 460
 461                     meshFacesOnThisSegment.shrink();
 462
 463                     DynamicList<label> meshEdgesOnThisSegment;
 464
 465                     DynamicList<label> meshPointsOnThisSegment;
 466
 467                     forAll(meshFacesOnThisSegment, mFOTS)
 468                     {
 469                         const label segFace = meshFacesOnThisSegment[mFOTS];
 470
 471                         const labelList& faceEdges = mesh.faceEdges()[segFace];
 472
 473                         forAll (faceEdges, fE)
 474                         {
 475                             const label faceEdge(faceEdges[fE]);
 476
 477                             if
 478                             (
 479                                 findIndex
 480                                 (
 481                                     meshEdgesOnThisSegment,
 482                                     faceEdge
 483                                 ) == -1
 484                             )
 485                             {
 486                                 meshEdgesOnThisSegment.append(faceEdge);
 487                             }
 488                         }
 489
 490                         const face& facePoints(mesh.faces()[segFace]);
 491
 492                         forAll (facePoints, fP)
 493                         {
 494                             const label facePoint(facePoints[fP]);
 495
 496                             if
 497                             (
 498                                 findIndex
 499                                 (
 500                                     meshPointsOnThisSegment,
 501                                     facePoint
 502                                 )
 503                              ==
 504                                 -1
 505                             )
 506                             {
 507                                 meshPointsOnThisSegment.append(facePoint);
 508                             }
 509                         }
 510                     }
 511
 512                     meshEdgesOnThisSegment.shrink();
 513
 514                     meshPointsOnThisSegment.shrink();
 515
 516                     if (iterationNo == 0)
 517                     {
 518                         // Assessing real cells in range is only required on
 519                         // the 1st iteration because they do not change from
 520                         // iteration to iteration.
 521
 522                         labelList realCellsFoundInRange
 523                         (
 524                             il_.realCellsInRangeOfSegment
 525                             (
 526                                 meshFacesOnThisSegment,
 527                                 meshEdgesOnThisSegment,
 528                                 meshPointsOnThisSegment
 529                             )
 530                         );
 531
 532                         forAll(realCellsFoundInRange,cFIR)
 533                         {
 534                             const label realCell = realCellsFoundInRange[cFIR];
 535
 536                             referredCell cellToRefer
 537                             (
 538                                 mesh,
 539                                 Pstream::myProcNo(),
 540                                 realCell,
 541                                 patch.faceCentres()[0],
 542                                 patch.faceCentres()[patch.size()/2],
 543                                 patch.faceNormals()[0],
 544                                 patch.faceNormals()[patch.size()/2]
 545                             );
 546
 547                             // Test all existing referred and real cells to
 548                             // check duplicates are not being made or cells
 549                             // aren't being referred back onto themselves
 550
 551                             bool addCellToRefer = true;
 552
 553                             // Check if cellToRefer is an existing referred cell
 554
 555                             forAll(referredInteractionList, rIL)
 556                             {
 557                                 if
 558                                 (
 559                                     cellToRefer.duplicate
 560                                     (
 561                                         referredInteractionList[rIL]
 562                                     )
 563                                 )
 564                                 {
 565                                     addCellToRefer = false;
 566
 567                                     break;
 568                                 }
 569                             }
 570
 571                             // Check for cellToRefer being referred back
 572                             // ontop of a real cell
 573
 574                             if
 575                             (
 576                                 cellToRefer.duplicate
 577                                 (
 578                                     Pstream::myProcNo(),
 579                                     mesh.nCells()
 580                                 )
 581                             )
 582                             {
 583                                 addCellToRefer = false;
 584                             }
 585
 586                             if (addCellToRefer)
 587                             {
 588                                 referredInteractionList.append(cellToRefer);
 589                             }
 590
 591                             // add real cells found in range of cyclic patch
 592                             // to whole mesh list
 593
 594                             if (findIndex (rCellsWRRP, realCell) == -1)
 595                             {
 596                                 rCellsWRRP.append(realCell);
 597                             }
 598                         }
 599                     }
 600
 601                     referredInteractionList.shrink();
 602
 603                     labelList referredCellsFoundInRange
 604                     (
 605                         il_.referredCellsInRangeOfSegment
 606                         (
 607                             referredInteractionList,
 608                             meshFacesOnThisSegment,
 609                             meshEdgesOnThisSegment,
 610                             meshPointsOnThisSegment
 611                         )
 612                     );
 613
 614                     forAll(referredCellsFoundInRange,cFIR)
 615                     {
 616                         referredCell& existingRefCell =
 617                             referredInteractionList
 618                             [
 619                                 referredCellsFoundInRange[cFIR]
 620                             ];
 621
 622                         referredCell cellToReRefer =
 623                             existingRefCell.reRefer
 624                             (
 625                                 patch.faceCentres()[0],
 626                                 patch.faceCentres()[patch.size()/2],
 627                                 patch.faceNormals()[0],
 628                                 patch.faceNormals()[patch.size()/2]
 629                             );
 630
 631                         // Test all existing referred and real cells to check
 632                         // duplicates are not being made or cells aren't being
 633                         // referred back onto themselves
 634
 635                         bool addCellToReRefer = true;
 636
 637                         // Check if cellToRefer is an existing referred cell
 638
 639                         forAll(referredInteractionList, rIL)
 640                         {
 641                             if
 642                             (
 643                                 cellToReRefer.duplicate
 644                                 (
 645                                     referredInteractionList[rIL]
 646                                 )
 647                             )
 648                             {
 649                                 addCellToReRefer = false;
 650
 651                                 break;
 652                             }
 653                         }
 654
 655                         // Check for cellToRefer being referred back
 656                         // ontop of a real cell
 657
 658                         if
 659                         (
 660                             cellToReRefer.duplicate
 661                             (
 662                                 Pstream::myProcNo(),
 663                                 mesh.nCells()
 664                             )
 665                         )
 666                         {
 667                             addCellToReRefer = false;
 668                         }
 669
 670                         if (addCellToReRefer)
 671                         {
 672                             referredInteractionList.append(cellToReRefer);
 673                         }
 674                     }
 675
 676                     // *********************************************************
 677                     // 2nd half of face list - 2nd side of boundary
 678                     // *********************************************************
 679
 680                     meshFacesOnThisSegment.clear();
 681
 682                     for (faceI = patch.size()/2; faceI < patch.size(); faceI++)
 683                     {
 684                         // unable to use the normal accessors for the polyPatch
 685                         // because points on separate halves need used
 686                         // separately
 687
 688                         meshFacesOnThisSegment.append(faceI + patch.start());
 689                     }
 690
 691                     meshFacesOnThisSegment.shrink();
 692
 693                     meshEdgesOnThisSegment.clear();
 694
 695                     meshPointsOnThisSegment.clear();
 696
 697                     forAll(meshFacesOnThisSegment, mFOTS)
 698                     {
 699                         const label segFace = meshFacesOnThisSegment[mFOTS];
 700
 701                         const labelList& faceEdges = mesh.faceEdges()[segFace];
 702
 703                         forAll (faceEdges, fE)
 704                         {
 705                             const label faceEdge(faceEdges[fE]);
 706
 707                             if
 708                             (
 709                                 findIndex
 710                                 (
 711                                     meshEdgesOnThisSegment,
 712                                     faceEdge
 713                                 )
 714                              ==
 715                                 -1
 716                             )
 717                             {
 718                                 meshEdgesOnThisSegment.append(faceEdge);
 719                             }
 720                         }
 721
 722                         const face& facePoints(mesh.faces()[segFace]);
 723
 724                         forAll (facePoints, fP)
 725                         {
 726                             const label facePoint(facePoints[fP]);
 727
 728                             if
 729                             (
 730                                 findIndex
 731                                 (
 732                                     meshPointsOnThisSegment,
 733                                     facePoint
 734                                 )
 735                              ==
 736                                 -1
 737                             )
 738                             {
 739                                 meshPointsOnThisSegment.append(facePoint);
 740                             }
 741                         }
 742                     }
 743
 744                     meshEdgesOnThisSegment.shrink();
 745
 746                     meshPointsOnThisSegment.shrink();
 747
 748                     if (iterationNo == 0)
 749                     {
 750                         // Assessing real cells in range is only required on
 751                         // the 1st iteration because they do not change from
 752                         // iteration to iteration.
 753
 754                         labelList realCellsFoundInRange
 755                         (
 756                             il_.realCellsInRangeOfSegment
 757                             (
 758                                 meshFacesOnThisSegment,
 759                                 meshEdgesOnThisSegment,
 760                                 meshPointsOnThisSegment
 761                             )
 762                         );
 763
 764                         forAll(realCellsFoundInRange,cFIR)
 765                         {
 766                             const label realCell = realCellsFoundInRange[cFIR];
 767
 768                             referredCell cellToRefer
 769                             (
 770                                 mesh,
 771                                 Pstream::myProcNo(),
 772                                 realCell,
 773                                 patch.faceCentres()[patch.size()/2],
 774                                 patch.faceCentres()[0],
 775                                 patch.faceNormals()[patch.size()/2],
 776                                 patch.faceNormals()[0]
 777                             );
 778
 779                             // Test all existing referred and real cells to
 780                             // check duplicates are not being made or cells
 781                             // aren't being referred back onto themselves
 782
 783                             bool addCellToRefer = true;
 784
 785                             // Check if cellToRefer is an existing referred cell
 786
 787                             forAll(referredInteractionList, rIL)
 788                             {
 789                                 if
 790                                 (
 791                                     cellToRefer.duplicate
 792                                     (
 793                                         referredInteractionList[rIL]
 794                                     )
 795                                 )
 796                                 {
 797                                     addCellToRefer = false;
 798
 799                                     break;
 800                                 }
 801                             }
 802
 803                             // Check for cellToRefer being referred back
 804                             // ontop of a real cell
 805
 806                             if
 807                             (
 808                                 cellToRefer.duplicate
 809                                 (
 810                                     Pstream::myProcNo(),
 811                                     mesh.nCells()
 812                                 )
 813                             )
 814                             {
 815                                 addCellToRefer = false;
 816                             }
 817
 818                             if (addCellToRefer)
 819                             {
 820                                 referredInteractionList.append(cellToRefer);
 821                             }
 822
 823                             // add real cells found in range of cyclic patch
 824                             // to whole mesh list
 825
 826                             if (findIndex (rCellsWRRP, realCell) == -1)
 827                             {
 828                                 rCellsWRRP.append(realCell);
 829                             }
 830                         }
 831                     }
 832
 833                     referredInteractionList.shrink();
 834
 835                     referredCellsFoundInRange =
 836                         il_.referredCellsInRangeOfSegment
 837                         (
 838                             referredInteractionList,
 839                             meshFacesOnThisSegment,
 840                             meshEdgesOnThisSegment,
 841                             meshPointsOnThisSegment
 842                         );
 843
 844                     forAll(referredCellsFoundInRange,cFIR)
 845                     {
 846                         referredCell& existingRefCell =
 847                             referredInteractionList
 848                             [
 849                                 referredCellsFoundInRange[cFIR]
 850                             ];
 851
 852                         referredCell cellToReRefer =
 853                             existingRefCell.reRefer
 854                             (
 855                                 patch.faceCentres()[patch.size()/2],
 856                                 patch.faceCentres()[0],
 857                                 patch.faceNormals()[patch.size()/2],
 858                                 patch.faceNormals()[0]
 859                             );
 860
 861                         // Test all existing referred and real cells to check
 862                         // duplicates are not being made or cells aren't being
 863                         // referred back onto themselves
 864
 865                         bool addCellToReRefer = true;
 866
 867                         // Check if cellToRefer is an existing referred cell
 868
 869                         forAll(referredInteractionList, rIL)
 870                         {
 871                             if
 872                             (
 873                                 cellToReRefer.duplicate
 874                                 (
 875                                     referredInteractionList[rIL]
 876                                 )
 877                             )
 878                             {
 879                                 addCellToReRefer = false;
 880
 881                                 break;
 882                             }
 883                         }
 884
 885                         // Check for cellToRefer being referred back
 886                         // ontop of a real cell
 887
 888                         if
 889                         (
 890                             cellToReRefer.duplicate
 891                             (
 892                                 Pstream::myProcNo(),
 893                                 mesh.nCells()
 894                             )
 895                         )
 896                         {
 897                             addCellToReRefer = false;
 898                         }
 899
 900                         if (addCellToReRefer)
 901                         {
 902                             referredInteractionList.append(cellToReRefer);
 903                         }
 904                     }
 905                 }
 906             }
 907         }
 908
 909         if (Pstream::parRun())
 910         {
 911             labelList procPatches(mesh.globalData().processorPatches());
 912
 913             forAll(procPatches,pP)
 914             {
 915                 const processorPolyPatch& patch =
 916                     refCast<const processorPolyPatch>
 917                     (
 918                         mesh.boundaryMesh()[procPatches[pP]]
 919                     );
 920
 921                 DynamicList<referredCell> referredCellsToTransfer;
 922
 923                 const vectorList& neighbFaceCentres =
 924                     allNeighbourFaceCentres[pP];
 925
 926                 const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
 927
 928                 label nUP;
 929
 930                 for
 931                 (
 932                     nUP = -nUndecomposedPatches;
 933                     nUP <= nUndecomposedPatches;
 934                     nUP++
 935                 )
 936                 {
 937                     // faceT is used to specify one face on this patch segment
 938                     // that will be used to calculate the transformation values.
 939                     // All faces are guaranteed to produce the same transform
 940                     // because of the checks carried out at the start of the
 941                     // function.  Setting to -1 until the 1st face on this
 942                     // segment is found.
 943
 944                     label faceT = -1;
 945
 946                     DynamicList<label> meshFacesOnThisSegment;
 947
 948                     forAll (patch, faceI)
 949                     {
 950                         if (processorPatchSegmentMapping[pP][faceI] == nUP)
 951                         {
 952                             if (faceT == -1)
 953                             {
 954                                 faceT = faceI;
 955                             }
 956
 957                             meshFacesOnThisSegment.append
 958                             (
 959                                 faceI + patch.start()
 960                             );
 961                         }
 962                     }
 963
 964                     meshFacesOnThisSegment.shrink();
 965
 966                     DynamicList<label> meshEdgesOnThisSegment;
 967
 968                     DynamicList<label> meshPointsOnThisSegment;
 969
 970                     forAll(meshFacesOnThisSegment, mFOTS)
 971                     {
 972                         const label segFace = meshFacesOnThisSegment[mFOTS];
 973
 974                         const labelList& faceEdges = mesh.faceEdges()[segFace];
 975
 976                         forAll (faceEdges, fE)
 977                         {
 978                             const label faceEdge(faceEdges[fE]);
 979
 980                             if
 981                             (
 982                                 findIndex
 983                                 (
 984                                     meshEdgesOnThisSegment,
 985                                     faceEdge
 986                                 )
 987                              ==
 988                                 -1
 989                             )
 990                             {
 991                                 meshEdgesOnThisSegment.append(faceEdge);
 992                             }
 993                         }
 994
 995                         const face& facePoints(mesh.faces()[segFace]);
 996
 997                         forAll (facePoints, fP)
 998                         {
 999                             const label facePoint(facePoints[fP]);
1000
1001                             if
1002                             (
1003                                 findIndex
1004                                 (
1005                                     meshPointsOnThisSegment,
1006                                     facePoint
1007                                 )
1008                              ==
1009                                 -1
1010                             )
1011                             {
1012                                 meshPointsOnThisSegment.append(facePoint);
1013                             }
1014                         }
1015                     }
1016
1017                     meshEdgesOnThisSegment.shrink();
1018
1019                     meshPointsOnThisSegment.shrink();
1020
1021                     if (meshFacesOnThisSegment.size())
1022                     {
1023                         if (faceT == -1)
1024                         {
1025                             FatalErrorIn ("referredCellList.C")
1026                                 << nl << "faceT == -1 encountered but "
1027                                 << meshFacesOnThisSegment.size()
1028                                 << " faces found on patch segment."
1029                                 << abort(FatalError);
1030                         }
1031
1032                         if (iterationNo == 0)
1033                         {
1034                             // Assessing real cells in range is only required on
1035                             // the 1st iteration because they do not change from
1036                             // iteration to iteration.
1037
1038                             labelList realCellsFoundInRange
1039                             (
1040                                 il_.realCellsInRangeOfSegment
1041                                 (
1042                                     meshFacesOnThisSegment,
1043                                     meshEdgesOnThisSegment,
1044                                     meshPointsOnThisSegment
1045                                 )
1046                             );
1047
1048                             forAll(realCellsFoundInRange,cFIR)
1049                             {
1050                                 const label realCell =
1051                                     realCellsFoundInRange[cFIR];
1052
1053                                 referredCell cellToRefer
1054                                 (
1055                                     mesh,
1056                                     Pstream::myProcNo(),
1057                                     realCell,
1058                                     patch.faceCentres()[faceT],
1059                                     neighbFaceCentres[faceT],
1060                                     patch.faceNormals()[faceT],
1061                                     neighbFaceAreas[faceT]
1062                                     /(mag(neighbFaceAreas[faceT]) + VSMALL)
1063                                 );
1064
1065                                 referredCellsToTransfer.append(cellToRefer);
1066
1067                                 // add real cells found in range of processor
1068                                 // patch to whole mesh list
1069
1070                                 if (findIndex (rCellsWRRP, realCell) == -1)
1071                                 {
1072                                     rCellsWRRP.append(realCell);
1073                                 }
1074                             }
1075                         }
1076
1077                         referredInteractionList.shrink();
1078
1079                         labelList referredCellsFoundInRange
1080                         (
1081                             il_.referredCellsInRangeOfSegment
1082                             (
1083                                 referredInteractionList,
1084                                 meshFacesOnThisSegment,
1085                                 meshEdgesOnThisSegment,
1086                                 meshPointsOnThisSegment
1087                             )
1088                         );
1089
1090                         forAll(referredCellsFoundInRange,cFIR)
1091                         {
1092                             referredCell& existingRefCell =
1093                                 referredInteractionList
1094                                 [
1095                                     referredCellsFoundInRange[cFIR]
1096                                 ];
1097
1098                             referredCell cellToReRefer =
1099                                 existingRefCell.reRefer
1100                                 (
1101                                     patch.faceCentres()[faceT],
1102                                     neighbFaceCentres[faceT],
1103                                     patch.faceNormals()[faceT],
1104                                     neighbFaceAreas[faceT]
1105                                     /(mag(neighbFaceAreas[faceT]) + VSMALL)
1106                                 );
1107
1108                             referredCellsToTransfer.append(cellToReRefer);
1109                         }
1110                     }
1111                 }
1112
1113                 referredCellsToTransfer.shrink();
1114
1115                 // Send these cells to the neighbouring processor.
1116
1117                 {
1118                     OPstream toNeighbProc
1119                     (
1120                         Pstream::blocking,
1121                         patch.neighbProcNo()
1122                     );
1123
1124                     toNeighbProc << referredCellsToTransfer;
1125                 }
1126             }
1127
1128             forAll(procPatches,pP)
1129             {
1130                 const processorPolyPatch& patch =
1131                 refCast<const processorPolyPatch>
1132                 (
1133                     mesh.boundaryMesh()[procPatches[pP]]
1134                 );
1135
1136                 // Receive the cells from neighbour
1137
1138                 List<referredCell> referredCellsFromNeighbour(patch.size());
1139
1140                 {
1141                     IPstream fromNeighbProc
1142                     (
1143                         Pstream::blocking,
1144                         patch.neighbProcNo()
1145                     );
1146
1147                     fromNeighbProc >> referredCellsFromNeighbour;
1148                 }
1149
1150                 // Check to see if they are duplicates, if not append
1151                 // them to the referredInteractionList
1152
1153                 forAll(referredCellsFromNeighbour,rCFN)
1154                 {
1155                     referredCell& cellToRefer =
1156                     referredCellsFromNeighbour[rCFN];
1157
1158                     // Test all existing referred and real cells to check
1159                     // duplicates are not being made or cells aren't being
1160                     // referred back onto themselves
1161
1162                     bool addCellToRefer = true;
1163
1164                     // Check if cellToRefer is an existing referred cell
1165
1166                     forAll(referredInteractionList, rIL)
1167                     {
1168                         if (cellToRefer.duplicate(referredInteractionList[rIL]))
1169                         {
1170                             addCellToRefer = false;
1171
1172                             break;
1173                         }
1174                     }
1175
1176                     // Check for cellToRefer being referred back ontop of a real
1177                     // cell
1178
1179                     if
1180                     (
1181                         cellToRefer.duplicate
1182                         (
1183                             Pstream::myProcNo(),
1184                             mesh.nCells()
1185                         )
1186                     )
1187                     {
1188                         addCellToRefer = false;
1189                     }
1190
1191                     if (addCellToRefer)
1192                     {
1193                         referredInteractionList.append(cellToRefer);
1194                     }
1195                 }
1196             }
1197         }
1198
1199         if (iterationNo == 0)
1200         {
1201             // record all real cells in range of any referring patch (cyclic or
1202             // processor) on the first iteration when the real cells are
1203             // evaluated.
1204
1205             rCellsWRRP.shrink();
1206
1207             // construct {faces, edges, points}WithinRCutMaxOfAnyReferringPatch
1208
1209             forAll(rCellsWRRP, rCWR)
1210             {
1211                 const label realCell(rCellsWRRP[rCWR]);
1212
1213                 const labelList& rCFaces
1214                 (
1215                     mesh.cells()[realCell]
1216                 );
1217
1218                 forAll(rCFaces, rCF)
1219                 {
1220                     const label f(rCFaces[rCF]);
1221
1222                     if (findIndex(rFacesWRRP,f) == -1)
1223                     {
1224                         rFacesWRRP.append(f);
1225                     }
1226                 }
1227
1228                 const labelList& rCEdges
1229                 (
1230                     mesh.cellEdges()[realCell]
1231                 );
1232
1233                 forAll(rCEdges, rCE)
1234                 {
1235                     const label e(rCEdges[rCE]);
1236
1237                     if (findIndex(rEdgesWRRP,e) == -1)
1238                     {
1239                         rEdgesWRRP.append(e);
1240                     }
1241                 }
1242
1243                 const labelList& rCPoints
1244                 (
1245                     mesh.cellPoints()[realCell]
1246                 );
1247
1248                 forAll(rCPoints, rCP)
1249                 {
1250                     const label p(rCPoints[rCP]);
1251
1252                     if (findIndex(rPointsWRRP,p) == -1)
1253                     {
1254                         rPointsWRRP.append(p);
1255                     }
1256                 }
1257             }
1258
1259             rFacesWRRP.shrink();
1260
1261             rEdgesWRRP.shrink();
1262
1263             rPointsWRRP.shrink();
1264         }
1265
1266         iterationNo++;
1267
1268         cellsReferredThisIteration =
1269             referredInteractionList.size() - refIntListStartSize;
1270
1271         reduce(cellsReferredThisIteration, sumOp<label>());
1272
1273         Info<< tab << "Cells added this iteration: "
1274             << cellsReferredThisIteration << endl;
1275     }
1276
1277     referredInteractionList.shrink();
1278
1279     //     Info<< "referredInteractionList.size() = "
1280     //         << referredInteractionList.size() << endl;
1281
1282     //     forAll(referredInteractionList,rIL)
1283     //     {
1284     //         Info<< referredInteractionList[rIL];
1285     //     }
1286
1287     (*this).setSize
1288     (
1289         referredInteractionList.size()
1290     );
1291
1292     forAll(referredInteractionList, rIL)
1293     {
1294         (*this)[rIL] = referredInteractionList[rIL];
1295     }
1296
1297     Info<< nl << "Finding real cells in range of referred cells" << endl;
1298
1299     forAll(*this, rC)
1300     {
1301         const polyMesh& mesh(il_.mesh());
1302
1303         referredCell& refCell = (*this)[rC];
1304
1305         DynamicList<label> realCellsFoundInRange;
1306
1307         const vectorList& refCellPoints = refCell.vertexPositions();
1308
1309         forAll(rFacesWRRP, rCF)
1310         {
1311             const label f(rFacesWRRP[rCF]);
1312
1313             if (il_.testPointFaceDistance(refCellPoints,f))
1314             {
1315                 const label cellO(mesh.faceOwner()[f]);
1316
1317                 if (findIndex(realCellsFoundInRange, cellO) == -1)
1318                 {
1319                     realCellsFoundInRange.append(cellO);
1320                 }
1321
1322                 if (mesh.isInternalFace(f))
1323                 {
1324                     // boundary faces will not have neighbour information
1325
1326                     const label cellN(mesh.faceNeighbour()[f]);
1327
1328                     if (findIndex(realCellsFoundInRange, cellN) == -1)
1329                     {
1330                         realCellsFoundInRange.append(cellN);
1331                     }
1332                 }
1333             }
1334         }
1335
1336         forAll(rPointsWRRP, rCP)
1337         {
1338             const label p(rPointsWRRP[rCP]);
1339
1340             if (il_.testPointFaceDistance(p,refCell))
1341             {
1342                 const labelList& pCells(mesh.pointCells()[p]);
1343
1344                 forAll(pCells, pC)
1345                 {
1346                     const label cellI(pCells[pC]);
1347
1348                     if (findIndex(realCellsFoundInRange, cellI) == -1)
1349                     {
1350                         realCellsFoundInRange.append(cellI);
1351                     }
1352                 }
1353             }
1354         }
1355
1356
1357         const edgeList& refCellEdges = refCell.edges();
1358
1359         forAll(rEdgesWRRP, rCE)
1360         {
1361             const label edgeIIndex(rEdgesWRRP[rCE]);
1362
1363             const edge& eI(mesh.edges()[edgeIIndex]);
1364
1365             forAll(refCellEdges, rCE)
1366             {
1367                 const edge& eJ(refCellEdges[rCE]);
1368
1369                 if
1370                 (
1371                     il_.testEdgeEdgeDistance
1372                     (
1373                         eI,
1374                         refCellPoints[eJ.start()],
1375                         refCellPoints[eJ.end()]
1376                     )
1377                 )
1378                 {
1379                     const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
1380
1381                     forAll(eICells, eIC)
1382                     {
1383                         const label cellI(eICells[eIC]);
1384
1385                         if (findIndex(realCellsFoundInRange, cellI) == -1)
1386                         {
1387                             realCellsFoundInRange.append(cellI);
1388                         }
1389                     }
1390                 }
1391             }
1392         }
1393
1394 //         scalar rCutMaxSqr = molCloud_.rCutMax()*molCloud_.rCutMax();
1395 //
1396 //         forAll (molCloud_.mesh().points(), pointIIndex)
1397 //         {
1398 //             const point& ptI
1399 //             (
1400 //                 molCloud_.mesh().points()[pointIIndex]
1401 //             );
1402 //
1403 //             forAll(refCellPoints, rCP)
1404 //             {
1405 //                 if (magSqr(ptI - refCellPoints[rCP]) <= rCutMaxSqr)
1406 //                 {
1407 //                     const labelList& ptICells
1408 //                     (
1409 //                         molCloud_.mesh().pointCells()[pointIIndex]
1410 //                     );
1411 //
1412 //                     forAll(ptICells, pIC)
1413 //                     {
1414 //                         const label cellI(ptICells[pIC]);
1415 //
1416 //                         if (findIndex(realCellsFoundInRange, cellI) == -1)
1417 //                         {
1418 //                             realCellsFoundInRange.append(cellI);
1419 //                         }
1420 //                     }
1421 //                 }
1422 //             }
1423 //         }
1424
1425         refCell.realCells() = realCellsFoundInRange.shrink();
1426     }
1427 }
1428
1429
1430 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1431
1432 Foam::referredCellList::referredCellList
1433 (
1434     interactionLists& il,
1435     bool pointPointListBuild
1436 )
1437 :
1438     List<referredCell>(),
1439     il_(il)
1440 {
1441     buildReferredCellList(pointPointListBuild);
1442 }
1443
1444
1445 Foam::referredCellList::referredCellList(interactionLists& il)
1446 :
1447     List<referredCell>(),
1448     il_(il)
1449 {
1450     Info<< "Read referredCellList from disk not implemented" << endl;
1451 }
1452
1453
1454 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1455
1456 Foam::referredCellList::~referredCellList()
1457 {}
1458
1459
1460 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1461
1462 void Foam::referredCellList::referMolecules
1463 (
1464     const List<DynamicList<molecule*> >& cellOccupancy
1465 )
1466 {
1467     // Create referred molecules for sending using cell occupancy and
1468     // cellSendingReferralLists
1469
1470     forAll(il_.cellSendingReferralLists(), cSRL)
1471     {
1472         const sendingReferralList& sRL
1473         (
1474             il_.cellSendingReferralLists()[cSRL]
1475         );
1476
1477         List<DynamicList<referredMolecule> > molsToReferOut(sRL.size());
1478
1479         forAll(sRL, sRLI)
1480         {
1481             List<molecule*> realMols = cellOccupancy[sRL[sRLI]];
1482
1483             forAll (realMols, rM)
1484             {
1485                 molecule* mol = realMols[rM];
1486
1487                 molsToReferOut[sRLI].append
1488                 (
1489                     referredMolecule
1490                     (
1491                         mol->id(),
1492                         mol->position(),
1493                         mol->sitePositions()
1494                     )
1495                 );
1496             }
1497
1498             molsToReferOut[sRLI].shrink();
1499         }
1500
1501         // Send lists of referred molecules to other processors
1502
1503         if (sRL.destinationProc() != Pstream::myProcNo())
1504         {
1505             if (Pstream::parRun())
1506             {
1507                 OPstream toInteractingProc
1508                 (
1509                     Pstream::blocking,
1510                     sRL.destinationProc()
1511                 );
1512
1513                 toInteractingProc << molsToReferOut;
1514             }
1515         }
1516         else
1517         {
1518             // Refer molecules directly for referred cells on the same
1519             // processor.
1520
1521             const receivingReferralList& rRL
1522             (
1523                 il_.cellReceivingReferralLists()[cSRL]
1524             );
1525
1526             forAll(rRL, rRLI)
1527             {
1528                 forAll(rRL[rRLI], rC)
1529                 {
1530                     referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
1531
1532                     refCellToRefMolsTo.referInMols(molsToReferOut[rRLI]);
1533                 }
1534             }
1535         }
1536     }
1537
1538     // Receive referred molecule lists to and distribute to referredCells
1539     // according tocellReceivingReferralLists, referredCells deal with the
1540     // transformations themselves
1541
1542     forAll(il_.cellReceivingReferralLists(), cRRL)
1543     {
1544         const receivingReferralList& rRL
1545         (
1546             il_.cellReceivingReferralLists()[cRRL]
1547         );
1548
1549         List<List<referredMolecule> > molsToReferIn(rRL.size());
1550
1551         if (rRL.sourceProc() != Pstream::myProcNo())
1552         {
1553             if (Pstream::parRun())
1554             {
1555                 IPstream fromInteractingProc
1556                 (
1557                     Pstream::blocking,
1558                     rRL.sourceProc()
1559                 );
1560
1561                 fromInteractingProc >> molsToReferIn;
1562             }
1563
1564             forAll(rRL, rRLI)
1565             {
1566                 forAll(rRL[rRLI], rC)
1567                 {
1568                     referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
1569
1570                     refCellToRefMolsTo.referInMols(molsToReferIn[rRLI]);
1571                 }
1572             }
1573         }
1574     }
1575 }
1576
1577
1578 // ************************************************************************* //