src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "hexRef8.H"
  28
  29 #include "polyMesh.H"
  30 #include "polyTopoChange.H"
  31 #include "meshTools.H"
  32 #include "polyAddFace.H"
  33 #include "polyAddPoint.H"
  34 #include "polyAddCell.H"
  35 #include "polyModifyFace.H"
  36 #include "syncTools.H"
  37 #include "faceSet.H"
  38 #include "cellSet.H"
  39 #include "pointSet.H"
  40 #include "OFstream.H"
  41 #include "Time.H"
  42 #include "FaceCellWave.H"
  43 #include "mapDistributePolyMesh.H"
  44 #include "refinementData.H"
  45 #include "refinementDistanceData.H"
  46
  47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  48
  49 namespace Foam
  50 {
  51     defineTypeNameAndDebug(hexRef8, 0);
  52
  53     //- Reduction class. If x and y are not equal assign value.
  54     template< int value >
  55     class ifEqEqOp
  56     {
  57         public:
  58         void operator()( label& x, const label& y ) const
  59         {
  60             x = (x==y) ? x:value;
  61         }
  62     };
  63 }
  64
  65
  66 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  67
  68 void Foam::hexRef8::reorder
  69 (
  70     const labelList& map,
  71     const label len,
  72     const label null,
  73     labelList& elems
  74 )
  75 {
  76     labelList newElems(len, null);
  77
  78     forAll(elems, i)
  79     {
  80         label newI = map[i];
  81
  82         if (newI >= len)
  83         {
  84             FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
  85         }
  86
  87         if (newI >= 0)
  88         {
  89             newElems[newI] = elems[i];
  90         }
  91     }
  92
  93     elems.transfer(newElems);
  94 }
  95
  96
  97 void Foam::hexRef8::getFaceInfo
  98 (
  99     const label faceI,
 100     label& patchID,
 101     label& zoneID,
 102     label& zoneFlip
 103 ) const
 104 {
 105     patchID = -1;
 106
 107     if (!mesh_.isInternalFace(faceI))
 108     {
 109         patchID = mesh_.boundaryMesh().whichPatch(faceI);
 110     }
 111
 112     zoneID = mesh_.faceZones().whichZone(faceI);
 113
 114     zoneFlip = false;
 115
 116     if (zoneID >= 0)
 117     {
 118         const faceZone& fZone = mesh_.faceZones()[zoneID];
 119
 120         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 121     }
 122 }
 123
 124
 125 // Adds a face on top of existing faceI.
 126 Foam::label Foam::hexRef8::addFace
 127 (
 128     polyTopoChange& meshMod,
 129     const label faceI,
 130     const face& newFace,
 131     const label own,
 132     const label nei
 133 ) const
 134 {
 135     label patchID, zoneID, zoneFlip;
 136
 137     getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 138
 139     label newFaceI = -1;
 140
 141     if ((nei == -1) || (own < nei))
 142     {
 143         // Ordering ok.
 144         newFaceI = meshMod.setAction
 145         (
 146             polyAddFace
 147             (
 148                 newFace,                    // face
 149                 own,                        // owner
 150                 nei,                        // neighbour
 151                 -1,                         // master point
 152                 -1,                         // master edge
 153                 faceI,                      // master face for addition
 154                 false,                      // flux flip
 155                 patchID,                    // patch for face
 156                 zoneID,                     // zone for face
 157                 zoneFlip                    // face zone flip
 158             )
 159         );
 160     }
 161     else
 162     {
 163         // Reverse owner/neighbour
 164         newFaceI = meshMod.setAction
 165         (
 166             polyAddFace
 167             (
 168                 newFace.reverseFace(),      // face
 169                 nei,                        // owner
 170                 own,                        // neighbour
 171                 -1,                         // master point
 172                 -1,                         // master edge
 173                 faceI,                      // master face for addition
 174                 false,                      // flux flip
 175                 patchID,                    // patch for face
 176                 zoneID,                     // zone for face
 177                 zoneFlip                    // face zone flip
 178             )
 179         );
 180     }
 181     return newFaceI;
 182 }
 183
 184
 185 // Adds an internal face from an edge. Assumes orientation correct.
 186 // Problem is that the face is between four new vertices. So what do we provide
 187 // as master? The only existing mesh item we have is the edge we have split.
 188 // Have to be careful in only using it if it has internal faces since otherwise
 189 // polyMeshMorph will complain (because it cannot generate a sensible mapping
 190 // for the face)
 191 Foam::label Foam::hexRef8::addInternalFace
 192 (
 193     polyTopoChange& meshMod,
 194     const label meshFaceI,
 195     const label meshPointI,
 196     const face& newFace,
 197     const label own,
 198     const label nei
 199 ) const
 200 {
 201     if (mesh_.isInternalFace(meshFaceI))
 202     {
 203         return meshMod.setAction
 204         (
 205             polyAddFace
 206             (
 207                 newFace,                    // face
 208                 own,                        // owner
 209                 nei,                        // neighbour
 210                 -1,                         // master point
 211                 -1,                         // master edge
 212                 meshFaceI,                  // master face for addition
 213                 false,                      // flux flip
 214                 -1,                         // patch for face
 215                 -1,                         // zone for face
 216                 false                       // face zone flip
 217             )
 218         );
 219     }
 220     else
 221     {
 222         // Two choices:
 223         // - append (i.e. create out of nothing - will not be mapped)
 224         //   problem: field does not get mapped.
 225         // - inflate from point.
 226         //   problem: does interpolative mapping which constructs full
 227         //   volPointInterpolation!
 228
 229         // For now create out of nothing
 230
 231         return meshMod.setAction
 232         (
 233             polyAddFace
 234             (
 235                 newFace,                    // face
 236                 own,                        // owner
 237                 nei,                        // neighbour
 238                 -1,                         // master point
 239                 -1,                         // master edge
 240                 -1,                         // master face for addition
 241                 false,                      // flux flip
 242                 -1,                         // patch for face
 243                 -1,                         // zone for face
 244                 false                       // face zone flip
 245             )
 246         );
 247
 248
 249         ////- Inflate-from-point:
 250         //// Check if point has any internal faces we can use.
 251         //label masterPointI = -1;
 252         //
 253         //const labelList& pFaces = mesh_.pointFaces()[meshPointI];
 254         //
 255         //forAll(pFaces, i)
 256         //{
 257         //    if (mesh_.isInternalFace(pFaces[i]))
 258         //    {
 259         //        // meshPoint uses internal faces so ok to inflate from it
 260         //        masterPointI = meshPointI;
 261         //
 262         //        break;
 263         //    }
 264         //}
 265         //
 266         //return meshMod.setAction
 267         //(
 268         //    polyAddFace
 269         //    (
 270         //        newFace,                    // face
 271         //        own,                        // owner
 272         //        nei,                        // neighbour
 273         //        masterPointI,               // master point
 274         //        -1,                         // master edge
 275         //        -1,                         // master face for addition
 276         //        false,                      // flux flip
 277         //        -1,                         // patch for face
 278         //        -1,                         // zone for face
 279         //        false                       // face zone flip
 280         //    )
 281         //);
 282     }
 283 }
 284
 285
 286 // Modifies existing faceI for either new owner/neighbour or new face points.
 287 void Foam::hexRef8::modFace
 288 (
 289     polyTopoChange& meshMod,
 290     const label faceI,
 291     const face& newFace,
 292     const label own,
 293     const label nei
 294 ) const
 295 {
 296     label patchID, zoneID, zoneFlip;
 297
 298     getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 299
 300     if
 301     (
 302         (own != mesh_.faceOwner()[faceI])
 303      || (
 304             mesh_.isInternalFace(faceI)
 305          && (nei != mesh_.faceNeighbour()[faceI])
 306         )
 307      || (newFace != mesh_.faces()[faceI])
 308     )
 309     {
 310         if ((nei == -1) || (own < nei))
 311         {
 312             meshMod.setAction
 313             (
 314                 polyModifyFace
 315                 (
 316                     newFace,            // modified face
 317                     faceI,              // label of face being modified
 318                     own,                // owner
 319                     nei,                // neighbour
 320                     false,              // face flip
 321                     patchID,            // patch for face
 322                     false,              // remove from zone
 323                     zoneID,             // zone for face
 324                     zoneFlip            // face flip in zone
 325                 )
 326             );
 327         }
 328         else
 329         {
 330             meshMod.setAction
 331             (
 332                 polyModifyFace
 333                 (
 334                     newFace.reverseFace(),  // modified face
 335                     faceI,                  // label of face being modified
 336                     nei,                    // owner
 337                     own,                    // neighbour
 338                     false,                  // face flip
 339                     patchID,                // patch for face
 340                     false,                  // remove from zone
 341                     zoneID,                 // zone for face
 342                     zoneFlip                // face flip in zone
 343                 )
 344             );
 345         }
 346     }
 347 }
 348
 349
 350 // Bit complex way to determine the unrefined edge length.
 351 Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
 352 {
 353     // Determine minimum edge length per refinement level
 354     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 355
 356     const scalar GREAT2 = sqr(GREAT);
 357
 358     label nLevels = gMax(cellLevel_)+1;
 359
 360     scalarField typEdgeLenSqr(nLevels, GREAT2);
 361
 362
 363     // 1. Look only at edges surrounded by cellLevel cells only.
 364     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 365
 366     {
 367         // Per edge the cellLevel of connected cells. -1 if not set,
 368         // labelMax if different levels, otherwise levels of connected cells.
 369         labelList edgeLevel(mesh_.nEdges(), -1);
 370
 371         forAll(cellLevel_, cellI)
 372         {
 373             const label cLevel = cellLevel_[cellI];
 374
 375             const labelList& cEdges = mesh_.cellEdges()[cellI];
 376
 377             forAll(cEdges, i)
 378             {
 379                 label edgeI = cEdges[i];
 380
 381                 if (edgeLevel[edgeI] == -1)
 382                 {
 383                     edgeLevel[edgeI] = cLevel;
 384                 }
 385                 else if (edgeLevel[edgeI] == labelMax)
 386                 {
 387                     // Already marked as on different cellLevels
 388                 }
 389                 else if (edgeLevel[edgeI] != cLevel)
 390                 {
 391                     edgeLevel[edgeI] = labelMax;
 392                 }
 393             }
 394         }
 395
 396         // Make sure that edges with different levels on different processors
 397         // are also marked. Do the same test (edgeLevel != cLevel) on coupled
 398         // edges.
 399         syncTools::syncEdgeList
 400         (
 401             mesh_,
 402             edgeLevel,
 403             ifEqEqOp<labelMax>(),
 404             labelMin,
 405             false               // no separation
 406         );
 407
 408         // Now use the edgeLevel with a valid value to determine the
 409         // length per level.
 410         forAll(edgeLevel, edgeI)
 411         {
 412             const label eLevel = edgeLevel[edgeI];
 413
 414             if (eLevel >= 0 && eLevel < labelMax)
 415             {
 416                 const edge& e = mesh_.edges()[edgeI];
 417
 418                 scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
 419
 420                 typEdgeLenSqr[eLevel] = min(typEdgeLenSqr[eLevel], edgeLenSqr);
 421             }
 422         }
 423     }
 424
 425     // Get the minimum per level over all processors. Note minimum so if
 426     // cells are not cubic we use the smallest edge side.
 427     Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>());
 428     Pstream::listCombineScatter(typEdgeLenSqr);
 429
 430     if (debug)
 431     {
 432         Pout<< "hexRef8::getLevel0EdgeLength() :"
 433             << " After phase1: Edgelengths (squared) per refinementlevel:"
 434             << typEdgeLenSqr << endl;
 435     }
 436
 437
 438     // 2. For any levels where we haven't determined a valid length yet
 439     //    use any surrounding cell level. Here we use the max so we don't
 440     //    pick up levels between celllevel and higher celllevel (will have
 441     //    edges sized according to highest celllevel)
 442     //    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 443
 444     scalarField maxEdgeLenSqr(nLevels, -GREAT2);
 445
 446     forAll(cellLevel_, cellI)
 447     {
 448         const label cLevel = cellLevel_[cellI];
 449
 450         const labelList& cEdges = mesh_.cellEdges()[cellI];
 451
 452         forAll(cEdges, i)
 453         {
 454             const edge& e = mesh_.edges()[cEdges[i]];
 455
 456             scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
 457
 458             maxEdgeLenSqr[cLevel] = max(maxEdgeLenSqr[cLevel], edgeLenSqr);
 459         }
 460     }
 461
 462     Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>());
 463     Pstream::listCombineScatter(maxEdgeLenSqr);
 464
 465     if (debug)
 466     {
 467         Pout<< "hexRef8::getLevel0EdgeLength() :"
 468             << " Crappy Edgelengths (squared) per refinementlevel:"
 469             << maxEdgeLenSqr << endl;
 470     }
 471
 472
 473     // 3. Combine the two sets of lengths
 474     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 475
 476     forAll(typEdgeLenSqr, levelI)
 477     {
 478         if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
 479         {
 480             typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
 481         }
 482     }
 483
 484     if (debug)
 485     {
 486         Pout<< "hexRef8::getLevel0EdgeLength() :"
 487             << " Final Edgelengths (squared) per refinementlevel:"
 488             << typEdgeLenSqr << endl;
 489     }
 490
 491     // Find lowest level present
 492     scalar level0Size = -1;
 493
 494     forAll(typEdgeLenSqr, levelI)
 495     {
 496         scalar lenSqr = typEdgeLenSqr[levelI];
 497
 498         if (lenSqr < GREAT2)
 499         {
 500             level0Size = Foam::sqrt(lenSqr)*(1<<levelI);
 501
 502             if (debug)
 503             {
 504                 Pout<< "hexRef8::getLevel0EdgeLength() :"
 505                     << " For level:" << levelI
 506                     << " have edgeLen:" << Foam::sqrt(lenSqr)
 507                     << " with equivalent level0 len:" << level0Size
 508                     << endl;
 509             }
 510             break;
 511         }
 512     }
 513
 514     if (level0Size == -1)
 515     {
 516         FatalErrorIn("hexRef8::getLevel0EdgeLength()")
 517             << "Problem : typEdgeLenSqr:" << typEdgeLenSqr << abort(FatalError);
 518     }
 519
 520     return level0Size;
 521 }
 522
 523
 524 // Check whether pointI is an anchor on cellI.
 525 // If it is not check whether any other point on the face is an anchor cell.
 526 Foam::label Foam::hexRef8::getAnchorCell
 527 (
 528     const labelListList& cellAnchorPoints,
 529     const labelListList& cellAddedCells,
 530     const label cellI,
 531     const label faceI,
 532     const label pointI
 533 ) const
 534 {
 535     if (cellAnchorPoints[cellI].size() > 0)
 536     {
 537         label index = findIndex(cellAnchorPoints[cellI], pointI);
 538
 539         if (index != -1)
 540         {
 541             return cellAddedCells[cellI][index];
 542         }
 543
 544
 545         // pointI is not an anchor cell.
 546         // Maybe we are already a refined face so check all the face
 547         // vertices.
 548         const face& f = mesh_.faces()[faceI];
 549
 550         forAll(f, fp)
 551         {
 552             label index = findIndex(cellAnchorPoints[cellI], f[fp]);
 553
 554             if (index != -1)
 555             {
 556                 return cellAddedCells[cellI][index];
 557             }
 558         }
 559
 560         // Problem.
 561
 562         // Pick up points of the cell
 563         const labelList cPoints(cellPoints(cellI));
 564
 565         Perr<< "cell:" << cellI << " points:" << endl;
 566         forAll(cPoints, i)
 567         {
 568             label pointI = cPoints[i];
 569
 570             Perr<< "    " << pointI << " coord:" << mesh_.points()[pointI]
 571                 << nl;
 572         }
 573         Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI]
 574             << endl;
 575
 576         FatalErrorIn("hexRef8::getAnchorCell(..)")
 577             << "Could not find point " << pointI
 578             << " in the anchorPoints for cell " << cellI << endl
 579             << "Does your original mesh obey the 2:1 constraint and"
 580             << " did you use consistentRefinement to make your cells to refine"
 581             << " obey this constraint as well?"
 582             << abort(FatalError);
 583
 584         return -1;
 585     }
 586     else
 587     {
 588         return cellI;
 589     }
 590 }
 591
 592
 593 // Get new owner and neighbour
 594 void Foam::hexRef8::getFaceNeighbours
 595 (
 596     const labelListList& cellAnchorPoints,
 597     const labelListList& cellAddedCells,
 598     const label faceI,
 599     const label pointI,
 600
 601     label& own,
 602     label& nei
 603 ) const
 604 {
 605     // Is owner split?
 606     own = getAnchorCell
 607     (
 608         cellAnchorPoints,
 609         cellAddedCells,
 610         mesh_.faceOwner()[faceI],
 611         faceI,
 612         pointI
 613     );
 614
 615     if (mesh_.isInternalFace(faceI))
 616     {
 617         nei = getAnchorCell
 618         (
 619             cellAnchorPoints,
 620             cellAddedCells,
 621             mesh_.faceNeighbour()[faceI],
 622             faceI,
 623             pointI
 624         );
 625     }
 626     else
 627     {
 628         nei = -1;
 629     }
 630 }
 631
 632
 633 // Get point with the lowest pointLevel
 634 Foam::label Foam::hexRef8::findMinLevel(const labelList& f) const
 635 {
 636     label minLevel = labelMax;
 637     label minFp = -1;
 638
 639     forAll(f, fp)
 640     {
 641         label level = pointLevel_[f[fp]];
 642
 643         if (level < minLevel)
 644         {
 645             minLevel = level;
 646             minFp = fp;
 647         }
 648     }
 649
 650     return minFp;
 651 }
 652
 653
 654 // Get point with the highest pointLevel
 655 Foam::label Foam::hexRef8::findMaxLevel(const labelList& f) const
 656 {
 657     label maxLevel = labelMin;
 658     label maxFp = -1;
 659
 660     forAll(f, fp)
 661     {
 662         label level = pointLevel_[f[fp]];
 663
 664         if (level > maxLevel)
 665         {
 666             maxLevel = level;
 667             maxFp = fp;
 668         }
 669     }
 670
 671     return maxFp;
 672 }
 673
 674
 675 Foam::label Foam::hexRef8::countAnchors
 676 (
 677     const labelList& f,
 678     const label anchorLevel
 679 ) const
 680 {
 681     label nAnchors = 0;
 682
 683     forAll(f, fp)
 684     {
 685         if (pointLevel_[f[fp]] <= anchorLevel)
 686         {
 687             nAnchors++;
 688         }
 689     }
 690     return nAnchors;
 691 }
 692
 693
 694 // Find point with certain pointLevel. Skip any higher levels.
 695 Foam::label Foam::hexRef8::findLevel
 696 (
 697     const face& f,
 698     const label startFp,
 699     const bool searchForward,
 700     const label wantedLevel
 701 ) const
 702 {
 703     label fp = startFp;
 704
 705     forAll(f, i)
 706     {
 707         label pointI = f[fp];
 708
 709         if (pointLevel_[pointI] < wantedLevel)
 710         {
 711             FatalErrorIn("hexRef8::findLevel")
 712                 << "face:" << f
 713                 << " level:" << IndirectList<label>(pointLevel_, f)()
 714                 << " startFp:" << startFp
 715                 << " wantedLevel:" << wantedLevel
 716                 << abort(FatalError);
 717         }
 718         else if (pointLevel_[pointI] == wantedLevel)
 719         {
 720             return fp;
 721         }
 722
 723         if (searchForward)
 724         {
 725             fp = f.fcIndex(fp);
 726         }
 727         else
 728         {
 729             fp = f.rcIndex(fp);
 730         }
 731     }
 732
 733     FatalErrorIn("hexRef8::findLevel")
 734         << "face:" << f
 735         << " level:" << IndirectList<label>(pointLevel_, f)()
 736         << " startFp:" << startFp
 737         << " wantedLevel:" << wantedLevel
 738         << abort(FatalError);
 739
 740     return -1;
 741 }
 742
 743
 744 // Gets cell level such that the face has four points <= level.
 745 Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
 746 {
 747     const face& f = mesh_.faces()[faceI];
 748
 749     if (f.size() <= 4)
 750     {
 751         return pointLevel_[f[findMaxLevel(f)]];
 752     }
 753     else
 754     {
 755         label ownLevel = cellLevel_[mesh_.faceOwner()[faceI]];
 756
 757         if (countAnchors(f, ownLevel) == 4)
 758         {
 759             return ownLevel;
 760         }
 761         else if (countAnchors(f, ownLevel+1) == 4)
 762         {
 763             return ownLevel+1;
 764         }
 765         else
 766         {
 767             //FatalErrorIn("hexRef8::getAnchorLevel(const label) const")
 768             //    << "face:" << faceI
 769             //    << " centre:" << mesh_.faceCentres()[faceI]
 770             //    << " verts:" << f
 771             //    << " point levels:" << IndirectList<label>(pointLevel_, f)()
 772             //    << " own:" << mesh_.faceOwner()[faceI]
 773             //    << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]]
 774             //    << abort(FatalError);
 775
 776             return -1;
 777         }
 778     }
 779 }
 780
 781
 782 void Foam::hexRef8::checkInternalOrientation
 783 (
 784     polyTopoChange& meshMod,
 785     const label cellI,
 786     const label faceI,
 787     const point& ownPt,
 788     const point& neiPt,
 789     const face& newFace
 790 )
 791 {
 792     face compactFace(identity(newFace.size()));
 793     pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)());
 794
 795     vector n(compactFace.normal(compactPoints));
 796
 797     vector dir(neiPt - ownPt);
 798
 799     if ((dir & n) < 0)
 800     {
 801         FatalErrorIn("checkInternalOrientation(..)")
 802             << "cell:" << cellI << " old face:" << faceI
 803             << " newFace:" << newFace << endl
 804             << " coords:" << compactPoints
 805             << " ownPt:" << ownPt
 806             << " neiPt:" << neiPt
 807             << abort(FatalError);
 808     }
 809
 810     vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
 811
 812     scalar s = (fcToOwn&n) / (dir&n);
 813
 814     if (s < 0.1 || s > 0.9)
 815     {
 816         FatalErrorIn("checkInternalOrientation(..)")
 817             << "cell:" << cellI << " old face:" << faceI
 818             << " newFace:" << newFace << endl
 819             << " coords:" << compactPoints
 820             << " ownPt:" << ownPt
 821             << " neiPt:" << neiPt
 822             << " s:" << s
 823             << abort(FatalError);
 824     }
 825 }
 826
 827
 828 void Foam::hexRef8::checkBoundaryOrientation
 829 (
 830     polyTopoChange& meshMod,
 831     const label cellI,
 832     const label faceI,
 833     const point& ownPt,
 834     const point& boundaryPt,
 835     const face& newFace
 836 )
 837 {
 838     face compactFace(identity(newFace.size()));
 839     pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)());
 840
 841     vector n(compactFace.normal(compactPoints));
 842
 843     vector dir(boundaryPt - ownPt);
 844
 845     if ((dir & n) < 0)
 846     {
 847         FatalErrorIn("checkBoundaryOrientation(..)")
 848             << "cell:" << cellI << " old face:" << faceI
 849             << " newFace:" << newFace
 850             << " coords:" << compactPoints
 851             << " ownPt:" << ownPt
 852             << " boundaryPt:" << boundaryPt
 853             << abort(FatalError);
 854     }
 855
 856     vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
 857
 858     scalar s = (fcToOwn&dir) / magSqr(dir);
 859
 860     if (s < 0.7 || s > 1.3)
 861     {
 862         WarningIn("checkBoundaryOrientation(..)")
 863             << "cell:" << cellI << " old face:" << faceI
 864             << " newFace:" << newFace
 865             << " coords:" << compactPoints
 866             << " ownPt:" << ownPt
 867             << " boundaryPt:" << boundaryPt
 868             << " s:" << s
 869             << endl;
 870             //<< abort(FatalError);
 871     }
 872 }
 873
 874
 875 // If p0 and p1 are existing vertices check if edge is split and insert
 876 // splitPoint.
 877 void Foam::hexRef8::insertEdgeSplit
 878 (
 879     const labelList& edgeMidPoint,
 880     const label p0,
 881     const label p1,
 882     DynamicList<label>& verts
 883 ) const
 884 {
 885     if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
 886     {
 887         label edgeI = meshTools::findEdge(mesh_, p0, p1);
 888
 889         if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
 890         {
 891             verts.append(edgeMidPoint[edgeI]);
 892         }
 893     }
 894 }
 895
 896
 897 // Internal faces are one per edge between anchor points. So one per midPoint
 898 // between the anchor points. Here we store the information on the midPoint
 899 // and if we have enough information:
 900 // - two anchors
 901 // - two face mid points
 902 // we add the face. Note that this routine can get called anywhere from
 903 // two times (two unrefined faces) to four times (two refined faces) so
 904 // the first call that adds the information creates the face.
 905 Foam::label Foam::hexRef8::storeMidPointInfo
 906 (
 907     const labelListList& cellAnchorPoints,
 908     const labelListList& cellAddedCells,
 909     const labelList& cellMidPoint,
 910     const labelList& edgeMidPoint,
 911     const label cellI,
 912     const label faceI,
 913     const bool faceOrder,
 914     const label edgeMidPointI,
 915     const label anchorPointI,
 916     const label faceMidPointI,
 917
 918     Map<edge>& midPointToAnchors,
 919     Map<edge>& midPointToFaceMids,
 920     polyTopoChange& meshMod
 921 ) const
 922 {
 923     // See if need to store anchors.
 924
 925     bool changed = false;
 926     bool haveTwoAnchors = false;
 927
 928     Map<edge>::iterator edgeMidFnd = midPointToAnchors.find(edgeMidPointI);
 929
 930     if (edgeMidFnd == midPointToAnchors.end())
 931     {
 932         midPointToAnchors.insert(edgeMidPointI, edge(anchorPointI, -1));
 933     }
 934     else
 935     {
 936         edge& e = edgeMidFnd();
 937
 938         if (anchorPointI != e[0])
 939         {
 940             if (e[1] == -1)
 941             {
 942                 e[1] = anchorPointI;
 943                 changed = true;
 944             }
 945         }
 946
 947         if (e[0] != -1 && e[1] != -1)
 948         {
 949             haveTwoAnchors = true;
 950         }
 951     }
 952
 953     bool haveTwoFaceMids = false;
 954
 955     Map<edge>::iterator faceMidFnd = midPointToFaceMids.find(edgeMidPointI);
 956
 957     if (faceMidFnd == midPointToFaceMids.end())
 958     {
 959         midPointToFaceMids.insert(edgeMidPointI, edge(faceMidPointI, -1));
 960     }
 961     else
 962     {
 963         edge& e = faceMidFnd();
 964
 965         if (faceMidPointI != e[0])
 966         {
 967             if (e[1] == -1)
 968             {
 969                 e[1] = faceMidPointI;
 970                 changed = true;
 971             }
 972         }
 973
 974         if (e[0] != -1 && e[1] != -1)
 975         {
 976             haveTwoFaceMids = true;
 977         }
 978     }
 979
 980     // Check if this call of storeMidPointInfo is the one that completed all
 981     // the nessecary information.
 982
 983     if (changed && haveTwoAnchors && haveTwoFaceMids)
 984     {
 985         const edge& anchors = midPointToAnchors[edgeMidPointI];
 986         const edge& faceMids = midPointToFaceMids[edgeMidPointI];
 987
 988         label otherFaceMidPointI = faceMids.otherVertex(faceMidPointI);
 989
 990         // Create face consistent with anchorI being the owner.
 991         // Note that the edges between the edge mid point and the face mids
 992         // might be marked for splitting. Note that these edge splits cannot
 993         // be between cellMid and face mids.
 994
 995         DynamicList<label> newFaceVerts(4);
 996         if (faceOrder == (mesh_.faceOwner()[faceI] == cellI))
 997         {
 998             newFaceVerts.append(faceMidPointI);
 999
1000             // Check & insert edge split if any
1001             insertEdgeSplit
1002             (
1003                 edgeMidPoint,
1004                 faceMidPointI,  // edge between faceMid and
1005                 edgeMidPointI,  // edgeMid
1006                 newFaceVerts
1007             );
1008
1009             newFaceVerts.append(edgeMidPointI);
1010
1011             insertEdgeSplit
1012             (
1013                 edgeMidPoint,
1014                 edgeMidPointI,
1015                 otherFaceMidPointI,
1016                 newFaceVerts
1017             );
1018
1019             newFaceVerts.append(otherFaceMidPointI);
1020             newFaceVerts.append(cellMidPoint[cellI]);
1021         }
1022         else
1023         {
1024             newFaceVerts.append(otherFaceMidPointI);
1025
1026             insertEdgeSplit
1027             (
1028                 edgeMidPoint,
1029                 otherFaceMidPointI,
1030                 edgeMidPointI,
1031                 newFaceVerts
1032             );
1033
1034             newFaceVerts.append(edgeMidPointI);
1035
1036             insertEdgeSplit
1037             (
1038                 edgeMidPoint,
1039                 edgeMidPointI,
1040                 faceMidPointI,
1041                 newFaceVerts
1042             );
1043
1044             newFaceVerts.append(faceMidPointI);
1045             newFaceVerts.append(cellMidPoint[cellI]);
1046         }
1047
1048         face newFace;
1049         newFace.transfer(newFaceVerts.shrink());
1050         newFaceVerts.clear();
1051
1052         label anchorCell0 = getAnchorCell
1053         (
1054             cellAnchorPoints,
1055             cellAddedCells,
1056             cellI,
1057             faceI,
1058             anchorPointI
1059         );
1060         label anchorCell1 = getAnchorCell
1061         (
1062             cellAnchorPoints,
1063             cellAddedCells,
1064             cellI,
1065             faceI,
1066             anchors.otherVertex(anchorPointI)
1067         );
1068
1069
1070         label own, nei;
1071         point ownPt, neiPt;
1072
1073         if (anchorCell0 < anchorCell1)
1074         {
1075             own = anchorCell0;
1076             nei = anchorCell1;
1077
1078             ownPt = mesh_.points()[anchorPointI];
1079             neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1080
1081         }
1082         else
1083         {
1084             own = anchorCell1;
1085             nei = anchorCell0;
1086             newFace = newFace.reverseFace();
1087
1088             ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1089             neiPt = mesh_.points()[anchorPointI];
1090         }
1091
1092         if (debug)
1093         {
1094             point ownPt, neiPt;
1095
1096             if (anchorCell0 < anchorCell1)
1097             {
1098                 ownPt = mesh_.points()[anchorPointI];
1099                 neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1100             }
1101             else
1102             {
1103                 ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1104                 neiPt = mesh_.points()[anchorPointI];
1105             }
1106
1107             checkInternalOrientation
1108             (
1109                 meshMod,
1110                 cellI,
1111                 faceI,
1112                 ownPt,
1113                 neiPt,
1114                 newFace
1115             );
1116         }
1117
1118         return addInternalFace
1119         (
1120             meshMod,
1121             faceI,
1122             anchorPointI,
1123             newFace,
1124             own,
1125             nei
1126         );
1127     }
1128     else
1129     {
1130         return -1;
1131     }
1132 }
1133
1134
1135 // Creates all the 12 internal faces for cellI.
1136 void Foam::hexRef8::createInternalFaces
1137 (
1138     const labelListList& cellAnchorPoints,
1139     const labelListList& cellAddedCells,
1140     const labelList& cellMidPoint,
1141     const labelList& faceMidPoint,
1142     const labelList& faceAnchorLevel,
1143     const labelList& edgeMidPoint,
1144     const label cellI,
1145
1146     polyTopoChange& meshMod
1147 ) const
1148 {
1149     // Find in every face the cellLevel+1 points (from edge subdivision)
1150     // and the anchor points.
1151
1152     const cell& cFaces = mesh_.cells()[cellI];
1153     const label cLevel = cellLevel_[cellI];
1154
1155     // From edge mid to anchor points
1156     Map<edge> midPointToAnchors(24);
1157     // From edge mid to face mids
1158     Map<edge> midPointToFaceMids(24);
1159
1160     // Running count of number of internal faces added so far.
1161     label nFacesAdded = 0;
1162
1163     forAll(cFaces, i)
1164     {
1165         label faceI = cFaces[i];
1166
1167         const face& f = mesh_.faces()[faceI];
1168         const labelList& fEdges = mesh_.faceEdges()[faceI];
1169
1170         // We are on the cellI side of face f. The face will have 1 or 4
1171         // cLevel points and lots of higher numbered ones.
1172
1173         label faceMidPointI = -1;
1174
1175         label nAnchors = countAnchors(f, cLevel);
1176
1177         if (nAnchors == 1)
1178         {
1179             // Only one anchor point. So the other side of the face has already
1180             // been split using cLevel+1 and cLevel+2 points.
1181
1182             // Find the one anchor.
1183             label anchorFp = -1;
1184
1185             forAll(f, fp)
1186             {
1187                 if (pointLevel_[f[fp]] <= cLevel)
1188                 {
1189                     anchorFp = fp;
1190                     break;
1191                 }
1192             }
1193
1194             // Now the face mid point is the second cLevel+1 point
1195             label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1);
1196             label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1);
1197
1198             faceMidPointI = f[faceMid];
1199         }
1200         else if (nAnchors == 4)
1201         {
1202             // There is no face middle yet but the face will be marked for
1203             // splitting.
1204
1205             faceMidPointI = faceMidPoint[faceI];
1206         }
1207         else
1208         {
1209             FatalErrorIn("createInternalFaces")
1210                 << "nAnchors:" << nAnchors
1211                 << " faceI:" << faceI
1212                 << abort(FatalError);
1213         }
1214
1215
1216
1217         // Now loop over all the anchors (might be just one) and store
1218         // the edge mids connected to it. storeMidPointInfo will collect
1219         // all the info and combine it all.
1220
1221         forAll(f, fp0)
1222         {
1223             label point0 = f[fp0];
1224
1225             if (pointLevel_[point0] <= cLevel)
1226             {
1227                 // Anchor.
1228
1229                 // Walk forward
1230                 // ~~~~~~~~~~~~
1231                 // to cLevel+1 or edgeMidPoint of this level.
1232
1233
1234                 label edgeMidPointI = -1;
1235
1236                 label fp1 = f.fcIndex(fp0);
1237
1238                 if (pointLevel_[f[fp1]] <= cLevel)
1239                 {
1240                     // Anchor. Edge will be split.
1241                     label edgeI = fEdges[fp0];
1242
1243                     edgeMidPointI = edgeMidPoint[edgeI];
1244
1245                     if (edgeMidPointI == -1)
1246                     {
1247                         const labelList cPoints(cellPoints(cellI));
1248
1249                         FatalErrorIn("createInternalFaces(..)")
1250                             << "cell:" << cellI << " cLevel:" << cLevel
1251                             << " cell points:" << cPoints
1252                             << " pointLevel:"
1253                             << IndirectList<label>(pointLevel_, cPoints)()
1254                             << " face:" << faceI
1255                             << " f:" << f
1256                             << " pointLevel:"
1257                             << IndirectList<label>(pointLevel_, f)()
1258                             << " faceAnchorLevel:" << faceAnchorLevel[faceI]
1259                             << " faceMidPoint:" << faceMidPoint[faceI]
1260                             << " faceMidPointI:" << faceMidPointI
1261                             << " fp:" << fp0
1262                             << abort(FatalError);
1263                     }
1264                 }
1265                 else
1266                 {
1267                     // Search forward in face to clevel+1
1268                     label edgeMid = findLevel(f, fp1, true, cLevel+1);
1269
1270                     edgeMidPointI = f[edgeMid];
1271                 }
1272
1273                 label newFaceI = storeMidPointInfo
1274                 (
1275                     cellAnchorPoints,
1276                     cellAddedCells,
1277                     cellMidPoint,
1278                     edgeMidPoint,
1279
1280                     cellI,
1281                     faceI,
1282                     true,                   // mid point after anchor
1283                     edgeMidPointI,          // edgemid
1284                     point0,                 // anchor
1285                     faceMidPointI,
1286
1287                     midPointToAnchors,
1288                     midPointToFaceMids,
1289                     meshMod
1290                 );
1291
1292                 if (newFaceI != -1)
1293                 {
1294                     nFacesAdded++;
1295
1296                     if (nFacesAdded == 12)
1297                     {
1298                         break;
1299                     }
1300                 }
1301
1302
1303
1304                 // Walk backward
1305                 // ~~~~~~~~~~~~~
1306
1307                 label fpMin1 = f.rcIndex(fp0);
1308
1309                 if (pointLevel_[f[fpMin1]] <= cLevel)
1310                 {
1311                     // Anchor. Edge will be split.
1312                     label edgeI = fEdges[fpMin1];
1313
1314                     edgeMidPointI = edgeMidPoint[edgeI];
1315
1316                     if (edgeMidPointI == -1)
1317                     {
1318                         const labelList cPoints(cellPoints(cellI));
1319
1320                         FatalErrorIn("createInternalFaces(..)")
1321                             << "cell:" << cellI << " cLevel:" << cLevel
1322                             << " cell points:" << cPoints
1323                             << " pointLevel:"
1324                             << IndirectList<label>(pointLevel_, cPoints)()
1325                             << " face:" << faceI
1326                             << " f:" << f
1327                             << " pointLevel:"
1328                             << IndirectList<label>(pointLevel_, f)()
1329                             << " faceAnchorLevel:" << faceAnchorLevel[faceI]
1330                             << " faceMidPoint:" << faceMidPoint[faceI]
1331                             << " faceMidPointI:" << faceMidPointI
1332                             << " fp:" << fp0
1333                             << abort(FatalError);
1334                     }
1335                 }
1336                 else
1337                 {
1338                     // Search back to clevel+1
1339                     label edgeMid = findLevel(f, fpMin1, false, cLevel+1);
1340
1341                     edgeMidPointI = f[edgeMid];
1342                 }
1343
1344                 newFaceI = storeMidPointInfo
1345                 (
1346                     cellAnchorPoints,
1347                     cellAddedCells,
1348                     cellMidPoint,
1349                     edgeMidPoint,
1350
1351                     cellI,
1352                     faceI,
1353                     false,                  // mid point before anchor
1354                     edgeMidPointI,          // edgemid
1355                     point0,                 // anchor
1356                     faceMidPointI,
1357
1358                     midPointToAnchors,
1359                     midPointToFaceMids,
1360                     meshMod
1361                 );
1362
1363                 if (newFaceI != -1)
1364                 {
1365                     nFacesAdded++;
1366
1367                     if (nFacesAdded == 12)
1368                     {
1369                         break;
1370                     }
1371                 }
1372             }   // done anchor
1373         }   // done face
1374
1375         if (nFacesAdded == 12)
1376         {
1377             break;
1378         }
1379     }
1380 }
1381
1382
1383 void Foam::hexRef8::walkFaceToMid
1384 (
1385     const labelList& edgeMidPoint,
1386     const label cLevel,
1387     const label faceI,
1388     const label startFp,
1389     DynamicList<label>& faceVerts
1390 ) const
1391 {
1392     const face& f = mesh_.faces()[faceI];
1393     const labelList& fEdges = mesh_.faceEdges()[faceI];
1394
1395     label fp = startFp;
1396
1397     // Starting from fp store all (1 or 2) vertices until where the face
1398     // gets split
1399
1400     while (true)
1401     {
1402         if (edgeMidPoint[fEdges[fp]] >= 0)
1403         {
1404             faceVerts.append(edgeMidPoint[fEdges[fp]]);
1405         }
1406
1407         fp = f.fcIndex(fp);
1408
1409         if (pointLevel_[f[fp]] <= cLevel)
1410         {
1411             // Next anchor. Have already append split point on edge in code
1412             // above.
1413             return;
1414         }
1415         else if (pointLevel_[f[fp]] == cLevel+1)
1416         {
1417             // Mid level
1418             faceVerts.append(f[fp]);
1419
1420             return;
1421         }
1422         else if (pointLevel_[f[fp]] == cLevel+2)
1423         {
1424             // Store and continue to cLevel+1.
1425             faceVerts.append(f[fp]);
1426         }
1427     }
1428 }
1429
1430
1431 // Same as walkFaceToMid but now walk back.
1432 void Foam::hexRef8::walkFaceFromMid
1433 (
1434     const labelList& edgeMidPoint,
1435     const label cLevel,
1436     const label faceI,
1437     const label startFp,
1438     DynamicList<label>& faceVerts
1439 ) const
1440 {
1441     const face& f = mesh_.faces()[faceI];
1442     const labelList& fEdges = mesh_.faceEdges()[faceI];
1443
1444     label fp = f.rcIndex(startFp);
1445
1446     while (true)
1447     {
1448         if (pointLevel_[f[fp]] <= cLevel)
1449         {
1450             // anchor.
1451             break;
1452         }
1453         else if (pointLevel_[f[fp]] == cLevel+1)
1454         {
1455             // Mid level
1456             faceVerts.append(f[fp]);
1457             break;
1458         }
1459         else if (pointLevel_[f[fp]] == cLevel+2)
1460         {
1461             // Continue to cLevel+1.
1462         }
1463         fp = f.rcIndex(fp);
1464     }
1465
1466     // Store
1467     while (true)
1468     {
1469         if (edgeMidPoint[fEdges[fp]] >= 0)
1470         {
1471             faceVerts.append(edgeMidPoint[fEdges[fp]]);
1472         }
1473
1474         fp = f.fcIndex(fp);
1475
1476         if (fp == startFp)
1477         {
1478             break;
1479         }
1480         faceVerts.append(f[fp]);
1481     }
1482 }
1483
1484
1485 // Updates refineCell (cells marked for refinement) so across all faces
1486 // there will be 2:1 consistency after refinement.
1487 Foam::label Foam::hexRef8::faceConsistentRefinement
1488 (
1489     const bool maxSet,
1490     PackedList<1>& refineCell
1491 ) const
1492 {
1493     label nChanged = 0;
1494
1495     // Internal faces.
1496     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1497     {
1498         label own = mesh_.faceOwner()[faceI];
1499         label ownLevel = cellLevel_[own] + refineCell.get(own);
1500
1501         label nei = mesh_.faceNeighbour()[faceI];
1502         label neiLevel = cellLevel_[nei] + refineCell.get(nei);
1503
1504         if (ownLevel > (neiLevel+1))
1505         {
1506             if (maxSet)
1507             {
1508                 refineCell.set(nei, 1);
1509             }
1510             else
1511             {
1512                 refineCell.set(own, 0);
1513             }
1514             nChanged++;
1515         }
1516         else if (neiLevel > (ownLevel+1))
1517         {
1518             if (maxSet)
1519             {
1520                 refineCell.set(own, 1);
1521             }
1522             else
1523             {
1524                 refineCell.set(nei, 0);
1525             }
1526             nChanged++;
1527         }
1528     }
1529
1530
1531     // Coupled faces. Swap owner level to get neighbouring cell level.
1532     // (only boundary faces of neiLevel used)
1533     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1534
1535     forAll(neiLevel, i)
1536     {
1537         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1538
1539         neiLevel[i] = cellLevel_[own] + refineCell.get(own);
1540     }
1541
1542     // Swap to neighbour
1543     syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
1544
1545     // Now we have neighbour value see which cells need refinement
1546     forAll(neiLevel, i)
1547     {
1548         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1549         label ownLevel = cellLevel_[own] + refineCell.get(own);
1550
1551         if (ownLevel > (neiLevel[i]+1))
1552         {
1553             if (!maxSet)
1554             {
1555                 refineCell.set(own, 0);
1556                 nChanged++;
1557             }
1558         }
1559         else if (neiLevel[i] > (ownLevel+1))
1560         {
1561             if (maxSet)
1562             {
1563                 refineCell.set(own, 1);
1564                 nChanged++;
1565             }
1566         }
1567     }
1568
1569     return nChanged;
1570 }
1571
1572
1573 // Debug: check if wanted refinement is compatible with 2:1
1574 void Foam::hexRef8::checkWantedRefinementLevels
1575 (
1576     const labelList& cellsToRefine
1577 ) const
1578 {
1579     PackedList<1> refineCell(mesh_.nCells(), 0);
1580     forAll(cellsToRefine, i)
1581     {
1582         refineCell.set(cellsToRefine[i], 1);
1583     }
1584
1585     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1586     {
1587         label own = mesh_.faceOwner()[faceI];
1588         label ownLevel = cellLevel_[own] + refineCell.get(own);
1589
1590         label nei = mesh_.faceNeighbour()[faceI];
1591         label neiLevel = cellLevel_[nei] + refineCell.get(nei);
1592
1593         if (mag(ownLevel-neiLevel) > 1)
1594         {
1595             FatalErrorIn
1596             (
1597                 "hexRef8::checkWantedRefinementLevels(const labelList&)"
1598             )   << "cell:" << own
1599                 << " current level:" << cellLevel_[own]
1600                 << " level after refinement:" << ownLevel
1601                 << nl
1602                 << "neighbour cell:" << nei
1603                 << " current level:" << cellLevel_[nei]
1604                 << " level after refinement:" << neiLevel
1605                 << nl
1606                 << "which does not satisfy 2:1 constraints anymore."
1607                 << abort(FatalError);
1608         }
1609     }
1610
1611     // Coupled faces. Swap owner level to get neighbouring cell level.
1612     // (only boundary faces of neiLevel used)
1613     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1614
1615     forAll(neiLevel, i)
1616     {
1617         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1618
1619         neiLevel[i] = cellLevel_[own] + refineCell.get(own);
1620     }
1621
1622     // Swap to neighbour
1623     syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
1624
1625     // Now we have neighbour value see which cells need refinement
1626     forAll(neiLevel, i)
1627     {
1628         label faceI = i + mesh_.nInternalFaces();
1629
1630         label own = mesh_.faceOwner()[faceI];
1631         label ownLevel = cellLevel_[own] + refineCell.get(own);
1632
1633         if (mag(ownLevel - neiLevel[i]) > 1)
1634         {
1635             label patchI = mesh_.boundaryMesh().whichPatch(faceI);
1636
1637             FatalErrorIn
1638             (
1639                 "hexRef8::checkWantedRefinementLevels(const labelList&)"
1640             )   << "Celllevel does not satisfy 2:1 constraint."
1641                 << " On coupled face "
1642                 << faceI
1643                 << " on patch " << patchI << " "
1644                 << mesh_.boundaryMesh()[patchI].name()
1645                 << " owner cell " << own
1646                 << " current level:" << cellLevel_[own]
1647                 << " level after refinement:" << ownLevel
1648                 << nl
1649                 << " (coupled) neighbour cell will get refinement "
1650                 << neiLevel[i]
1651                 << abort(FatalError);
1652         }
1653     }
1654 }
1655
1656
1657 // Set instance for mesh files
1658 void Foam::hexRef8::setInstance(const fileName& inst)
1659 {
1660     if (debug)
1661     {
1662         Pout<< "hexRef8::setInstance(const fileName& inst) : "
1663             << "Resetting file instance to " << inst << endl;
1664     }
1665
1666     cellLevel_.instance() = inst;
1667     pointLevel_.instance() = inst;
1668     history_.instance() = inst;
1669 }
1670
1671
1672 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1673
1674 // Construct from mesh, read refinement data
1675 Foam::hexRef8::hexRef8(const polyMesh& mesh)
1676 :
1677     mesh_(mesh),
1678     cellLevel_
1679     (
1680         IOobject
1681         (
1682             "cellLevel",
1683             mesh_.facesInstance(),
1684             polyMesh::meshSubDir,
1685             mesh_,
1686             IOobject::READ_IF_PRESENT,
1687             IOobject::AUTO_WRITE
1688         ),
1689         labelList(mesh_.nCells(), 0)
1690     ),
1691     pointLevel_
1692     (
1693         IOobject
1694         (
1695             "pointLevel",
1696             mesh_.facesInstance(),
1697             polyMesh::meshSubDir,
1698             mesh_,
1699             IOobject::READ_IF_PRESENT,
1700             IOobject::AUTO_WRITE
1701         ),
1702         labelList(mesh_.nPoints(), 0)
1703     ),
1704     level0Edge_(getLevel0EdgeLength()),
1705     history_
1706     (
1707         IOobject
1708         (
1709             "refinementHistory",
1710             mesh_.facesInstance(),
1711             polyMesh::meshSubDir,
1712             mesh_,
1713             IOobject::READ_IF_PRESENT,
1714             IOobject::AUTO_WRITE
1715         ),
1716         mesh_.nCells()    // All cells visible if could not be read
1717     ),
1718     faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
1719     savedPointLevel_(0),
1720     savedCellLevel_(0)
1721 {
1722     if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1723     {
1724         FatalErrorIn
1725         (
1726             "hexRef8::hexRef8(const polyMesh&)"
1727         )   << "History enabled but number of visible cells in it "
1728             << history_.visibleCells().size()
1729             << " is not equal to the number of cells in the mesh "
1730             << mesh_.nCells()
1731             << abort(FatalError);
1732     }
1733
1734     if
1735     (
1736         cellLevel_.size() != mesh_.nCells()
1737      || pointLevel_.size() != mesh_.nPoints()
1738     )
1739     {
1740         FatalErrorIn
1741         (
1742             "hexRef8::hexRef8(const polyMesh&)"
1743         )   << "Restarted from inconsistent cellLevel or pointLevel files."
1744             << endl
1745             << "Number of cells in mesh:" << mesh_.nCells()
1746             << " does not equal size of cellLevel:" << cellLevel_.size() << endl
1747             << "Number of points in mesh:" << mesh_.nPoints()
1748             << " does not equal size of pointLevel:" << pointLevel_.size()
1749             << abort(FatalError);
1750     }
1751
1752
1753     // Check refinement levels for consistency
1754     checkRefinementLevels(-1, labelList(0));
1755
1756
1757     // Check initial mesh for consistency
1758
1759     //if (debug)
1760     {
1761         checkMesh();
1762     }
1763 }
1764
1765
1766 // Construct from components
1767 Foam::hexRef8::hexRef8
1768 (
1769     const polyMesh& mesh,
1770     const labelList& cellLevel,
1771     const labelList& pointLevel,
1772     const refinementHistory& history
1773 )
1774 :
1775     mesh_(mesh),
1776     cellLevel_
1777     (
1778         IOobject
1779         (
1780             "cellLevel",
1781             mesh_.facesInstance(),
1782             polyMesh::meshSubDir,
1783             mesh_,
1784             IOobject::NO_READ,
1785             IOobject::AUTO_WRITE
1786         ),
1787         cellLevel
1788     ),
1789     pointLevel_
1790     (
1791         IOobject
1792         (
1793             "pointLevel",
1794             mesh_.facesInstance(),
1795             polyMesh::meshSubDir,
1796             mesh_,
1797             IOobject::NO_READ,
1798             IOobject::AUTO_WRITE
1799         ),
1800         pointLevel
1801     ),
1802     level0Edge_(getLevel0EdgeLength()),
1803     history_
1804     (
1805         IOobject
1806         (
1807             "refinementHistory",
1808             mesh_.facesInstance(),
1809             polyMesh::meshSubDir,
1810             mesh_,
1811             IOobject::NO_READ,
1812             IOobject::AUTO_WRITE
1813         ),
1814         history
1815     ),
1816     faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
1817     savedPointLevel_(0),
1818     savedCellLevel_(0)
1819 {
1820     if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1821     {
1822         FatalErrorIn
1823         (
1824             "hexRef8::hexRef8(const polyMesh&, const labelList&"
1825             ", const labelList&, const refinementHistory&)"
1826         )   << "History enabled but number of visible cells in it "
1827             << history_.visibleCells().size()
1828             << " is not equal to the number of cells in the mesh "
1829             << mesh_.nCells() << abort(FatalError);
1830     }
1831
1832     if
1833     (
1834         cellLevel_.size() != mesh_.nCells()
1835      || pointLevel_.size() != mesh_.nPoints()
1836     )
1837     {
1838         FatalErrorIn
1839         (
1840             "hexRef8::hexRef8(const polyMesh&, const labelList&"
1841             ", const labelList&, const refinementHistory&)"
1842         )   << "Incorrect cellLevel or pointLevel size." << endl
1843             << "Number of cells in mesh:" << mesh_.nCells()
1844             << " does not equal size of cellLevel:" << cellLevel_.size() << endl
1845             << "Number of points in mesh:" << mesh_.nPoints()
1846             << " does not equal size of pointLevel:" << pointLevel_.size()
1847             << abort(FatalError);
1848     }
1849
1850     // Check refinement levels for consistency
1851     checkRefinementLevels(-1, labelList(0));
1852
1853
1854     // Check initial mesh for consistency
1855
1856     //if (debug)
1857     {
1858         checkMesh();
1859     }
1860 }
1861
1862
1863 // Construct from components
1864 Foam::hexRef8::hexRef8
1865 (
1866     const polyMesh& mesh,
1867     const labelList& cellLevel,
1868     const labelList& pointLevel
1869 )
1870 :
1871     mesh_(mesh),
1872     cellLevel_
1873     (
1874         IOobject
1875         (
1876             "cellLevel",
1877             mesh_.facesInstance(),
1878             polyMesh::meshSubDir,
1879             mesh_,
1880             IOobject::NO_READ,
1881             IOobject::AUTO_WRITE
1882         ),
1883         cellLevel
1884     ),
1885     pointLevel_
1886     (
1887         IOobject
1888         (
1889             "pointLevel",
1890             mesh_.facesInstance(),
1891             polyMesh::meshSubDir,
1892             mesh_,
1893             IOobject::NO_READ,
1894             IOobject::AUTO_WRITE
1895         ),
1896         pointLevel
1897     ),
1898     level0Edge_(getLevel0EdgeLength()),
1899     history_
1900     (
1901         IOobject
1902         (
1903             "refinementHistory",
1904             mesh_.facesInstance(),
1905             polyMesh::meshSubDir,
1906             mesh_,
1907             IOobject::NO_READ,
1908             IOobject::AUTO_WRITE
1909         ),
1910         List<refinementHistory::splitCell8>(0),
1911         labelList(0)
1912     ),
1913     faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
1914     savedPointLevel_(0),
1915     savedCellLevel_(0)
1916 {
1917     if
1918     (
1919         cellLevel_.size() != mesh_.nCells()
1920      || pointLevel_.size() != mesh_.nPoints()
1921     )
1922     {
1923         FatalErrorIn
1924         (
1925             "hexRef8::hexRef8(const polyMesh&, const labelList&"
1926             ", const labelList&)"
1927         )   << "Incorrect cellLevel or pointLevel size." << endl
1928             << "Number of cells in mesh:" << mesh_.nCells()
1929             << " does not equal size of cellLevel:" << cellLevel_.size() << endl
1930             << "Number of points in mesh:" << mesh_.nPoints()
1931             << " does not equal size of pointLevel:" << pointLevel_.size()
1932             << abort(FatalError);
1933     }
1934
1935     // Check refinement levels for consistency
1936     checkRefinementLevels(-1, labelList(0));
1937
1938     // Check initial mesh for consistency
1939
1940     //if (debug)
1941     {
1942         checkMesh();
1943     }
1944 }
1945
1946
1947 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1948
1949 //- Get points of a cell (without using cellPoints addressing)
1950 Foam::labelList Foam::hexRef8::cellPoints(const label cellI) const
1951 {
1952     // Pick up points of the cell
1953     const cell& cFaces = mesh_.cells()[cellI];
1954
1955     labelHashSet cPoints(4*cFaces.size());
1956
1957     forAll(cFaces, i)
1958     {
1959         const face& f = mesh_.faces()[cFaces[i]];
1960
1961         forAll(f, fp)
1962         {
1963             cPoints.insert(f[fp]);
1964         }
1965     }
1966     return cPoints.toc();
1967 }
1968
1969
1970 Foam::labelList Foam::hexRef8::consistentRefinement
1971 (
1972     const labelList& cellsToRefine,
1973     const bool maxSet
1974 ) const
1975 {
1976     // Loop, modifying cellsToRefine, until no more changes to due to 2:1
1977     // conflicts.
1978     // maxSet = false : unselect cells to refine
1979     // maxSet = true  : select cells to refine
1980
1981     // Go to straight boolList.
1982     PackedList<1> refineCell(mesh_.nCells(), 0);
1983     forAll(cellsToRefine, i)
1984     {
1985         refineCell.set(cellsToRefine[i], 1);
1986     }
1987
1988     while (true)
1989     {
1990         label nChanged = faceConsistentRefinement(maxSet, refineCell);
1991
1992         reduce(nChanged, sumOp<label>());
1993
1994         if (debug)
1995         {
1996             Pout<< "hexRef8::consistentRefinement : Changed " << nChanged
1997                 << " refinement levels due to 2:1 conflicts."
1998                 << endl;
1999         }
2000
2001         if (nChanged == 0)
2002         {
2003             break;
2004         }
2005     }
2006
2007
2008     // Convert back to labelList.
2009     label nRefined = 0;
2010
2011     forAll(refineCell, cellI)
2012     {
2013         if (refineCell.get(cellI) == 1)
2014         {
2015             nRefined++;
2016         }
2017     }
2018
2019     labelList newCellsToRefine(nRefined);
2020     nRefined = 0;
2021
2022     forAll(refineCell, cellI)
2023     {
2024         if (refineCell.get(cellI) == 1)
2025         {
2026             newCellsToRefine[nRefined++] = cellI;
2027         }
2028     }
2029
2030     if (debug)
2031     {
2032         checkWantedRefinementLevels(newCellsToRefine);
2033     }
2034
2035     return newCellsToRefine;
2036 }
2037
2038
2039 // Given a list of cells to refine determine additional cells to refine
2040 // such that the overall refinement:
2041 // - satisfies maxFaceDiff (e.g. 2:1) across neighbouring faces
2042 // - satisfies maxPointDiff (e.g. 4:1) across selected point connected
2043 //   cells. This is used to ensure that e.g. cells on the surface are not
2044 //   point connected to cells which are 8 times smaller.
2045 Foam::labelList Foam::hexRef8::consistentSlowRefinement
2046 (
2047     const label maxFaceDiff,
2048     const labelList& cellsToRefine,
2049     const labelList& facesToCheck,
2050     const label maxPointDiff,
2051     const labelList& pointsToCheck
2052 ) const
2053 {
2054     const labelList& faceOwner = mesh_.faceOwner();
2055     const labelList& faceNeighbour = mesh_.faceNeighbour();
2056
2057
2058     if (maxFaceDiff <= 0)
2059     {
2060         FatalErrorIn
2061         (
2062             "hexRef8::consistentSlowRefinement"
2063             "(const label, const labelList&, const labelList&"
2064             ", const label, const labelList&)"
2065         )   << "Illegal maxFaceDiff " << maxFaceDiff << nl
2066             << "Value should be >= 1" << exit(FatalError);
2067     }
2068
2069
2070     // Bit tricky. Say we want a distance of three cells between two
2071     // consecutive refinement levels. This is done by using FaceCellWave to
2072     // transport out the new refinement level. It gets decremented by one
2073     // every cell it crosses so if we initialize it to maxFaceDiff
2074     // we will get a field everywhere that tells us whether an unselected cell
2075     // needs refining as well.
2076
2077
2078     // Initial information about (distance to) cellLevel on all cells
2079     List<refinementData> allCellInfo(mesh_.nCells());
2080
2081     // Initial information about (distance to) cellLevel on all faces
2082     List<refinementData> allFaceInfo(mesh_.nFaces());
2083
2084     forAll(allCellInfo, cellI)
2085     {
2086         // maxFaceDiff since refinementData counts both
2087         // faces and cells.
2088         allCellInfo[cellI] = refinementData
2089         (
2090             maxFaceDiff*(cellLevel_[cellI]+1),// when cell is to be refined
2091             maxFaceDiff*cellLevel_[cellI]     // current level
2092         );
2093     }
2094
2095     // Cells to be refined will have cellLevel+1
2096     forAll(cellsToRefine, i)
2097     {
2098         label cellI = cellsToRefine[i];
2099
2100         allCellInfo[cellI].count() = allCellInfo[cellI].refinementCount();
2101     }
2102
2103
2104     // Labels of seed faces
2105     DynamicList<label> seedFaces(mesh_.nFaces()/100);
2106     // refinementLevel data on seed faces
2107     DynamicList<refinementData> seedFacesInfo(mesh_.nFaces()/100);
2108
2109
2110     // Additional buffer layer thickness by changing initial count. Usually
2111     // this happens on boundary faces. Bit tricky. Use allFaceInfo to mark
2112     // off thus marked faces so they're skipped in the next loop.
2113     forAll(facesToCheck, i)
2114     {
2115         label faceI = facesToCheck[i];
2116
2117         if (allFaceInfo[faceI].valid())
2118         {
2119             // Can only occur if face has already gone through loop below.
2120             FatalErrorIn
2121             (
2122                 "hexRef8::consistentSlowRefinement"
2123                 "(const label, const labelList&, const labelList&"
2124                 ", const label, const labelList&)"
2125             )   << "Argument facesToCheck seems to have duplicate entries!"
2126                 << endl
2127                 << "face:" << faceI << " occurs at positions "
2128                 << findIndices(facesToCheck, faceI)
2129                 << abort(FatalError);
2130         }
2131
2132
2133         const refinementData& ownData = allCellInfo[faceOwner[faceI]];
2134
2135         if (mesh_.isInternalFace(faceI))
2136         {
2137             // Seed face if neighbouring cell (after possible refinement)
2138             // will be refined one more than the current owner or neighbour.
2139
2140             const refinementData& neiData = allCellInfo[faceNeighbour[faceI]];
2141
2142             label faceCount;
2143             label faceRefineCount;
2144             if (neiData.count() > ownData.count())
2145             {
2146                 faceCount = neiData.count() + maxFaceDiff;
2147                 faceRefineCount = faceCount + maxFaceDiff;
2148             }
2149             else
2150             {
2151                 faceCount = ownData.count() + maxFaceDiff;
2152                 faceRefineCount = faceCount + maxFaceDiff;
2153             }
2154
2155             seedFaces.append(faceI);
2156             seedFacesInfo.append
2157             (
2158                 refinementData
2159                 (
2160                     faceRefineCount,
2161                     faceCount
2162                 )
2163             );
2164             allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2165         }
2166         else
2167         {
2168             label faceCount = ownData.count() + maxFaceDiff;
2169             label faceRefineCount = faceCount + maxFaceDiff;
2170
2171             seedFaces.append(faceI);
2172             seedFacesInfo.append
2173             (
2174                 refinementData
2175                 (
2176                     faceRefineCount,
2177                     faceCount
2178                 )
2179             );
2180             allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2181         }
2182     }
2183
2184
2185     // Just seed with all faces inbetween different refinement levels for now
2186     // (alternatively only seed faces on cellsToRefine but that gives problems
2187     //  if no cells to refine)
2188     forAll(faceNeighbour, faceI)
2189     {
2190         // Check if face already handled in loop above
2191         if (!allFaceInfo[faceI].valid())
2192         {
2193             label own = faceOwner[faceI];
2194             label nei = faceNeighbour[faceI];
2195
2196             // Seed face with transported data from highest cell.
2197
2198             if (allCellInfo[own].count() > allCellInfo[nei].count())
2199             {
2200                 allFaceInfo[faceI].updateFace
2201                 (
2202                     mesh_,
2203                     faceI,
2204                     own,
2205                     allCellInfo[own],
2206                     FaceCellWave<refinementData>::propagationTol()
2207                 );
2208                 seedFaces.append(faceI);
2209                 seedFacesInfo.append(allFaceInfo[faceI]);
2210             }
2211             else if (allCellInfo[own].count() < allCellInfo[nei].count())
2212             {
2213                 allFaceInfo[faceI].updateFace
2214                 (
2215                     mesh_,
2216                     faceI,
2217                     nei,
2218                     allCellInfo[nei],
2219                     FaceCellWave<refinementData>::propagationTol()
2220                 );
2221                 seedFaces.append(faceI);
2222                 seedFacesInfo.append(allFaceInfo[faceI]);
2223             }
2224         }
2225     }
2226
2227     // Seed all boundary faces with owner value. This is to make sure that
2228     // they are visited (probably only important for coupled faces since
2229     // these need to be visited from both sides)
2230     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2231     {
2232         // Check if face already handled in loop above
2233         if (!allFaceInfo[faceI].valid())
2234         {
2235             label own = faceOwner[faceI];
2236
2237             // Seed face with transported data from owner.
2238             refinementData faceData;
2239             faceData.updateFace
2240             (
2241                 mesh_,
2242                 faceI,
2243                 own,
2244                 allCellInfo[own],
2245                 FaceCellWave<refinementData>::propagationTol()
2246             );
2247             seedFaces.append(faceI);
2248             seedFacesInfo.append(faceData);
2249         }
2250     }
2251
2252
2253     // face-cell-face transport engine
2254     FaceCellWave<refinementData> levelCalc
2255     (
2256         mesh_,
2257         allFaceInfo,
2258         allCellInfo
2259     );
2260
2261     while(true)
2262     {
2263         if (debug)
2264         {
2265             Pout<< "hexRef8::consistentSlowRefinement : Seeded "
2266                 << seedFaces.size() << " faces between cells with different"
2267                 << " refinement level." << endl;
2268         }
2269
2270         // Set seed faces
2271         levelCalc.setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2272         seedFaces.clear();
2273         seedFacesInfo.clear();
2274
2275         // Iterate until no change. Now 2:1 face difference should be satisfied
2276         levelCalc.iterate(mesh_.globalData().nTotalFaces());
2277
2278
2279         // Now check point-connected cells (face-connected cells already ok):
2280         // - get per point max of connected cells
2281         // - sync across coupled points
2282         // - check cells against above point max
2283
2284         if (maxPointDiff == -1)
2285         {
2286             // No need to do any point checking.
2287             break;
2288         }
2289
2290         // Determine per point the max cell level. (done as count, not
2291         // as cell level purely for ease)
2292         labelList maxPointCount(mesh_.nPoints(), 0);
2293
2294         const labelListList& pointCells = mesh_.pointCells();
2295
2296         forAll(pointCells, pointI)
2297         {
2298             label& pLevel = maxPointCount[pointI];
2299
2300             const labelList& pCells = pointCells[pointI];
2301
2302             forAll(pCells, i)
2303             {
2304                 pLevel = max(pLevel, allCellInfo[pCells[i]].count());
2305             }
2306         }
2307
2308         // Sync maxPointCount to neighbour
2309         syncTools::syncPointList
2310         (
2311             mesh_,
2312             maxPointCount,
2313             maxEqOp<label>(),
2314             labelMin,           // null value
2315             false               // no separation
2316         );
2317
2318         // Update allFaceInfo from maxPointCount for all points to check
2319         // (usually on boundary faces)
2320
2321         // Per face the new refinement data
2322         Map<refinementData> changedFacesInfo(pointsToCheck.size());
2323
2324         forAll(pointsToCheck, i)
2325         {
2326             label pointI = pointsToCheck[i];
2327
2328             // Loop over all cells using the point and check whether their
2329             // refinement level is much less than the maximum.
2330
2331             const labelList& pCells = pointCells[pointI];
2332
2333             forAll(pCells, pCellI)
2334             {
2335                 label cellI = pCells[pCellI];
2336
2337                 refinementData& cellInfo = allCellInfo[cellI];
2338
2339                 if
2340                 (
2341                    !cellInfo.isRefined()
2342                  && (
2343                         maxPointCount[pointI]
2344                       > cellInfo.count() + maxFaceDiff*maxPointDiff
2345                     )
2346                 )
2347                 {
2348                     // Mark cell for refinement
2349                     cellInfo.count() = cellInfo.refinementCount();
2350
2351                     // Insert faces of cell as seed faces.
2352                     const cell& cFaces = mesh_.cells()[cellI];
2353
2354                     forAll(cFaces, cFaceI)
2355                     {
2356                         label faceI = cFaces[cFaceI];
2357
2358                         refinementData faceData;
2359                         faceData.updateFace
2360                         (
2361                             mesh_,
2362                             faceI,
2363                             cellI,
2364                             cellInfo,
2365                             FaceCellWave<refinementData>::propagationTol()
2366                         );
2367
2368                         if (faceData.count() > allFaceInfo[faceI].count())
2369                         {
2370                             changedFacesInfo.insert(faceI, faceData);
2371                         }
2372                     }
2373                 }
2374             }
2375         }
2376
2377         label nChanged = changedFacesInfo.size();
2378         reduce(nChanged, sumOp<label>());
2379
2380         if (nChanged == 0)
2381         {
2382             break;
2383         }
2384
2385
2386         // Transfer into seedFaces, seedFacesInfo
2387         seedFaces.setSize(changedFacesInfo.size());
2388         seedFacesInfo.setSize(changedFacesInfo.size());
2389
2390         forAllConstIter(Map<refinementData>, changedFacesInfo, iter)
2391         {
2392             seedFaces.append(iter.key());
2393             seedFacesInfo.append(iter());
2394         }
2395     }
2396
2397
2398     if (debug)
2399     {
2400         for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2401         {
2402             label own = mesh_.faceOwner()[faceI];
2403             label ownLevel =
2404                 cellLevel_[own]
2405               + (allCellInfo[own].isRefined() ? 1 : 0);
2406
2407             label nei = mesh_.faceNeighbour()[faceI];
2408             label neiLevel =
2409                 cellLevel_[nei]
2410               + (allCellInfo[nei].isRefined() ? 1 : 0);
2411
2412             if (mag(ownLevel-neiLevel) > 1)
2413             {
2414                 FatalErrorIn
2415                 (
2416                     "hexRef8::consistentSlowRefinement"
2417                 )   << "cell:" << own
2418                     << " current level:" << cellLevel_[own]
2419                     << " current refData:" << allCellInfo[own]
2420                     << " level after refinement:" << ownLevel
2421                     << nl
2422                     << "neighbour cell:" << nei
2423                     << " current level:" << cellLevel_[nei]
2424                     << " current refData:" << allCellInfo[nei]
2425                     << " level after refinement:" << neiLevel
2426                     << nl
2427                     << "which does not satisfy 2:1 constraints anymore." << nl
2428                     << "face:" << faceI << " faceRefData:" << allFaceInfo[faceI]
2429                     << abort(FatalError);
2430             }
2431         }
2432
2433
2434         // Coupled faces. Swap owner level to get neighbouring cell level.
2435         // (only boundary faces of neiLevel used)
2436
2437         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2438         labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
2439         labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
2440
2441         forAll(neiLevel, i)
2442         {
2443             label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2444             neiLevel[i] = cellLevel_[own];
2445             neiCount[i] = allCellInfo[own].count();
2446             neiRefCount[i] = allCellInfo[own].refinementCount();
2447         }
2448
2449         // Swap to neighbour
2450         syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
2451         syncTools::swapBoundaryFaceList(mesh_, neiCount, false);
2452         syncTools::swapBoundaryFaceList(mesh_, neiRefCount, false);
2453
2454         // Now we have neighbour value see which cells need refinement
2455         forAll(neiLevel, i)
2456         {
2457             label faceI = i+mesh_.nInternalFaces();
2458
2459             label own = mesh_.faceOwner()[faceI];
2460             label ownLevel =
2461                 cellLevel_[own]
2462               + (allCellInfo[own].isRefined() ? 1 : 0);
2463
2464             label nbrLevel =
2465                 neiLevel[i]
2466               + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2467
2468             if (mag(ownLevel - nbrLevel) > 1)
2469             {
2470                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
2471
2472                 FatalErrorIn
2473                 (
2474                     "hexRef8::consistentSlowRefinement"
2475                     "(const label, const labelList&, const labelList&"
2476                     ", const label, const labelList&)"
2477                 )   << "Celllevel does not satisfy 2:1 constraint."
2478                     << " On coupled face "
2479                     << faceI
2480                     << " refData:" << allFaceInfo[faceI]
2481                     << " on patch " << patchI << " "
2482                     << mesh_.boundaryMesh()[patchI].name() << nl
2483                     << "owner cell " << own
2484                     << " current level:" << cellLevel_[own]
2485                     << " current count:" << allCellInfo[own].count()
2486                     << " current refCount:"
2487                     << allCellInfo[own].refinementCount()
2488                     << " level after refinement:" << ownLevel
2489                     << nl
2490                     << "(coupled) neighbour cell"
2491                     << " has current level:" << neiLevel[i]
2492                     << " current count:" << neiCount[i]
2493                     << " current refCount:" << neiRefCount[i]
2494                     << " level after refinement:" << nbrLevel
2495                     << abort(FatalError);
2496             }
2497         }
2498     }
2499
2500     // Convert back to labelList of cells to refine.
2501
2502     label nRefined = 0;
2503
2504     forAll(allCellInfo, cellI)
2505     {
2506         if (allCellInfo[cellI].isRefined())
2507         {
2508             nRefined++;
2509         }
2510     }
2511
2512     // Updated list of cells to refine
2513     labelList newCellsToRefine(nRefined);
2514     nRefined = 0;
2515
2516     forAll(allCellInfo, cellI)
2517     {
2518         if (allCellInfo[cellI].isRefined())
2519         {
2520             newCellsToRefine[nRefined++] = cellI;
2521         }
2522     }
2523
2524     if (debug)
2525     {
2526         Pout<< "hexRef8::consistentSlowRefinement : From "
2527             << cellsToRefine.size() << " to " << newCellsToRefine.size()
2528             << " cells to refine." << endl;
2529     }
2530
2531     return newCellsToRefine;
2532 }
2533
2534
2535 Foam::labelList Foam::hexRef8::consistentSlowRefinement2
2536 (
2537     const label maxFaceDiff,
2538     const labelList& cellsToRefine,
2539     const labelList& facesToCheck
2540 ) const
2541 {
2542     const labelList& faceOwner = mesh_.faceOwner();
2543     const labelList& faceNeighbour = mesh_.faceNeighbour();
2544
2545     if (maxFaceDiff <= 0)
2546     {
2547         FatalErrorIn
2548         (
2549             "hexRef8::consistentSlowRefinement2"
2550             "(const label, const labelList&, const labelList&)"
2551         )   << "Illegal maxFaceDiff " << maxFaceDiff << nl
2552             << "Value should be >= 1" << exit(FatalError);
2553     }
2554
2555     const scalar level0Size = 2*maxFaceDiff*level0Edge_;
2556
2557
2558     // Bit tricky. Say we want a distance of three cells between two
2559     // consecutive refinement levels. This is done by using FaceCellWave to
2560     // transport out the 'refinement shell'. Anything inside the refinement
2561     // shell (given by a distance) gets marked for refinement.
2562
2563     // Initial information about (distance to) cellLevel on all cells
2564     List<refinementDistanceData> allCellInfo(mesh_.nCells());
2565
2566     // Initial information about (distance to) cellLevel on all faces
2567     List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
2568
2569
2570     // Mark cells with wanted refinement level
2571     forAll(cellsToRefine, i)
2572     {
2573         label cellI = cellsToRefine[i];
2574
2575         allCellInfo[cellI] = refinementDistanceData
2576         (
2577             level0Size,
2578             mesh_.cellCentres()[cellI],
2579             cellLevel_[cellI]+1             // wanted refinement
2580         );
2581     }
2582     // Mark all others with existing refinement level
2583     forAll(allCellInfo, cellI)
2584     {
2585         if (!allCellInfo[cellI].valid())
2586         {
2587             allCellInfo[cellI] = refinementDistanceData
2588             (
2589                 level0Size,
2590                 mesh_.cellCentres()[cellI],
2591                 cellLevel_[cellI]           // wanted refinement
2592             );
2593         }
2594     }
2595
2596
2597     // Labels of seed faces
2598     DynamicList<label> seedFaces(mesh_.nFaces()/100);
2599     // refinementLevel data on seed faces
2600     DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
2601
2602
2603     const pointField& cc = mesh_.cellCentres();
2604
2605     forAll(facesToCheck, i)
2606     {
2607         label faceI = facesToCheck[i];
2608
2609         if (allFaceInfo[faceI].valid())
2610         {
2611             // Can only occur if face has already gone through loop below.
2612             FatalErrorIn
2613             (
2614                 "hexRef8::consistentSlowRefinement2"
2615                 "(const label, const labelList&, const labelList&)"
2616             )   << "Argument facesToCheck seems to have duplicate entries!"
2617                 << endl
2618                 << "face:" << faceI << " occurs at positions "
2619                 << findIndices(facesToCheck, faceI)
2620                 << abort(FatalError);
2621         }
2622
2623         label own = faceOwner[faceI];
2624
2625         label ownLevel =
2626         (
2627             allCellInfo[own].valid()
2628           ? allCellInfo[own].originLevel()
2629           : cellLevel_[own]
2630         );
2631
2632         if (!mesh_.isInternalFace(faceI))
2633         {
2634             // Do as if boundary face would have neighbour with one higher
2635             // refinement level.
2636             const point& fc = mesh_.faceCentres()[faceI];
2637
2638             refinementDistanceData neiData
2639             (
2640                 level0Size,
2641                 2*fc - cc[own],    // est'd cell centre
2642                 ownLevel+1
2643             );
2644
2645             allFaceInfo[faceI].updateFace
2646             (
2647                 mesh_,
2648                 faceI,
2649                 own,        // not used (should be nei)
2650                 neiData,
2651                 FaceCellWave<refinementDistanceData>::propagationTol()
2652             );
2653         }
2654         else
2655         {
2656             label nei = faceNeighbour[faceI];
2657
2658             label neiLevel =
2659             (
2660                 allCellInfo[nei].valid()
2661               ? allCellInfo[nei].originLevel()
2662               : cellLevel_[nei]
2663             );
2664
2665             if (ownLevel == neiLevel)
2666             {
2667                 // Fake as if nei>own or own>nei (whichever one 'wins')
2668                 allFaceInfo[faceI].updateFace
2669                 (
2670                     mesh_,
2671                     faceI,
2672                     nei,
2673                     refinementDistanceData(level0Size, cc[nei], neiLevel+1),
2674                     FaceCellWave<refinementDistanceData>::propagationTol()
2675                 );
2676                 allFaceInfo[faceI].updateFace
2677                 (
2678                     mesh_,
2679                     faceI,
2680                     own,
2681                     refinementDistanceData(level0Size, cc[own], ownLevel+1),
2682                     FaceCellWave<refinementDistanceData>::propagationTol()
2683                 );
2684             }
2685             else
2686             {
2687                 // Difference in level anyway.
2688                 allFaceInfo[faceI].updateFace
2689                 (
2690                     mesh_,
2691                     faceI,
2692                     nei,
2693                     refinementDistanceData(level0Size, cc[nei], neiLevel),
2694                     FaceCellWave<refinementDistanceData>::propagationTol()
2695                 );
2696                 allFaceInfo[faceI].updateFace
2697                 (
2698                     mesh_,
2699                     faceI,
2700                     own,
2701                     refinementDistanceData(level0Size, cc[own], ownLevel),
2702                     FaceCellWave<refinementDistanceData>::propagationTol()
2703                 );
2704             }
2705         }
2706         seedFaces.append(faceI);
2707         seedFacesInfo.append(allFaceInfo[faceI]);
2708     }
2709
2710
2711     // Create some initial seeds to start walking from. This is only if there
2712     // are no facesToCheck.
2713     // Just seed with all faces inbetween different refinement levels for now
2714     forAll(faceNeighbour, faceI)
2715     {
2716         // Check if face already handled in loop above
2717         if (!allFaceInfo[faceI].valid())
2718         {
2719             label own = faceOwner[faceI];
2720
2721             label ownLevel =
2722             (
2723                 allCellInfo[own].valid()
2724               ? allCellInfo[own].originLevel()
2725               : cellLevel_[own]
2726             );
2727
2728             label nei = faceNeighbour[faceI];
2729
2730             label neiLevel =
2731             (
2732                 allCellInfo[nei].valid()
2733               ? allCellInfo[nei].originLevel()
2734               : cellLevel_[nei]
2735             );
2736
2737             if (ownLevel > neiLevel)
2738             {
2739                 // Set face to owner data. (since face not yet would be copy)
2740                 seedFaces.append(faceI);
2741                 allFaceInfo[faceI].updateFace
2742                 (
2743                     mesh_,
2744                     faceI,
2745                     own,
2746                     refinementDistanceData(level0Size, cc[own], ownLevel),
2747                     FaceCellWave<refinementDistanceData>::propagationTol()
2748                 );
2749                 seedFacesInfo.append(allFaceInfo[faceI]);
2750             }
2751             else if (neiLevel > ownLevel)
2752             {
2753                 seedFaces.append(faceI);
2754                 allFaceInfo[faceI].updateFace
2755                 (
2756                     mesh_,
2757                     faceI,
2758                     nei,
2759                     refinementDistanceData(level0Size, cc[nei], neiLevel),
2760                     FaceCellWave<refinementDistanceData>::propagationTol()
2761                 );
2762                 seedFacesInfo.append(allFaceInfo[faceI]);
2763             }
2764         }
2765     }
2766
2767     seedFaces.shrink();
2768     seedFacesInfo.shrink();
2769
2770     // face-cell-face transport engine
2771     FaceCellWave<refinementDistanceData> levelCalc
2772     (
2773         mesh_,
2774         seedFaces,
2775         seedFacesInfo,
2776         allFaceInfo,
2777         allCellInfo,
2778         mesh_.globalData().nTotalCells()
2779     );
2780
2781
2782     //if (debug)
2783     //{
2784     //    // Dump wanted level
2785     //    volScalarField wantedLevel
2786     //    (
2787     //        IOobject
2788     //        (
2789     //            "wantedLevel",
2790     //            fMesh.time().timeName(),
2791     //            fMesh,
2792     //            IOobject::NO_READ,
2793     //            IOobject::AUTO_WRITE,
2794     //            false
2795     //        ),
2796     //        fMesh,
2797     //        dimensionedScalar("zero", dimless, 0)
2798     //    );
2799     //
2800     //    forAll(wantedLevel, cellI)
2801     //    {
2802     //        wantedLevel[cellI] = allCellInfo[cellI].wantedLevel(cc[cellI]);
2803     //    }
2804     //
2805     //    Pout<< "Writing " << wantedLevel.objectPath() << endl;
2806     //    wantedLevel.write();
2807     //}
2808
2809
2810     // Convert back to labelList of cells to refine.
2811     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2812
2813     // 1. Force original refinement cells to be picked up by setting the
2814     // originLevel of input cells to be a very large level (but within range
2815     // of 1<< shift inside refinementDistanceData::wantedLevel)
2816     forAll(cellsToRefine, i)
2817     {
2818         label cellI = cellsToRefine[i];
2819
2820         allCellInfo[cellI].originLevel() = sizeof(label)*8-2;
2821         allCellInfo[cellI].origin() = cc[cellI];
2822     }
2823
2824     // 2. Extend to 2:1. For non-cube cells the scalar distance does not work
2825     // so make sure it at least provides 2:1.
2826     PackedList<1> refineCell(mesh_.nCells(), 0);
2827     forAll(allCellInfo, cellI)
2828     {
2829         label wanted = allCellInfo[cellI].wantedLevel(cc[cellI]);
2830
2831         if (wanted > cellLevel_[cellI]+1)
2832         {
2833             refineCell.set(cellI, 1);
2834         }
2835     }
2836
2837     while (true)
2838     {
2839         label nChanged = faceConsistentRefinement(true, refineCell);
2840
2841         reduce(nChanged, sumOp<label>());
2842
2843         if (debug)
2844         {
2845             Pout<< "hexRef8::consistentSlowRefinement2 : Changed " << nChanged
2846                 << " refinement levels due to 2:1 conflicts."
2847                 << endl;
2848         }
2849
2850         if (nChanged == 0)
2851         {
2852             break;
2853         }
2854     }
2855
2856     // 3. Convert back to labelList.
2857     label nRefined = 0;
2858
2859     forAll(refineCell, cellI)
2860     {
2861         if (refineCell.get(cellI) == 1)
2862         {
2863             nRefined++;
2864         }
2865     }
2866
2867     labelList newCellsToRefine(nRefined);
2868     nRefined = 0;
2869
2870     forAll(refineCell, cellI)
2871     {
2872         if (refineCell.get(cellI) == 1)
2873         {
2874             newCellsToRefine[nRefined++] = cellI;
2875         }
2876     }
2877
2878     if (debug)
2879     {
2880         Pout<< "hexRef8::consistentSlowRefinement2 : From "
2881             << cellsToRefine.size() << " to " << newCellsToRefine.size()
2882             << " cells to refine." << endl;
2883
2884         // Check that newCellsToRefine obeys at least 2:1.
2885
2886         {
2887             cellSet cellsIn(mesh_, "cellsToRefineIn", cellsToRefine);
2888             Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2889                 << cellsIn.size() << " to cellSet "
2890                 << cellsIn.objectPath() << endl;
2891             cellsIn.write();
2892         }
2893         {
2894             cellSet cellsOut(mesh_, "cellsToRefineOut", newCellsToRefine);
2895             Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2896                 << cellsOut.size() << " to cellSet "
2897                 << cellsOut.objectPath() << endl;
2898             cellsOut.write();
2899         }
2900
2901         // Extend to 2:1
2902         PackedList<1> refineCell(mesh_.nCells(), 0);
2903         forAll(newCellsToRefine, i)
2904         {
2905             refineCell.set(newCellsToRefine[i], 1);
2906         }
2907         const PackedList<1> savedRefineCell(refineCell);
2908
2909         label nChanged = faceConsistentRefinement(true, refineCell);
2910
2911         {
2912             cellSet cellsOut2
2913             (
2914                 mesh_, "cellsToRefineOut2", newCellsToRefine.size()
2915             );
2916             forAll(refineCell, cellI)
2917             {
2918                 if (refineCell.get(cellI) == 1)
2919                 {
2920                     cellsOut2.insert(cellI);
2921                 }
2922             }
2923             Pout<< "hexRef8::consistentSlowRefinement2 : writing "
2924                 << cellsOut2.size() << " to cellSet "
2925                 << cellsOut2.objectPath() << endl;
2926             cellsOut2.write();
2927         }
2928
2929         if (nChanged > 0)
2930         {
2931             forAll(refineCell, cellI)
2932             {
2933                 if
2934                 (
2935                     refineCell.get(cellI) == 1
2936                  && savedRefineCell.get(cellI) == 0
2937                 )
2938                 {
2939                     FatalErrorIn
2940                     (
2941                         "hexRef8::consistentSlowRefinement2"
2942                         "(const label, const labelList&, const labelList&)"
2943                     )   << "Cell:" << cellI << " cc:"
2944                         << mesh_.cellCentres()[cellI]
2945                         << " was not marked for refinement but does not obey"
2946                         << " 2:1 constraints."
2947                         << abort(FatalError);
2948                 }
2949             }
2950         }
2951     }
2952
2953     return newCellsToRefine;
2954 }
2955
2956
2957 // Top level driver to insert topo changes to do all refinement.
2958 Foam::labelListList Foam::hexRef8::setRefinement
2959 (
2960     const labelList& cellLabels,
2961     polyTopoChange& meshMod
2962 )
2963 {
2964     if (debug)
2965     {
2966         Pout<< "hexRef8::setRefinement :"
2967             << " Checking initial mesh just to make sure" << endl;
2968
2969         checkMesh();
2970         checkRefinementLevels(-1, labelList(0));
2971     }
2972
2973     // Clear any saved point/cell data.
2974     savedPointLevel_.clear();
2975     savedCellLevel_.clear();
2976
2977
2978     // New point/cell level. Copy of pointLevel for existing points.
2979     DynamicList<label> newCellLevel(cellLevel_.size());
2980     forAll(cellLevel_, cellI)
2981     {
2982         newCellLevel.append(cellLevel_[cellI]);
2983     }
2984     DynamicList<label> newPointLevel(pointLevel_.size());
2985     forAll(pointLevel_, pointI)
2986     {
2987         newPointLevel.append(pointLevel_[pointI]);
2988     }
2989
2990
2991     if (debug)
2992     {
2993         Pout<< "hexRef8::setRefinement :"
2994             << " Allocating " << cellLabels.size() << " cell midpoints."
2995             << endl;
2996     }
2997
2998
2999     // Mid point per refined cell.
3000     // -1 : not refined
3001     // >=0: label of mid point.
3002     labelList cellMidPoint(mesh_.nCells(), -1);
3003
3004     forAll(cellLabels, i)
3005     {
3006         label cellI = cellLabels[i];
3007
3008         label anchorPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0];
3009
3010         cellMidPoint[cellI] = meshMod.setAction
3011         (
3012             polyAddPoint
3013             (
3014                 mesh_.cellCentres()[cellI],     // point
3015                 anchorPointI,                   // master point
3016                 -1,                             // zone for point
3017                 true                            // supports a cell
3018             )
3019         );
3020
3021         newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI]+1;
3022     }
3023
3024
3025     if (debug)
3026     {
3027         cellSet splitCells(mesh_, "splitCells", cellLabels.size());
3028
3029         forAll(cellMidPoint, cellI)
3030         {
3031             if (cellMidPoint[cellI] >= 0)
3032             {
3033                 splitCells.insert(cellI);
3034             }
3035         }
3036
3037         Pout<< "hexRef8::setRefinement : Dumping " << splitCells.size()
3038             << " cells to split to cellSet " << splitCells.objectPath()
3039             << endl;
3040
3041         splitCells.write();
3042     }
3043
3044
3045
3046     // Split edges
3047     // ~~~~~~~~~~~
3048
3049     if (debug)
3050     {
3051         Pout<< "hexRef8::setRefinement :"
3052             << " Allocating edge midpoints."
3053             << endl;
3054     }
3055
3056     // Unrefined edges are ones between cellLevel or lower points.
3057     // If any cell using this edge gets split then the edge needs to be split.
3058
3059     // -1  : no need to split edge
3060     // >=0 : label of introduced mid point
3061     labelList edgeMidPoint(mesh_.nEdges(), -1);
3062
3063     // Note: Loop over cells to be refined or edges?
3064     forAll(cellMidPoint, cellI)
3065     {
3066         if (cellMidPoint[cellI] >= 0)
3067         {
3068             const labelList& cEdges = mesh_.cellEdges()[cellI];
3069
3070             forAll(cEdges, i)
3071             {
3072                 label edgeI = cEdges[i];
3073
3074                 const edge& e = mesh_.edges()[edgeI];
3075
3076                 if
3077                 (
3078                     pointLevel_[e[0]] <= cellLevel_[cellI]
3079                  && pointLevel_[e[1]] <= cellLevel_[cellI]
3080                 )
3081                 {
3082                     edgeMidPoint[edgeI] = 12345;    // mark need for splitting
3083                 }
3084             }
3085         }
3086     }
3087
3088     // Synchronize edgeMidPoint across coupled patches. Take max so that
3089     // any split takes precedence.
3090     syncTools::syncEdgeList
3091     (
3092         mesh_,
3093         edgeMidPoint,
3094         maxEqOp<label>(),
3095         labelMin,
3096         false               // no separation
3097     );
3098
3099
3100     // Introduce edge points
3101     // ~~~~~~~~~~~~~~~~~~~~~
3102
3103     {
3104         // Phase 1: calculate midpoints and sync.
3105         // This needs doing for if people do not write binary and we slowly
3106         // get differences.
3107
3108         pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT));
3109
3110         forAll(edgeMidPoint, edgeI)
3111         {
3112             if (edgeMidPoint[edgeI] >= 0)
3113             {
3114                 // Edge marked to be split.
3115                 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3116             }
3117         }
3118         syncTools::syncEdgeList
3119         (
3120             mesh_,
3121             edgeMids,
3122             maxEqOp<vector>(),
3123             point(-GREAT, -GREAT, -GREAT),
3124             true               // apply separation
3125         );
3126
3127
3128         // Phase 2: introduce points at the synced locations.
3129         forAll(edgeMidPoint, edgeI)
3130         {
3131             if (edgeMidPoint[edgeI] >= 0)
3132             {
3133                 // Edge marked to be split. Replace edgeMidPoint with actual
3134                 // point label.
3135
3136                 const edge& e = mesh_.edges()[edgeI];
3137
3138                 edgeMidPoint[edgeI] = meshMod.setAction
3139                 (
3140                     polyAddPoint
3141                     (
3142                         edgeMids[edgeI],            // point
3143                         e[0],                       // master point
3144                         -1,                         // zone for point
3145                         true                        // supports a cell
3146                     )
3147                 );
3148
3149                 newPointLevel(edgeMidPoint[edgeI]) =
3150                     max
3151                     (
3152                         pointLevel_[e[0]],
3153                         pointLevel_[e[1]]
3154                     )
3155                   + 1;
3156             }
3157         }
3158     }
3159
3160     if (debug)
3161     {
3162         OFstream str(mesh_.time().path()/"edgeMidPoint.obj");
3163
3164         forAll(edgeMidPoint, edgeI)
3165         {
3166             if (edgeMidPoint[edgeI] >= 0)
3167             {
3168                 const edge& e = mesh_.edges()[edgeI];
3169
3170                 meshTools::writeOBJ(str, e.centre(mesh_.points()));
3171             }
3172         }
3173
3174         Pout<< "hexRef8::setRefinement :"
3175             << " Dumping edge centres to split to file " << str.name() << endl;
3176     }
3177
3178
3179     // Calculate face level
3180     // ~~~~~~~~~~~~~~~~~~~~
3181     // (after splitting)
3182
3183     if (debug)
3184     {
3185         Pout<< "hexRef8::setRefinement :"
3186             << " Allocating face midpoints."
3187             << endl;
3188     }
3189
3190     // Face anchor level. There are guaranteed 4 points with level
3191     // <= anchorLevel. These are the corner points.
3192     labelList faceAnchorLevel(mesh_.nFaces());
3193
3194     for (label faceI = 0; faceI < mesh_.nFaces(); faceI++)
3195     {
3196         faceAnchorLevel[faceI] = getAnchorLevel(faceI);
3197     }
3198
3199     // -1  : no need to split face
3200     // >=0 : label of introduced mid point
3201     labelList faceMidPoint(mesh_.nFaces(), -1);
3202
3203
3204     // Internal faces: look at cells on both sides. Uniquely determined since
3205     // face itself guaranteed to be same level as most refined neighbour.
3206     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3207     {
3208         if (faceAnchorLevel[faceI] >= 0)
3209         {
3210             label own = mesh_.faceOwner()[faceI];
3211             label ownLevel = cellLevel_[own];
3212             label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3213
3214             label nei = mesh_.faceNeighbour()[faceI];
3215             label neiLevel = cellLevel_[nei];
3216             label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3217
3218             if
3219             (
3220                 newOwnLevel > faceAnchorLevel[faceI]
3221              || newNeiLevel > faceAnchorLevel[faceI]
3222             )
3223             {
3224                 faceMidPoint[faceI] = 12345;    // mark to be split
3225             }
3226         }
3227     }
3228
3229     // Coupled patches handled like internal faces except now all information
3230     // from neighbour comes from across processor.
3231     // Boundary faces are more complicated since the boundary face can
3232     // be more refined than its owner (or neighbour for coupled patches)
3233     // (does not happen if refining/unrefining only, but does e.g. when
3234     //  refinining and subsetting)
3235
3236     {
3237         labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
3238
3239         forAll(newNeiLevel, i)
3240         {
3241             label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3242             label ownLevel = cellLevel_[own];
3243             label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3244
3245             newNeiLevel[i] = newOwnLevel;
3246         }
3247
3248         // Swap.
3249         syncTools::swapBoundaryFaceList(mesh_, newNeiLevel, false);
3250
3251         // So now we have information on the neighbour.
3252
3253         forAll(newNeiLevel, i)
3254         {
3255             label faceI = i+mesh_.nInternalFaces();
3256
3257             if (faceAnchorLevel[faceI] >= 0)
3258             {
3259                 label own = mesh_.faceOwner()[faceI];
3260                 label ownLevel = cellLevel_[own];
3261                 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3262
3263                 if
3264                 (
3265                     newOwnLevel > faceAnchorLevel[faceI]
3266                  || newNeiLevel[i] > faceAnchorLevel[faceI]
3267                 )
3268                 {
3269                     faceMidPoint[faceI] = 12345;    // mark to be split
3270                 }
3271             }
3272         }
3273     }
3274
3275
3276     // Synchronize faceMidPoint across coupled patches. (logical or)
3277     syncTools::syncFaceList
3278     (
3279         mesh_,
3280         faceMidPoint,
3281         maxEqOp<label>(),
3282         false
3283     );
3284
3285
3286
3287     // Introduce face points
3288     // ~~~~~~~~~~~~~~~~~~~~~
3289
3290     {
3291         // Phase 1: determine mid points and sync. See comment for edgeMids
3292         // above
3293         pointField bFaceMids
3294         (
3295             mesh_.nFaces()-mesh_.nInternalFaces(),
3296             point(-GREAT, -GREAT, -GREAT)
3297         );
3298
3299         forAll(bFaceMids, i)
3300         {
3301             label faceI = i+mesh_.nInternalFaces();
3302
3303             if (faceMidPoint[faceI] >= 0)
3304             {
3305                 bFaceMids[i] = mesh_.faceCentres()[faceI];
3306             }
3307         }
3308         syncTools::syncBoundaryFaceList
3309         (
3310             mesh_,
3311             bFaceMids,
3312             maxEqOp<vector>(),
3313             true               // apply separation
3314         );
3315
3316         forAll(faceMidPoint, faceI)
3317         {
3318             if (faceMidPoint[faceI] >= 0)
3319             {
3320                 // Face marked to be split. Replace faceMidPoint with actual
3321                 // point label.
3322
3323                 const face& f = mesh_.faces()[faceI];
3324
3325                 faceMidPoint[faceI] = meshMod.setAction
3326                 (
3327                     polyAddPoint
3328                     (
3329                         (
3330                             faceI < mesh_.nInternalFaces()
3331                           ? mesh_.faceCentres()[faceI]
3332                           : bFaceMids[faceI-mesh_.nInternalFaces()]
3333                         ),                          // point
3334                         f[0],                       // master point
3335                         -1,                         // zone for point
3336                         true                        // supports a cell
3337                     )
3338                 );
3339
3340                 // Determine the level of the corner points and midpoint will
3341                 // be one higher.
3342                 newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI]+1;
3343             }
3344         }
3345     }
3346
3347     if (debug)
3348     {
3349         faceSet splitFaces(mesh_, "splitFaces", cellLabels.size());
3350
3351         forAll(faceMidPoint, faceI)
3352         {
3353             if (faceMidPoint[faceI] >= 0)
3354             {
3355                 splitFaces.insert(faceI);
3356             }
3357         }
3358
3359         Pout<< "hexRef8::setRefinement : Dumping " << splitFaces.size()
3360             << " faces to split to faceSet " << splitFaces.objectPath() << endl;
3361
3362         splitFaces.write();
3363     }
3364
3365
3366     // Information complete
3367     // ~~~~~~~~~~~~~~~~~~~~
3368     // At this point we have all the information we need. We should no
3369     // longer reference the cellLabels to refine. All the information is:
3370     // - cellMidPoint >= 0 : cell needs to be split
3371     // - faceMidPoint >= 0 : face needs to be split
3372     // - edgeMidPoint >= 0 : edge needs to be split
3373
3374
3375
3376     // Get the corner/anchor points
3377     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3378
3379     if (debug)
3380     {
3381         Pout<< "hexRef8::setRefinement :"
3382             << " Finding cell anchorPoints (8 per cell)"
3383             << endl;
3384     }
3385
3386     // There will always be 8 points on the hex that have were introduced
3387     // with the hex and will have the same or lower refinement level.
3388
3389     // Per cell the 8 corner points.
3390     labelListList cellAnchorPoints(mesh_.nCells());
3391
3392     {
3393         labelList nAnchorPoints(mesh_.nCells(), 0);
3394
3395         forAll(cellMidPoint, cellI)
3396         {
3397             if (cellMidPoint[cellI] >= 0)
3398             {
3399                 cellAnchorPoints[cellI].setSize(8);
3400             }
3401         }
3402
3403         forAll(pointLevel_, pointI)
3404         {
3405             const labelList& pCells = mesh_.pointCells()[pointI];
3406
3407             forAll(pCells, pCellI)
3408             {
3409                 label cellI = pCells[pCellI];
3410
3411                 if
3412                 (
3413                     cellMidPoint[cellI] >= 0
3414                  && pointLevel_[pointI] <= cellLevel_[cellI]
3415                 )
3416                 {
3417                     if (nAnchorPoints[cellI] == 8)
3418                     {
3419                         FatalErrorIn
3420                         (
3421                             "hexRef8::setRefinement(const labelList&"
3422                             ", polyTopoChange&)"
3423                         )   << "cell " << cellI
3424                             << " of level " << cellLevel_[cellI]
3425                             << " uses more than 8 points of equal or"
3426                             << " lower level" << nl
3427                             << "Points so far:" << cellAnchorPoints[cellI]
3428                             << abort(FatalError);
3429                     }
3430
3431                     cellAnchorPoints[cellI][nAnchorPoints[cellI]++] = pointI;
3432                 }
3433             }
3434         }
3435
3436         forAll(cellMidPoint, cellI)
3437         {
3438             if (cellMidPoint[cellI] >= 0)
3439             {
3440                 if (nAnchorPoints[cellI] != 8)
3441                 {
3442                     const labelList cPoints(cellPoints(cellI));
3443
3444                     FatalErrorIn
3445                     (
3446                         "hexRef8::setRefinement(const labelList&"
3447                         ", polyTopoChange&)"
3448                     )   << "cell " << cellI
3449                         << " of level " << cellLevel_[cellI]
3450                         << " does not seem to have 8 points of equal or"
3451                         << " lower level" << endl
3452                         << "cellPoints:" << cPoints << endl
3453                         << "pointLevels:"
3454                         << IndirectList<label>(pointLevel_, cPoints)() << endl
3455                         << abort(FatalError);
3456                 }
3457             }
3458         }
3459     }
3460
3461
3462     // Add the cells
3463     // ~~~~~~~~~~~~~
3464
3465     if (debug)
3466     {
3467         Pout<< "hexRef8::setRefinement :"
3468             << " Adding cells (1 per anchorPoint)"
3469             << endl;
3470     }
3471
3472     // Per cell the 7 added cells (+ original cell)
3473     labelListList cellAddedCells(mesh_.nCells());
3474
3475     forAll(cellAnchorPoints, cellI)
3476     {
3477         const labelList& cAnchors = cellAnchorPoints[cellI];
3478
3479         if (cAnchors.size() == 8)
3480         {
3481             labelList& cAdded = cellAddedCells[cellI];
3482             cAdded.setSize(8);
3483
3484             // Original cell at 0
3485             cAdded[0] = cellI;
3486
3487             // Update cell level
3488             newCellLevel[cellI] = cellLevel_[cellI]+1;
3489
3490
3491             for (label i = 1; i < 8; i++)
3492             {
3493                 cAdded[i] = meshMod.setAction
3494                 (
3495                     polyAddCell
3496                     (
3497                         -1,                                 // master point
3498                         -1,                                 // master edge
3499                         -1,                                 // master face
3500                         cellI,                              // master cell
3501                         mesh_.cellZones().whichZone(cellI)  // zone for cell
3502                     )
3503                 );
3504
3505                 newCellLevel(cAdded[i]) = cellLevel_[cellI]+1;
3506             }
3507         }
3508     }
3509
3510
3511     // Faces
3512     // ~~~~~
3513     // 1. existing faces that get split (into four always)
3514     // 2. existing faces that do not get split but only edges get split
3515     // 3. existing faces that do not get split but get new owner/neighbour
3516     // 4. new internal faces inside split cells.
3517
3518     if (debug)
3519     {
3520         Pout<< "hexRef8::setRefinement :"
3521             << " Marking faces to be handled"
3522             << endl;
3523     }
3524
3525     // Get all affected faces.
3526     PackedList<1> affectedFace(mesh_.nFaces(), 0);
3527
3528     {
3529         forAll(cellMidPoint, cellI)
3530         {
3531             if (cellMidPoint[cellI] >= 0)
3532             {
3533                 const cell& cFaces = mesh_.cells()[cellI];
3534
3535                 forAll(cFaces, i)
3536                 {
3537                     affectedFace.set(cFaces[i], 1);
3538                 }
3539             }
3540         }
3541
3542         forAll(faceMidPoint, faceI)
3543         {
3544             if (faceMidPoint[faceI] >= 0)
3545             {
3546                 affectedFace.set(faceI, 1);
3547             }
3548         }
3549
3550         forAll(edgeMidPoint, edgeI)
3551         {
3552             if (edgeMidPoint[edgeI] >= 0)
3553             {
3554                 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
3555
3556                 forAll(eFaces, i)
3557                 {
3558                     affectedFace.set(eFaces[i], 1);
3559                 }
3560             }
3561         }
3562     }
3563
3564
3565     // 1. Faces that get split
3566     // ~~~~~~~~~~~~~~~~~~~~~~~
3567
3568     if (debug)
3569     {
3570         Pout<< "hexRef8::setRefinement : Splitting faces" << endl;
3571     }
3572
3573     forAll(faceMidPoint, faceI)
3574     {
3575         if (faceMidPoint[faceI] >= 0 && affectedFace.get(faceI) == 1)
3576         {
3577             // Face needs to be split and hasn't yet been done in some way
3578             // (affectedFace - is impossible since this is first change but
3579             //  just for completeness)
3580
3581             const face& f = mesh_.faces()[faceI];
3582
3583             // Has original faceI been used (three faces added, original gets
3584             // modified)
3585             bool modifiedFace = false;
3586
3587             label anchorLevel = faceAnchorLevel[faceI];
3588
3589             face newFace(4);
3590
3591             forAll(f, fp)
3592             {
3593                 label pointI = f[fp];
3594
3595                 if (pointLevel_[pointI] <= anchorLevel)
3596                 {
3597                     // point is anchor. Start collecting face.
3598
3599                     DynamicList<label> faceVerts(4);
3600
3601                     faceVerts.append(pointI);
3602
3603                     // Walk forward to mid point.
3604                     // - if next is +2 midpoint is +1
3605                     // - if next is +1 it is midpoint
3606                     // - if next is +0 there has to be edgeMidPoint
3607
3608                     walkFaceToMid
3609                     (
3610                         edgeMidPoint,
3611                         anchorLevel,
3612                         faceI,
3613                         fp,
3614                         faceVerts
3615                     );
3616
3617                     faceVerts.append(faceMidPoint[faceI]);
3618
3619                     walkFaceFromMid
3620                     (
3621                         edgeMidPoint,
3622                         anchorLevel,
3623                         faceI,
3624                         fp,
3625                         faceVerts
3626                     );
3627
3628                     // Convert dynamiclist to face.
3629                     newFace.transfer(faceVerts.shrink());
3630                     faceVerts.clear();
3631
3632                     //Pout<< "Split face:" << faceI << " verts:" << f
3633                     //    << " into quad:" << newFace << endl;
3634
3635                     // Get new owner/neighbour
3636                     label own, nei;
3637                     getFaceNeighbours
3638                     (
3639                         cellAnchorPoints,
3640                         cellAddedCells,
3641                         faceI,
3642                         pointI,          // Anchor point
3643
3644                         own,
3645                         nei
3646                     );
3647
3648
3649                     if (debug)
3650                     {
3651                         if (mesh_.isInternalFace(faceI))
3652                         {
3653                             label oldOwn = mesh_.faceOwner()[faceI];
3654                             label oldNei = mesh_.faceNeighbour()[faceI];
3655
3656                             checkInternalOrientation
3657                             (
3658                                 meshMod,
3659                                 oldOwn,
3660                                 faceI,
3661                                 mesh_.cellCentres()[oldOwn],
3662                                 mesh_.cellCentres()[oldNei],
3663                                 newFace
3664                             );
3665                         }
3666                         else
3667                         {
3668                             label oldOwn = mesh_.faceOwner()[faceI];
3669
3670                             checkBoundaryOrientation
3671                             (
3672                                 meshMod,
3673                                 oldOwn,
3674                                 faceI,
3675                                 mesh_.cellCentres()[oldOwn],
3676                                 mesh_.faceCentres()[faceI],
3677                                 newFace
3678                             );
3679                         }
3680                     }
3681
3682
3683                     if (!modifiedFace)
3684                     {
3685                         modifiedFace = true;
3686
3687                         modFace(meshMod, faceI, newFace, own, nei);
3688                     }
3689                     else
3690                     {
3691                         addFace(meshMod, faceI, newFace, own, nei);
3692                     }
3693                 }
3694             }
3695
3696             // Mark face as having been handled
3697             affectedFace.set(faceI, 0);
3698         }
3699     }
3700
3701
3702     // 2. faces that do not get split but use edges that get split
3703     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3704
3705     if (debug)
3706     {
3707         Pout<< "hexRef8::setRefinement :"
3708             << " Adding edge splits to unsplit faces"
3709             << endl;
3710     }
3711
3712     forAll(edgeMidPoint, edgeI)
3713     {
3714         if (edgeMidPoint[edgeI] >= 0)
3715         {
3716             // Split edge. Check that face not already handled above.
3717
3718             const labelList& eFaces = mesh_.edgeFaces()[edgeI];
3719
3720             forAll(eFaces, i)
3721             {
3722                 label faceI = eFaces[i];
3723
3724                 if (faceMidPoint[faceI] < 0 && affectedFace.get(faceI) == 1)
3725                 {
3726                     // Unsplit face. Add edge splits to face.
3727
3728                     const face& f = mesh_.faces()[faceI];
3729                     const labelList& fEdges = mesh_.faceEdges()[faceI];
3730
3731                     DynamicList<label> newFaceVerts(f.size());
3732
3733                     forAll(f, fp)
3734                     {
3735                         newFaceVerts.append(f[fp]);
3736
3737                         label edgeI = fEdges[fp];
3738
3739                         if (edgeMidPoint[edgeI] >= 0)
3740                         {
3741                             newFaceVerts.append(edgeMidPoint[edgeI]);
3742                         }
3743                     }
3744
3745                     face newFace;
3746                     newFace.transfer(newFaceVerts.shrink());
3747
3748
3749                     // The point with the lowest level should be an anchor
3750                     // point of the neighbouring cells.
3751                     label anchorFp = findMinLevel(f);
3752
3753                     label own, nei;
3754                     getFaceNeighbours
3755                     (
3756                         cellAnchorPoints,
3757                         cellAddedCells,
3758                         faceI,
3759                         f[anchorFp],          // Anchor point
3760
3761                         own,
3762                         nei
3763                     );
3764
3765
3766                     if (debug)
3767                     {
3768                         if (mesh_.isInternalFace(faceI))
3769                         {
3770                             label oldOwn = mesh_.faceOwner()[faceI];
3771                             label oldNei = mesh_.faceNeighbour()[faceI];
3772
3773                             checkInternalOrientation
3774                             (
3775                                 meshMod,
3776                                 oldOwn,
3777                                 faceI,
3778                                 mesh_.cellCentres()[oldOwn],
3779                                 mesh_.cellCentres()[oldNei],
3780                                 newFace
3781                             );
3782                         }
3783                         else
3784                         {
3785                             label oldOwn = mesh_.faceOwner()[faceI];
3786
3787                             checkBoundaryOrientation
3788                             (
3789                                 meshMod,
3790                                 oldOwn,
3791                                 faceI,
3792                                 mesh_.cellCentres()[oldOwn],
3793                                 mesh_.faceCentres()[faceI],
3794                                 newFace
3795                             );
3796                         }
3797                     }
3798
3799                     modFace(meshMod, faceI, newFace, own, nei);
3800
3801                     // Mark face as having been handled
3802                     affectedFace.set(faceI, 0);
3803                 }
3804             }
3805         }
3806     }
3807
3808
3809     // 3. faces that do not get split but whose owner/neighbour change
3810     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3811
3812     if (debug)
3813     {
3814         Pout<< "hexRef8::setRefinement :"
3815             << " Changing owner/neighbour for otherwise unaffected faces"
3816             << endl;
3817     }
3818
3819     forAll(affectedFace, faceI)
3820     {
3821         if (affectedFace.get(faceI) == 1)
3822         {
3823             const face& f = mesh_.faces()[faceI];
3824
3825             // The point with the lowest level should be an anchor
3826             // point of the neighbouring cells.
3827             label anchorFp = findMinLevel(f);
3828
3829             label own, nei;
3830             getFaceNeighbours
3831             (
3832                 cellAnchorPoints,
3833                 cellAddedCells,
3834                 faceI,
3835                 f[anchorFp],          // Anchor point
3836
3837                 own,
3838                 nei
3839             );
3840
3841             modFace(meshMod, faceI, f, own, nei);
3842
3843             // Mark face as having been handled
3844             affectedFace.set(faceI, 0);
3845         }
3846     }
3847
3848
3849     // 4. new internal faces inside split cells.
3850     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3851
3852
3853     // This is the hard one. We have to find the splitting points between
3854     // the anchor points. But the edges between the anchor points might have
3855     // been split (into two,three or four edges).
3856
3857     if (debug)
3858     {
3859         Pout<< "hexRef8::setRefinement :"
3860             << " Create new internal faces for split cells"
3861             << endl;
3862     }
3863
3864     forAll(cellMidPoint, cellI)
3865     {
3866         if (cellMidPoint[cellI] >= 0)
3867         {
3868             createInternalFaces
3869             (
3870                 cellAnchorPoints,
3871                 cellAddedCells,
3872                 cellMidPoint,
3873                 faceMidPoint,
3874                 faceAnchorLevel,
3875                 edgeMidPoint,
3876                 cellI,
3877                 meshMod
3878             );
3879         }
3880     }
3881
3882     // Extend pointLevels and cellLevels for the new cells. Could also be done
3883     // in updateMesh but saves passing cellAddedCells out of this routine.
3884
3885     // Check
3886     if (debug)
3887     {
3888         label minPointI = labelMax;
3889         label maxPointI = labelMin;
3890
3891         forAll(cellMidPoint, cellI)
3892         {
3893             if (cellMidPoint[cellI] >= 0)
3894             {
3895                 minPointI = min(minPointI, cellMidPoint[cellI]);
3896                 maxPointI = max(maxPointI, cellMidPoint[cellI]);
3897             }
3898         }
3899         forAll(faceMidPoint, faceI)
3900         {
3901             if (faceMidPoint[faceI] >= 0)
3902             {
3903                 minPointI = min(minPointI, faceMidPoint[faceI]);
3904                 maxPointI = max(maxPointI, faceMidPoint[faceI]);
3905             }
3906         }
3907         forAll(edgeMidPoint, edgeI)
3908         {
3909             if (edgeMidPoint[edgeI] >= 0)
3910             {
3911                 minPointI = min(minPointI, edgeMidPoint[edgeI]);
3912                 maxPointI = max(maxPointI, edgeMidPoint[edgeI]);
3913             }
3914         }
3915
3916         if (minPointI != labelMax && minPointI != mesh_.nPoints())
3917         {
3918             FatalErrorIn("hexRef8::setRefinement(..)")
3919                 << "Added point labels not consecutive to existing mesh points."
3920                 << nl
3921                 << "mesh_.nPoints():" << mesh_.nPoints()
3922                 << " minPointI:" << minPointI
3923                 << " maxPointI:" << maxPointI
3924                 << abort(FatalError);
3925         }
3926     }
3927
3928     pointLevel_.transfer(newPointLevel.shrink());
3929     newPointLevel.clear();
3930     cellLevel_.transfer(newCellLevel.shrink());
3931     newCellLevel.clear();
3932
3933     // Mark files as changed
3934     setInstance(mesh_.facesInstance());
3935
3936
3937     // Update the live split cells tree.
3938     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3939
3940     // New unrefinement structure
3941     if (history_.active())
3942     {
3943         if (debug)
3944         {
3945             Pout<< "hexRef8::setRefinement :"
3946                 << " Updating refinement history to " << cellLevel_.size()
3947                 << " cells" << endl;
3948         }
3949
3950         // Extend refinement history for new cells
3951         history_.resize(cellLevel_.size());
3952
3953         forAll(cellAddedCells, cellI)
3954         {
3955             const labelList& addedCells = cellAddedCells[cellI];
3956
3957             if (addedCells.size() > 0)
3958             {
3959                 // Cell was split.
3960                 history_.storeSplit(cellI, addedCells);
3961             }
3962         }
3963     }
3964
3965     // Compact cellAddedCells.
3966
3967     labelListList refinedCells(cellLabels.size());
3968
3969     forAll(cellLabels, i)
3970     {
3971         label cellI = cellLabels[i];
3972
3973         refinedCells[i].transfer(cellAddedCells[cellI]);
3974     }
3975
3976     return refinedCells;
3977 }
3978
3979
3980 void Foam::hexRef8::storeData
3981 (
3982     const labelList& pointsToStore,
3983     const labelList& facesToStore,
3984     const labelList& cellsToStore
3985 )
3986 {
3987     savedPointLevel_.resize(2*pointsToStore.size());
3988     forAll(pointsToStore, i)
3989     {
3990         label pointI = pointsToStore[i];
3991         savedPointLevel_.insert(pointI, pointLevel_[pointI]);
3992     }
3993
3994     savedCellLevel_.resize(2*cellsToStore.size());
3995     forAll(cellsToStore, i)
3996     {
3997         label cellI = cellsToStore[i];
3998         savedCellLevel_.insert(cellI, cellLevel_[cellI]);
3999     }
4000 }
4001
4002
4003 // Gets called after the mesh change. setRefinement will already have made
4004 // sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4005 // only need to account for reordering.
4006 void Foam::hexRef8::updateMesh(const mapPolyMesh& map)
4007 {
4008     Map<label> dummyMap(0);
4009
4010     updateMesh(map, dummyMap, dummyMap, dummyMap);
4011 }
4012
4013
4014 // Gets called after the mesh change. setRefinement will already have made
4015 // sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
4016 // only need to account for reordering.
4017 void Foam::hexRef8::updateMesh
4018 (
4019     const mapPolyMesh& map,
4020     const Map<label>& pointsToRestore,
4021     const Map<label>& facesToRestore,
4022     const Map<label>& cellsToRestore
4023 )
4024 {
4025     // Update celllevel
4026     if (debug)
4027     {
4028         Pout<< "hexRef8::updateMesh :"
4029             << " Updating various lists"
4030             << endl;
4031     }
4032
4033     {
4034         const labelList& reverseCellMap = map.reverseCellMap();
4035
4036         if (debug)
4037         {
4038             Pout<< "hexRef8::updateMesh :"
4039                 << " reverseCellMap:" << map.reverseCellMap().size()
4040                 << " cellMap:" << map.cellMap().size()
4041                 << " nCells:" << mesh_.nCells()
4042                 << " nOldCells:" << map.nOldCells()
4043                 << " cellLevel_:" << cellLevel_.size()
4044                 << " reversePointMap:" << map.reversePointMap().size()
4045                 << " pointMap:" << map.pointMap().size()
4046                 << " nPoints:" << mesh_.nPoints()
4047                 << " nOldPoints:" << map.nOldPoints()
4048                 << " pointLevel_:" << pointLevel_.size()
4049                 << endl;
4050         }
4051
4052         if (reverseCellMap.size() == cellLevel_.size())
4053         {
4054             // Assume it is after hexRef8 that this routine is called.
4055             // Just account for reordering. We cannot use cellMap since
4056             // then cells created from cells would get cellLevel_ of
4057             // cell they were created from.
4058             reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4059         }
4060         else
4061         {
4062             // Map data
4063             const labelList& cellMap = map.cellMap();
4064
4065             labelList newCellLevel(cellMap.size());
4066             forAll(cellMap, newCellI)
4067             {
4068                 label oldCellI = cellMap[newCellI];
4069
4070                 if (oldCellI == -1)
4071                 {
4072                     //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4073                     //    << "Problem : cell " << newCellI
4074                     //    << " at " << mesh_.cellCentres()[newCellI]
4075                     //    << " does not originate from another cell"
4076                     //    << " (i.e. is inflated)." << nl
4077                     //    << "Hence we cannot determine the new cellLevel"
4078                     //    << " for it." << abort(FatalError);
4079                     newCellLevel[newCellI] = -1;
4080                 }
4081                 else
4082                 {
4083                     newCellLevel[newCellI] = cellLevel_[oldCellI];
4084                 }
4085             }
4086             cellLevel_.transfer(newCellLevel);
4087         }
4088
4089         // See if any cells to restore. This will be for some new cells
4090         // the corresponding old cell.
4091         forAllConstIter(Map<label>, cellsToRestore, iter)
4092         {
4093             label newCellI = iter.key();
4094             label storedCellI = iter();
4095
4096             Map<label>::iterator fnd = savedCellLevel_.find(storedCellI);
4097
4098             if (fnd == savedCellLevel_.end())
4099             {
4100                 FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4101                     << "Problem : trying to restore old value for new cell "
4102                     << newCellI << " but cannot find old cell " << storedCellI
4103                     << " in map of stored values " << savedCellLevel_
4104                     << abort(FatalError);
4105             }
4106             cellLevel_[newCellI] = fnd();
4107         }
4108
4109         //if (findIndex(cellLevel_, -1) != -1)
4110         //{
4111         //    WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
4112         //        << "Problem : "
4113         //        << "cellLevel_ contains illegal value -1 after mapping
4114         //        << " at cell " << findIndex(cellLevel_, -1) << endl
4115         //        << "This means that another program has inflated cells"
4116         //        << " (created cells out-of-nothing) and hence we don't know"
4117         //        << " their cell level. Continuing with illegal value."
4118         //        << abort(FatalError);
4119         //}
4120     }
4121
4122
4123     // Update pointlevel
4124     {
4125         const labelList& reversePointMap = map.reversePointMap();
4126
4127         if (reversePointMap.size() == pointLevel_.size())
4128         {
4129             // Assume it is after hexRef8 that this routine is called.
4130             reorder(reversePointMap, mesh_.nPoints(), -1,  pointLevel_);
4131         }
4132         else
4133         {
4134             // Map data
4135             const labelList& pointMap = map.pointMap();
4136
4137             labelList newPointLevel(pointMap.size());
4138
4139             forAll(pointMap, newPointI)
4140             {
4141                 label oldPointI = pointMap[newPointI];
4142
4143                 if (oldPointI == -1)
4144                 {
4145                     //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4146                     //    << "Problem : point " << newPointI
4147                     //    << " at " << mesh_.points()[newPointI]
4148                     //    << " does not originate from another point"
4149                     //    << " (i.e. is inflated)." << nl
4150                     //    << "Hence we cannot determine the new pointLevel"
4151                     //    << " for it." << abort(FatalError);
4152                     newPointLevel[newPointI] = -1;
4153                 }
4154                 else
4155                 {
4156                     newPointLevel[newPointI] = pointLevel_[oldPointI];
4157                 }
4158             }
4159             pointLevel_.transfer(newPointLevel);
4160         }
4161
4162         // See if any points to restore. This will be for some new points
4163         // the corresponding old point (the one from the call to storeData)
4164         forAllConstIter(Map<label>, pointsToRestore, iter)
4165         {
4166             label newPointI = iter.key();
4167             label storedPointI = iter();
4168
4169             Map<label>::iterator fnd = savedPointLevel_.find(storedPointI);
4170
4171             if (fnd == savedPointLevel_.end())
4172             {
4173                 FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
4174                     << "Problem : trying to restore old value for new point "
4175                     << newPointI << " but cannot find old point "
4176                     << storedPointI
4177                     << " in map of stored values " << savedPointLevel_
4178                     << abort(FatalError);
4179             }
4180             pointLevel_[newPointI] = fnd();
4181         }
4182
4183         //if (findIndex(pointLevel_, -1) != -1)
4184         //{
4185         //    WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
4186         //        << "Problem : "
4187         //        << "pointLevel_ contains illegal value -1 after mapping"
4188         //        << " at point" << findIndex(pointLevel_, -1) << endl
4189         //        << "This means that another program has inflated points"
4190         //        << " (created points out-of-nothing) and hence we don't know"
4191         //        << " their point level. Continuing with illegal value."
4192         //        //<< abort(FatalError);
4193         //}
4194     }
4195
4196     // Update refinement tree
4197     if (history_.active())
4198     {
4199         history_.updateMesh(map);
4200     }
4201
4202     // Mark files as changed
4203     setInstance(mesh_.facesInstance());
4204
4205     // Update face removal engine
4206     faceRemover_.updateMesh(map);
4207 }
4208
4209
4210 // Gets called after mesh subsetting. Maps are from new back to old.
4211 void Foam::hexRef8::subset
4212 (
4213     const labelList& pointMap,
4214     const labelList& faceMap,
4215     const labelList& cellMap
4216 )
4217 {
4218     // Update celllevel
4219     if (debug)
4220     {
4221         Pout<< "hexRef8::subset :"
4222             << " Updating various lists"
4223             << endl;
4224     }
4225
4226     if (history_.active())
4227     {
4228         WarningIn
4229         (
4230             "hexRef8::subset(const labelList&, const labelList&"
4231             ", const labelList&)"
4232         )   << "Subsetting will not work in combination with unrefinement."
4233             << nl
4234             << "Proceed at your own risk." << endl;
4235     }
4236
4237
4238     // Update celllevel
4239     {
4240         labelList newCellLevel(cellMap.size());
4241
4242         forAll(cellMap, newCellI)
4243         {
4244             newCellLevel[newCellI] = cellLevel_[cellMap[newCellI]];
4245         }
4246
4247         cellLevel_.transfer(newCellLevel);
4248
4249         if (findIndex(cellLevel_, -1) != -1)
4250         {
4251             FatalErrorIn("hexRef8::subset(..)")
4252                 << "Problem : "
4253                 << "cellLevel_ contains illegal value -1 after mapping:"
4254                 << cellLevel_
4255                 << abort(FatalError);
4256         }
4257     }
4258
4259     // Update pointlevel
4260     {
4261         labelList newPointLevel(pointMap.size());
4262
4263         forAll(pointMap, newPointI)
4264         {
4265             newPointLevel[newPointI] = pointLevel_[pointMap[newPointI]];
4266         }
4267
4268         pointLevel_.transfer(newPointLevel);
4269
4270         if (findIndex(pointLevel_, -1) != -1)
4271         {
4272             FatalErrorIn("hexRef8::subset(..)")
4273                 << "Problem : "
4274                 << "pointLevel_ contains illegal value -1 after mapping:"
4275                 << pointLevel_
4276                 << abort(FatalError);
4277         }
4278     }
4279
4280     // Update refinement tree
4281     if (history_.active())
4282     {
4283         history_.subset(pointMap, faceMap, cellMap);
4284     }
4285
4286     // Mark files as changed
4287     setInstance(mesh_.facesInstance());
4288
4289     // Nothing needs doing to faceRemover.
4290     //faceRemover_.subset(pointMap, faceMap, cellMap);
4291 }
4292
4293
4294 // Gets called after the mesh distribution
4295 void Foam::hexRef8::distribute(const mapDistributePolyMesh& map)
4296 {
4297     if (debug)
4298     {
4299         Pout<< "hexRef8::distribute :"
4300             << " Updating various lists"
4301             << endl;
4302     }
4303
4304     // Update celllevel
4305     map.distributeCellData(cellLevel_);
4306     // Update pointlevel
4307     map.distributePointData(pointLevel_);
4308
4309     // Update refinement tree
4310     if (history_.active())
4311     {
4312         history_.distribute(map);
4313     }
4314
4315     // Update face removal engine
4316     faceRemover_.distribute(map);
4317 }
4318
4319
4320 void Foam::hexRef8::checkMesh() const
4321 {
4322     const boundBox& meshBb = mesh_.globalData().bb();
4323     const scalar smallDim = 1E-6*mag(meshBb.max() - meshBb.min());
4324
4325     if (debug)
4326     {
4327         Pout<< "hexRef8::checkMesh : Using matching tolerance smallDim:"
4328             << smallDim << endl;
4329     }
4330
4331     // Check owner on coupled faces.
4332     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4333
4334     // There should be only one coupled face between two cells. Why? Since
4335     // otherwise mesh redistribution might cause multiple faces between two
4336     // cells
4337     {
4338         labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
4339         forAll(nei, i)
4340         {
4341             nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4342         }
4343
4344         // Replace data on coupled patches with their neighbour ones.
4345         syncTools::swapBoundaryFaceList(mesh_, nei, false);
4346
4347         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
4348
4349         forAll(patches, patchI)
4350         {
4351             const polyPatch& pp = patches[patchI];
4352
4353             if (pp.coupled())
4354             {
4355                 // Check how many faces between owner and neighbour. Should
4356                 // be only one.
4357                 HashTable<label, labelPair, labelPair::Hash<> >
4358                     cellToFace(2*pp.size());
4359
4360                 label faceI = pp.start();
4361
4362                 forAll(pp, i)
4363                 {
4364                     label own = mesh_.faceOwner()[faceI];
4365                     label bFaceI = faceI-mesh_.nInternalFaces();
4366
4367                     if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
4368                     {
4369                         FatalErrorIn("hexRef8::checkMesh()")
4370                             << "Faces do not seem to be correct across coupled"
4371                             << " boundaries" << endl
4372                             << "Coupled face " << faceI
4373                             << " between owner " << own
4374                             << " on patch " << pp.name()
4375                             << " and coupled neighbour " << nei[bFaceI]
4376                             << " has two faces connected to it:"
4377                             << faceI << " and "
4378                             << cellToFace[labelPair(own, nei[bFaceI])]
4379                             << abort(FatalError);
4380                     }
4381
4382                     faceI++;
4383                 }
4384             }
4385         }
4386     }
4387
4388     // Check face areas.
4389     // ~~~~~~~~~~~~~~~~~
4390
4391     {
4392         scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
4393         forAll(neiFaceAreas, i)
4394         {
4395             neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4396         }
4397
4398         // Replace data on coupled patches with their neighbour ones.
4399         syncTools::swapBoundaryFaceList(mesh_, neiFaceAreas, false);
4400
4401         forAll(neiFaceAreas, i)
4402         {
4403             label faceI = i+mesh_.nInternalFaces();
4404
4405             const scalar magArea = mag(mesh_.faceAreas()[faceI]);
4406
4407             if (mag(magArea - neiFaceAreas[i]) > smallDim)
4408             {
4409                 const face& f = mesh_.faces()[faceI];
4410                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4411
4412                 FatalErrorIn("hexRef8::checkMesh()")
4413                     << "Faces do not seem to be correct across coupled"
4414                     << " boundaries" << endl
4415                     << "Coupled face " << faceI
4416                     << " on patch " << patchI
4417                     << " " << mesh_.boundaryMesh()[patchI].name()
4418                     << " coords:" << IndirectList<point>(mesh_.points(), f)()
4419                     << " has face area:" << magArea
4420                     << " (coupled) neighbour face area differs:"
4421                     << neiFaceAreas[i]
4422                     << " to within tolerance " << smallDim
4423                     << abort(FatalError);
4424             }
4425         }
4426     }
4427
4428
4429     // Check number of points on faces.
4430     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4431     {
4432         labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
4433
4434         forAll(nVerts, i)
4435         {
4436             nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4437         }
4438
4439         // Replace data on coupled patches with their neighbour ones.
4440         syncTools::swapBoundaryFaceList(mesh_, nVerts, false);
4441
4442         forAll(nVerts, i)
4443         {
4444             label faceI = i+mesh_.nInternalFaces();
4445
4446             const face& f = mesh_.faces()[faceI];
4447
4448             if (f.size() != nVerts[i])
4449             {
4450                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4451
4452                 FatalErrorIn("hexRef8::checkMesh()")
4453                     << "Faces do not seem to be correct across coupled"
4454                     << " boundaries" << endl
4455                     << "Coupled face " << faceI
4456                     << " on patch " << patchI
4457                     << " " << mesh_.boundaryMesh()[patchI].name()
4458                     << " coords:" << IndirectList<point>(mesh_.points(), f)()
4459                     << " has size:" << f.size()
4460                     << " (coupled) neighbour face has size:"
4461                     << nVerts[i]
4462                     << abort(FatalError);
4463             }
4464         }
4465     }
4466
4467
4468     // Check points of face
4469     // ~~~~~~~~~~~~~~~~~~~~
4470     {
4471         // Anchor points.
4472         pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
4473
4474         forAll(anchorPoints, i)
4475         {
4476             label faceI = i+mesh_.nInternalFaces();
4477             const point& fc = mesh_.faceCentres()[faceI];
4478             const face& f = mesh_.faces()[faceI];
4479             const vector anchorVec(mesh_.points()[f[0]] - fc);
4480
4481             anchorPoints[i] = anchorVec;
4482         }
4483
4484         // Replace data on coupled patches with their neighbour ones. Apply
4485         // rotation transformation (but not separation since is relative vector
4486         // to point on same face.
4487         syncTools::swapBoundaryFaceList(mesh_, anchorPoints, false);
4488
4489         forAll(anchorPoints, i)
4490         {
4491             label faceI = i+mesh_.nInternalFaces();
4492             const point& fc = mesh_.faceCentres()[faceI];
4493             const face& f = mesh_.faces()[faceI];
4494             const vector anchorVec(mesh_.points()[f[0]] - fc);
4495
4496             if (mag(anchorVec - anchorPoints[i]) > smallDim)
4497             {
4498                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4499
4500                 FatalErrorIn("hexRef8::checkMesh()")
4501                     << "Faces do not seem to be correct across coupled"
4502                     << " boundaries" << endl
4503                     << "Coupled face " << faceI
4504                     << " on patch " << patchI
4505                     << " " << mesh_.boundaryMesh()[patchI].name()
4506                     << " coords:" << IndirectList<point>(mesh_.points(), f)()
4507                     << " has anchor vector:" << anchorVec
4508                     << " (coupled) neighbour face anchor vector differs:"
4509                     << anchorPoints[i]
4510                     << " to within tolerance " << smallDim
4511                     << abort(FatalError);
4512             }
4513         }
4514     }
4515
4516     if (debug)
4517     {
4518         Pout<< "hexRef8::checkMesh : Returning" << endl;
4519     }
4520 }
4521
4522
4523 void Foam::hexRef8::checkRefinementLevels
4524 (
4525     const label maxPointDiff,
4526     const labelList& pointsToCheck
4527 ) const
4528 {
4529     if (debug)
4530     {
4531         Pout<< "hexRef8::checkRefinementLevels :"
4532             << " Checking 2:1 refinement level" << endl;
4533     }
4534
4535     if
4536     (
4537         cellLevel_.size() != mesh_.nCells()
4538      || pointLevel_.size() != mesh_.nPoints()
4539     )
4540     {
4541         FatalErrorIn("hexRef8::checkRefinementLevels(const label)")
4542             << "cellLevel size should be number of cells"
4543             << " and pointLevel size should be number of points."<< nl
4544             << "cellLevel:" << cellLevel_.size()
4545             << " mesh.nCells():" << mesh_.nCells() << nl
4546             << "pointLevel:" << pointLevel_.size()
4547             << " mesh.nPoints():" << mesh_.nPoints()
4548             << abort(FatalError);
4549     }
4550
4551
4552     // Check 2:1 consistency.
4553     // ~~~~~~~~~~~~~~~~~~~~~~
4554
4555     {
4556         // Internal faces.
4557         for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
4558         {
4559             label own = mesh_.faceOwner()[faceI];
4560             label nei = mesh_.faceNeighbour()[faceI];
4561
4562             if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4563             {
4564                 FatalErrorIn
4565                 (
4566                     "hexRef8::checkRefinementLevels(const label)"
4567                 )   << "Celllevel does not satisfy 2:1 constraint." << nl
4568                     << "On face " << faceI << " owner cell " << own
4569                     << " has refinement " << cellLevel_[own]
4570                     << " neighbour cell " << nei << " has refinement "
4571                     << cellLevel_[nei]
4572                     << abort(FatalError);
4573             }
4574         }
4575
4576         // Coupled faces. Get neighbouring value
4577         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
4578
4579         forAll(neiLevel, i)
4580         {
4581             label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4582
4583             neiLevel[i] = cellLevel_[own];
4584         }
4585
4586         // No separation
4587         syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
4588
4589         forAll(neiLevel, i)
4590         {
4591             label faceI = i+mesh_.nInternalFaces();
4592
4593             label own = mesh_.faceOwner()[faceI];
4594
4595             if (mag(cellLevel_[own] - neiLevel[i]) > 1)
4596             {
4597                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4598
4599                 FatalErrorIn
4600                 (
4601                     "hexRef8::checkRefinementLevels(const label)"
4602                 )   << "Celllevel does not satisfy 2:1 constraint."
4603                     << " On coupled face " << faceI
4604                     << " on patch " << patchI << " "
4605                     << mesh_.boundaryMesh()[patchI].name()
4606                     << " owner cell " << own << " has refinement "
4607                     << cellLevel_[own]
4608                     << " (coupled) neighbour cell has refinement "
4609                     << neiLevel[i]
4610                     << abort(FatalError);
4611             }
4612         }
4613     }
4614
4615
4616     // Check pointLevel is synchronized
4617     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4618     {
4619         labelList syncPointLevel(pointLevel_);
4620
4621         // Get min level
4622         syncTools::syncPointList
4623         (
4624             mesh_,
4625             syncPointLevel,
4626             minEqOp<label>(),
4627             labelMax,
4628             false               // no separation
4629         );
4630
4631
4632         forAll(syncPointLevel, pointI)
4633         {
4634             if (pointLevel_[pointI] != syncPointLevel[pointI])
4635             {
4636                 FatalErrorIn
4637                 (
4638                     "hexRef8::checkRefinementLevels(const label)"
4639                 )   << "PointLevel is not consistent across coupled patches."
4640                     << endl
4641                     << "point:" << pointI << " coord:" << mesh_.points()[pointI]
4642                     << " has level " << pointLevel_[pointI]
4643                     << " whereas the coupled point has level "
4644                     << syncPointLevel[pointI]
4645                     << abort(FatalError);
4646             }
4647         }
4648     }
4649
4650
4651     // Check 2:1 across points (instead of faces)
4652     if (maxPointDiff != -1)
4653     {
4654         const labelListList& pointCells = mesh_.pointCells();
4655
4656         // Determine per point the max cell level.
4657         labelList maxPointLevel(mesh_.nPoints(), 0);
4658
4659         forAll(pointCells, pointI)
4660         {
4661             const labelList& pCells = pointCells[pointI];
4662
4663             label& pLevel = maxPointLevel[pointI];
4664
4665             forAll(pCells, i)
4666             {
4667                 pLevel = max(pLevel, cellLevel_[pCells[i]]);
4668             }
4669         }
4670
4671         // Sync maxPointLevel to neighbour
4672         syncTools::syncPointList
4673         (
4674             mesh_,
4675             maxPointLevel,
4676             maxEqOp<label>(),
4677             labelMin,           // null value
4678             false               // no separation
4679         );
4680
4681         // Check 2:1 across boundary points
4682         forAll(pointsToCheck, i)
4683         {
4684             label pointI = pointsToCheck[i];
4685
4686             const labelList& pCells = pointCells[pointI];
4687
4688             forAll(pCells, i)
4689             {
4690                 label cellI = pCells[i];
4691
4692                 if
4693                 (
4694                     mag(cellLevel_[cellI]-maxPointLevel[pointI])
4695                   > maxPointDiff
4696                 )
4697                 {
4698                     FatalErrorIn
4699                     (
4700                         "hexRef8::checkRefinementLevels(const label)"
4701                     )   << "Too big a difference between"
4702                         << " point-connected cells." << nl
4703                         << "cell:" << cellI
4704                         << " cellLevel:" << cellLevel_[cellI]
4705                         << " uses point:" << pointI
4706                         << " coord:" << mesh_.points()[pointI]
4707                         << " which is also used by a cell with level:"
4708                         << maxPointLevel[pointI]
4709                         << abort(FatalError);
4710                 }
4711             }
4712         }
4713     }
4714 }
4715
4716
4717
4718 //
4719 // Unrefinement
4720 // ~~~~~~~~~~~~
4721 //
4722
4723
4724 Foam::labelList Foam::hexRef8::getSplitPoints() const
4725 {
4726     if (debug)
4727     {
4728         checkRefinementLevels(-1, labelList(0));
4729     }
4730
4731     if (debug)
4732     {
4733         Pout<< "hexRef8::getSplitPoints :"
4734             << " Calculating unrefineable points" << endl;
4735     }
4736
4737
4738     if (!history_.active())
4739     {
4740         FatalErrorIn("hexRef8::getSplitPoints()")
4741             << "Only call if constructed with history capability"
4742             << abort(FatalError);
4743     }
4744
4745     // Master cell
4746     // -1 undetermined
4747     // -2 certainly not split point
4748     // >= label of master cell
4749     labelList splitMaster(mesh_.nPoints(), -1);
4750     labelList splitMasterLevel(mesh_.nPoints(), 0);
4751
4752     // Unmark all with not 8 cells
4753     const labelListList& pointCells = mesh_.pointCells();
4754
4755     forAll(pointCells, pointI)
4756     {
4757         const labelList& pCells = pointCells[pointI];
4758
4759         if (pCells.size() != 8)
4760         {
4761             splitMaster[pointI] = -2;
4762         }
4763     }
4764
4765     // Unmark all with different master cells
4766     const labelList& visibleCells = history_.visibleCells();
4767
4768     forAll(visibleCells, cellI)
4769     {
4770         //const labelList& cPoints = mesh_.cellPoints()[cellI];
4771         const labelList cPoints(cellPoints(cellI));
4772
4773         if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
4774         {
4775             label parentIndex = history_.parentIndex(cellI);
4776
4777             // Check same master.
4778             forAll(cPoints, i)
4779             {
4780                 label pointI = cPoints[i];
4781
4782                 label masterCellI = splitMaster[pointI];
4783
4784                 if (masterCellI == -1)
4785                 {
4786                     // First time visit of point. Store parent cell and
4787                     // level of the parent cell (with respect to cellI). This
4788                     // is additional guarantee that we're referring to the
4789                     // same master at the same refinement level.
4790
4791                     splitMaster[pointI] = parentIndex;
4792                     splitMasterLevel[pointI] = cellLevel_[cellI] - 1;
4793                 }
4794                 else if (masterCellI == -2)
4795                 {
4796                     // Already decided that point is not splitPoint
4797                 }
4798                 else if
4799                 (
4800                     (masterCellI != parentIndex)
4801                  || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1)
4802                 )
4803                 {
4804                     // Different masters so point is on two refinement
4805                     // patterns
4806                     splitMaster[pointI] = -2;
4807                 }
4808             }
4809         }
4810         else
4811         {
4812             // Either not visible or is unrefined cell
4813             forAll(cPoints, i)
4814             {
4815                 label pointI = cPoints[i];
4816
4817                 splitMaster[pointI] = -2;
4818             }
4819         }
4820     }
4821
4822     // Unmark boundary faces
4823     for
4824     (
4825         label faceI = mesh_.nInternalFaces();
4826         faceI < mesh_.nFaces();
4827         faceI++
4828     )
4829     {
4830         const face& f = mesh_.faces()[faceI];
4831
4832         forAll(f, fp)
4833         {
4834             splitMaster[f[fp]] = -2;
4835         }
4836     }
4837
4838
4839     // Collect into labelList
4840
4841     label nSplitPoints = 0;
4842
4843     forAll(splitMaster, pointI)
4844     {
4845         if (splitMaster[pointI] >= 0)
4846         {
4847             nSplitPoints++;
4848         }
4849     }
4850
4851     labelList splitPoints(nSplitPoints);
4852     nSplitPoints = 0;
4853
4854     forAll(splitMaster, pointI)
4855     {
4856         if (splitMaster[pointI] >= 0)
4857         {
4858             splitPoints[nSplitPoints++] = pointI;
4859         }
4860     }
4861
4862     return splitPoints;
4863 }
4864
4865
4866 //void Foam::hexRef8::markIndex
4867 //(
4868 //    const label maxLevel,
4869 //    const label level,
4870 //    const label index,
4871 //    const label markValue,
4872 //    labelList& indexValues
4873 //) const
4874 //{
4875 //    if (level < maxLevel && indexValues[index] == -1)
4876 //    {
4877 //        // Mark
4878 //        indexValues[index] = markValue;
4879 //
4880 //        // Mark parent
4881 //        const splitCell8& split = history_.splitCells()[index];
4882 //
4883 //        if (split.parent_ >= 0)
4884 //        {
4885 //            markIndex
4886 //            (
4887 //              maxLevel, level+1, split.parent_, markValue, indexValues);
4888 //            )
4889 //        }
4890 //    }
4891 //}
4892 //
4893 //
4894 //// Get all cells which (down to level) originate from the same cell.
4895 //// level=0 returns cell only, level=1 returns the 8 cells this cell
4896 //// originates from, level=2 returns 64 cells etc.
4897 //// If the cell does not originate from refinement returns just itself.
4898 //void Foam::hexRef8::markCellClusters
4899 //(
4900 //    const label maxLevel,
4901 //    labelList& cluster
4902 //) const
4903 //{
4904 //    cluster.setSize(mesh_.nCells());
4905 //    cluster = -1;
4906 //
4907 //    const DynamicList<splitCell8>& splitCells = history_.splitCells();
4908 //
4909 //    // Mark all splitCells down to level maxLevel with a cell originating from
4910 //    // it.
4911 //
4912 //    labelList indexLevel(splitCells.size(), -1);
4913 //
4914 //    forAll(visibleCells, cellI)
4915 //    {
4916 //        label index = visibleCells[cellI];
4917 //
4918 //        if (index >= 0)
4919 //        {
4920 //            markIndex(maxLevel, 0, index, cellI, indexLevel);
4921 //        }
4922 //    }
4923 //
4924 //    // Mark cells with splitCell
4925 //}
4926
4927
4928 Foam::labelList Foam::hexRef8::consistentUnrefinement
4929 (
4930     const labelList& pointsToUnrefine,
4931     const bool maxSet
4932 ) const
4933 {
4934     if (debug)
4935     {
4936         Pout<< "hexRef8::consistentUnrefinement :"
4937             << " Determining 2:1 consistent unrefinement" << endl;
4938     }
4939
4940     if (maxSet)
4941     {
4942         FatalErrorIn
4943         (
4944             "hexRef8::consistentUnrefinement(const labelList&, const bool"
4945         )   << "maxSet not implemented yet."
4946             << abort(FatalError);
4947     }
4948
4949     // Loop, modifying pointsToUnrefine, until no more changes to due to 2:1
4950     // conflicts.
4951     // maxSet = false : unselect points to refine
4952     // maxSet = true: select points to refine
4953
4954     // Maintain boolList for pointsToUnrefine and cellsToUnrefine
4955     PackedList<1> unrefinePoint(mesh_.nPoints(), 0);
4956
4957     forAll(pointsToUnrefine, i)
4958     {
4959         label pointI = pointsToUnrefine[i];
4960
4961         unrefinePoint.set(pointI, 1);
4962     }
4963
4964
4965     while (true)
4966     {
4967         // Construct cells to unrefine
4968         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
4969
4970         PackedList<1> unrefineCell(mesh_.nCells(), 0);
4971
4972         forAll(unrefinePoint, pointI)
4973         {
4974             if (unrefinePoint.get(pointI) == 1)
4975             {
4976                 const labelList& pCells = mesh_.pointCells()[pointI];
4977
4978                 forAll(pCells, j)
4979                 {
4980                     unrefineCell.set(pCells[j], 1);
4981                 }
4982             }
4983         }
4984
4985
4986         label nChanged = 0;
4987
4988
4989         // Check 2:1 consistency taking refinement into account
4990         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4991
4992         // Internal faces.
4993         for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
4994         {
4995             label own = mesh_.faceOwner()[faceI];
4996             label ownLevel = cellLevel_[own] - unrefineCell.get(own);
4997
4998             label nei = mesh_.faceNeighbour()[faceI];
4999             label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5000
5001             if (ownLevel < (neiLevel-1))
5002             {
5003                 // Since was 2:1 this can only occur if own is marked for
5004                 // unrefinement.
5005
5006                 if (maxSet)
5007                 {
5008                     unrefineCell.set(nei, 1);
5009                 }
5010                 else
5011                 {
5012                     if (unrefineCell.get(own) == 0)
5013                     {
5014                         FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5015                             << "problem" << abort(FatalError);
5016                     }
5017
5018                     unrefineCell.set(own, 0);
5019                 }
5020                 nChanged++;
5021             }
5022             else if (neiLevel < (ownLevel-1))
5023             {
5024                 if (maxSet)
5025                 {
5026                     unrefineCell.set(own, 1);
5027                 }
5028                 else
5029                 {
5030                     if (unrefineCell.get(nei) == 0)
5031                     {
5032                         FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5033                             << "problem" << abort(FatalError);
5034                     }
5035
5036                     unrefineCell.set(nei, 0);
5037                 }
5038                 nChanged++;
5039             }
5040         }
5041
5042
5043         // Coupled faces. Swap owner level to get neighbouring cell level.
5044         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
5045
5046         forAll(neiLevel, i)
5047         {
5048             label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5049
5050             neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5051         }
5052
5053         // Swap to neighbour
5054         syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
5055
5056         forAll(neiLevel, i)
5057         {
5058             label faceI = i+mesh_.nInternalFaces();
5059             label own = mesh_.faceOwner()[faceI];
5060             label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5061
5062             if (ownLevel < (neiLevel[i]-1))
5063             {
5064                 if (!maxSet)
5065                 {
5066                     if (unrefineCell.get(own) == 0)
5067                     {
5068                         FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5069                             << "problem" << abort(FatalError);
5070                     }
5071
5072                     unrefineCell.set(own, 0);
5073                     nChanged++;
5074                 }
5075             }
5076             else if (neiLevel[i] < (ownLevel-1))
5077             {
5078                 if (maxSet)
5079                 {
5080                     if (unrefineCell.get(own) == 1)
5081                     {
5082                         FatalErrorIn("hexRef8::consistentUnrefinement(..)")
5083                             << "problem" << abort(FatalError);
5084                     }
5085
5086                     unrefineCell.set(own, 1);
5087                     nChanged++;
5088                 }
5089             }
5090         }
5091
5092         reduce(nChanged, sumOp<label>());
5093
5094         if (debug)
5095         {
5096             Pout<< "hexRef8::consistentUnrefinement :"
5097                 << " Changed " << nChanged
5098                 << " refinement levels due to 2:1 conflicts."
5099                 << endl;
5100         }
5101
5102         if (nChanged == 0)
5103         {
5104             break;
5105         }
5106
5107
5108         // Convert cellsToUnrefine back into points to unrefine
5109         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5110
5111         // Knock out any point whose cell neighbour cannot be unrefined.
5112         forAll(unrefinePoint, pointI)
5113         {
5114             if (unrefinePoint.get(pointI) == 1)
5115             {
5116                 const labelList& pCells = mesh_.pointCells()[pointI];
5117
5118                 forAll(pCells, j)
5119                 {
5120                     if (unrefineCell.get(pCells[j]) == 0)
5121                     {
5122                         unrefinePoint.set(pointI, 0);
5123                         break;
5124                     }
5125                 }
5126             }
5127         }
5128     }
5129
5130
5131     // Convert back to labelList.
5132     label nSet = 0;
5133
5134     forAll(unrefinePoint, pointI)
5135     {
5136         if (unrefinePoint.get(pointI) == 1)
5137         {
5138             nSet++;
5139         }
5140     }
5141
5142     labelList newPointsToUnrefine(nSet);
5143     nSet = 0;
5144
5145     forAll(unrefinePoint, pointI)
5146     {
5147         if (unrefinePoint.get(pointI) == 1)
5148         {
5149             newPointsToUnrefine[nSet++] = pointI;
5150         }
5151     }
5152
5153     return newPointsToUnrefine;
5154 }
5155
5156
5157 void Foam::hexRef8::setUnrefinement
5158 (
5159     const labelList& splitPointLabels,
5160     polyTopoChange& meshMod
5161 )
5162 {
5163     if (!history_.active())
5164     {
5165         FatalErrorIn
5166         (
5167             "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5168         )   << "Only call if constructed with history capability"
5169             << abort(FatalError);
5170     }
5171
5172     if (debug)
5173     {
5174         Pout<< "hexRef8::setUnrefinement :"
5175             << " Checking initial mesh just to make sure" << endl;
5176
5177         checkMesh();
5178
5179         forAll(cellLevel_, cellI)
5180         {
5181             if (cellLevel_[cellI] < 0)
5182             {
5183                 FatalErrorIn
5184                 (
5185                     "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5186                 )   << "Illegal cell level " << cellLevel_[cellI]
5187                     << " for cell " << cellI
5188                     << abort(FatalError);
5189             }
5190         }
5191
5192
5193         // Write to sets.
5194         pointSet pSet(mesh_, "splitPoints", splitPointLabels);
5195         pSet.write();
5196
5197         cellSet cSet(mesh_, "splitPointCells", splitPointLabels.size());
5198
5199         forAll(splitPointLabels, i)
5200         {
5201             const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]];
5202
5203             forAll(pCells, j)
5204             {
5205                 cSet.insert(pCells[j]);
5206             }
5207         }
5208         cSet.write();
5209
5210         Pout<< "hexRef8::setRefinement : Dumping " << pSet.size()
5211             << " points and "
5212             << cSet.size() << " cells for unrefinement to" << nl
5213             << "    pointSet " << pSet.objectPath() << nl
5214             << "    cellSet " << cSet.objectPath()
5215             << endl;
5216     }
5217
5218
5219     labelList cellRegion;
5220     labelList cellRegionMaster;
5221     labelList facesToRemove;
5222
5223     {
5224         labelHashSet splitFaces(12*splitPointLabels.size());
5225
5226         forAll(splitPointLabels, i)
5227         {
5228             const labelList& pFaces = mesh_.pointFaces()[splitPointLabels[i]];
5229
5230             forAll(pFaces, j)
5231             {
5232                 splitFaces.insert(pFaces[j]);
5233             }
5234         }
5235
5236         // Check with faceRemover what faces will get removed. Note that this
5237         // can be more (but never less) than splitFaces provided.
5238         faceRemover_.compatibleRemoves
5239         (
5240             splitFaces.toc(),   // pierced faces
5241             cellRegion,         // per cell -1 or region it is merged into
5242             cellRegionMaster,   // per region the master cell
5243             facesToRemove       // new faces to be removed.
5244         );
5245
5246         if (facesToRemove.size() != splitFaces.size())
5247         {
5248             FatalErrorIn
5249             (
5250                 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5251             )   << "Ininitial set of split points to unrefine does not"
5252                 << " seem to be consistent or not mid points of refined cells"
5253                 << abort(FatalError);
5254         }
5255     }
5256
5257     // Redo the region master so it is consistent with our master.
5258     // This will guarantee that the new cell (for which faceRemover uses
5259     // the region master) is already compatible with our refinement structure.
5260
5261     forAll(splitPointLabels, i)
5262     {
5263         label pointI = splitPointLabels[i];
5264
5265         // Get original cell label
5266
5267         const labelList& pCells = mesh_.pointCells()[pointI];
5268
5269         // Check
5270         if (pCells.size() != 8)
5271         {
5272             FatalErrorIn
5273             (
5274                 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5275             )   << "splitPoint " << pointI
5276                 << " should have 8 cells using it. It has " << pCells
5277                 << abort(FatalError);
5278         }
5279
5280
5281         // Check that the lowest numbered pCells is the master of the region
5282         // (should be guaranteed by directRemoveFaces)
5283         //if (debug)
5284         {
5285             label masterCellI = min(pCells);
5286
5287             forAll(pCells, j)
5288             {
5289                 label cellI = pCells[j];
5290
5291                 label region = cellRegion[cellI];
5292
5293                 if (region == -1)
5294                 {
5295                     FatalErrorIn("hexRef8::setUnrefinement(..)")
5296                         << "Ininitial set of split points to unrefine does not"
5297                         << " seem to be consistent or not mid points"
5298                         << " of refined cells" << nl
5299                         << "cell:" << cellI << " on splitPoint " << pointI
5300                         << " has no region to be merged into"
5301                         << abort(FatalError);
5302                 }
5303
5304                 if (masterCellI != cellRegionMaster[region])
5305                 {
5306                     FatalErrorIn("hexRef8::setUnrefinement(..)")
5307                         << "cell:" << cellI << " on splitPoint:" << pointI
5308                         << " in region " << region
5309                         << " has master:" << cellRegionMaster[region]
5310                         << " which is not the lowest numbered cell"
5311                         << " among the pointCells:" << pCells
5312                         << abort(FatalError);
5313                 }
5314             }
5315         }
5316     }
5317
5318     // Insert all commands to combine cells. Never fails so don't have to
5319     // test for success.
5320     faceRemover_.setRefinement
5321     (
5322         facesToRemove,
5323         cellRegion,
5324         cellRegionMaster,
5325         meshMod
5326     );
5327
5328     // Remove the 8 cells that originated from merging around the split point
5329     // and adapt cell levels (not that pointLevels stay the same since points
5330     // either get removed or stay at the same position.
5331     forAll(splitPointLabels, i)
5332     {
5333         label pointI = splitPointLabels[i];
5334
5335         const labelList& pCells = mesh_.pointCells()[pointI];
5336
5337         label masterCellI = min(pCells);
5338
5339         forAll(pCells, j)
5340         {
5341             cellLevel_[pCells[j]]--;
5342         }
5343
5344         history_.combineCells(masterCellI, pCells);
5345     }
5346
5347     // Mark files as changed
5348     setInstance(mesh_.facesInstance());
5349
5350     // history_.updateMesh will take care of truncating.
5351 }
5352
5353
5354 // Write refinement to polyMesh directory.
5355 bool Foam::hexRef8::write() const
5356 {
5357     bool writeOk = cellLevel_.write() && pointLevel_.write();
5358
5359     if (history_.active())
5360     {
5361         writeOk = writeOk && history_.write();
5362     }
5363
5364     return writeOk;
5365 }
5366
5367
5368 // ************************************************************************* //