src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.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 "meshRefinement.H"
  28 #include "trackedParticle.H"
  29 #include "syncTools.H"
  30 #include "Time.H"
  31 #include "refinementSurfaces.H"
  32 #include "shellSurfaces.H"
  33 #include "faceSet.H"
  34 #include "decompositionMethod.H"
  35 #include "fvMeshDistribute.H"
  36 #include "polyTopoChange.H"
  37 #include "mapDistributePolyMesh.H"
  38 #include "featureEdgeMesh.H"
  39 #include "Cloud.H"
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 // Get faces (on the new mesh) that have in some way been affected by the
  44 // mesh change. Picks up all faces but those that are between two
  45 // unrefined faces. (Note that of an unchanged face the edge still might be
  46 // split but that does not change any face centre or cell centre.
  47 Foam::labelList Foam::meshRefinement::getChangedFaces
  48 (
  49     const mapPolyMesh& map,
  50     const labelList& oldCellsToRefine
  51 )
  52 {
  53     const polyMesh& mesh = map.mesh();
  54
  55     labelList changedFaces;
  56     {
  57         // Mark any face on a cell which has been added or changed
  58         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  59         // Note that refining a face changes the face centre (for a warped face)
  60         // which changes the cell centre. This again changes the cellcentre-
  61         // cellcentre edge across all faces of the cell.
  62         // Note: this does not happen for unwarped faces but unfortunately
  63         // we don't have this information.
  64
  65         const labelList& faceOwner = mesh.faceOwner();
  66         const labelList& faceNeighbour = mesh.faceNeighbour();
  67         const cellList& cells = mesh.cells();
  68         const label nInternalFaces = mesh.nInternalFaces();
  69
  70         // Mark refined cells on old mesh
  71         PackedList<1> oldRefineCell(map.nOldCells(), 0u);
  72
  73         forAll(oldCellsToRefine, i)
  74         {
  75             oldRefineCell.set(oldCellsToRefine[i], 1u);
  76         }
  77
  78         // Mark refined faces
  79         PackedList<1> refinedInternalFace(nInternalFaces, 0u);
  80
  81         // 1. Internal faces
  82
  83         for (label faceI = 0; faceI < nInternalFaces; faceI++)
  84         {
  85             label oldOwn = map.cellMap()[faceOwner[faceI]];
  86             label oldNei = map.cellMap()[faceNeighbour[faceI]];
  87
  88             if
  89             (
  90                 oldOwn >= 0
  91              && oldRefineCell.get(oldOwn) == 0u
  92              && oldNei >= 0
  93              && oldRefineCell.get(oldNei) == 0u
  94             )
  95             {
  96                 // Unaffected face since both neighbours were not refined.
  97             }
  98             else
  99             {
 100                 refinedInternalFace.set(faceI, 1u);
 101             }
 102         }
 103
 104
 105         // 2. Boundary faces
 106
 107         boolList refinedBoundaryFace(mesh.nFaces()-nInternalFaces, false);
 108
 109         forAll(mesh.boundaryMesh(), patchI)
 110         {
 111             const polyPatch& pp = mesh.boundaryMesh()[patchI];
 112
 113             label faceI = pp.start();
 114
 115             forAll(pp, i)
 116             {
 117                 label oldOwn = map.cellMap()[faceOwner[faceI]];
 118
 119                 if (oldOwn >= 0 && oldRefineCell.get(oldOwn) == 0u)
 120                 {
 121                     // owner did exist and wasn't refined.
 122                 }
 123                 else
 124                 {
 125                     refinedBoundaryFace[faceI-nInternalFaces] = true;
 126                 }
 127                 faceI++;
 128             }
 129         }
 130
 131         // Synchronise refined face status
 132         syncTools::syncBoundaryFaceList
 133         (
 134             mesh,
 135             refinedBoundaryFace,
 136             orEqOp<bool>(),
 137             false
 138         );
 139
 140
 141         // 3. Mark all faces affected by refinement. Refined faces are in
 142         //    - refinedInternalFace
 143         //    - refinedBoundaryFace
 144         boolList changedFace(mesh.nFaces(), false);
 145
 146         forAll(refinedInternalFace, faceI)
 147         {
 148             if (refinedInternalFace.get(faceI) == 1u)
 149             {
 150                 const cell& ownFaces = cells[faceOwner[faceI]];
 151                 forAll(ownFaces, ownI)
 152                 {
 153                     changedFace[ownFaces[ownI]] = true;
 154                 }
 155                 const cell& neiFaces = cells[faceNeighbour[faceI]];
 156                 forAll(neiFaces, neiI)
 157                 {
 158                     changedFace[neiFaces[neiI]] = true;
 159                 }
 160             }
 161         }
 162
 163         forAll(refinedBoundaryFace, i)
 164         {
 165             if (refinedBoundaryFace[i])
 166             {
 167                 const cell& ownFaces = cells[faceOwner[i+nInternalFaces]];
 168                 forAll(ownFaces, ownI)
 169                 {
 170                     changedFace[ownFaces[ownI]] = true;
 171                 }
 172             }
 173         }
 174
 175         syncTools::syncFaceList
 176         (
 177             mesh,
 178             changedFace,
 179             orEqOp<bool>(),
 180             false
 181         );
 182
 183
 184         // Now we have in changedFace marked all affected faces. Pack.
 185         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 186
 187         label nChanged = 0;
 188
 189         forAll(changedFace, faceI)
 190         {
 191             if (changedFace[faceI])
 192             {
 193                 nChanged++;
 194             }
 195         }
 196
 197         changedFaces.setSize(nChanged);
 198         nChanged = 0;
 199
 200         forAll(changedFace, faceI)
 201         {
 202             if (changedFace[faceI])
 203             {
 204                 changedFaces[nChanged++] = faceI;
 205             }
 206         }
 207     }
 208
 209     label nChangedFaces = changedFaces.size();
 210     reduce(nChangedFaces, sumOp<label>());
 211
 212     if (debug)
 213     {
 214         Pout<< "getChangedFaces : Detected "
 215             << " local:" << changedFaces.size()
 216             << " global:" << nChangedFaces
 217             << " changed faces out of " << mesh.globalData().nTotalFaces()
 218             << endl;
 219
 220         faceSet changedFacesSet(mesh, "changedFaces", changedFaces);
 221         Pout<< "getChangedFaces : Writing " << changedFaces.size()
 222             << " changed faces to faceSet " << changedFacesSet.name()
 223             << endl;
 224         changedFacesSet.write();
 225     }
 226
 227     return changedFaces;
 228 }
 229
 230
 231 // Mark cell for refinement (if not already marked). Return false if
 232 // refinelimit hit. Keeps running count (in nRefine) of cells marked for
 233 // refinement
 234 bool Foam::meshRefinement::markForRefine
 235 (
 236     const label markValue,
 237     const label nAllowRefine,
 238
 239     label& cellValue,
 240     label& nRefine
 241 )
 242 {
 243     if (cellValue == -1)
 244     {
 245         cellValue = markValue;
 246         nRefine++;
 247     }
 248
 249     return nRefine <= nAllowRefine;
 250 }
 251
 252
 253 // Calculates list of cells to refine based on intersection with feature edge.
 254 Foam::label Foam::meshRefinement::markFeatureRefinement
 255 (
 256     const point& keepPoint,
 257     const PtrList<featureEdgeMesh>& featureMeshes,
 258     const labelList& featureLevels,
 259     const label nAllowRefine,
 260
 261     labelList& refineCell,
 262     label& nRefine
 263 ) const
 264 {
 265     // We want to refine all cells containing a feature edge.
 266     // - don't want to search over whole mesh
 267     // - don't want to build octree for whole mesh
 268     // - so use tracking to follow the feature edges
 269     //
 270     // 1. find non-manifold points on feature edges (i.e. start of feature edge
 271     //    or 'knot')
 272     // 2. seed particle starting at keepPoint going to this non-manifold point
 273     // 3. track particles to their non-manifold point
 274     //
 275     // 4. track particles across their connected feature edges, marking all
 276     //    visited cells with their level (through trackingData)
 277     // 5. do 4 until all edges have been visited.
 278
 279
 280     // Find all start cells of features. Is done by tracking from keepPoint.
 281     Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
 282
 283     // Create particles on whichever processor holds the keepPoint.
 284     label cellI = mesh_.findCell(keepPoint);
 285
 286     if (cellI != -1)
 287     {
 288         forAll(featureMeshes, featI)
 289         {
 290             const featureEdgeMesh& featureMesh = featureMeshes[featI];
 291             const labelListList& pointEdges = featureMesh.pointEdges();
 292
 293             forAll(pointEdges, pointI)
 294             {
 295                 if (pointEdges[pointI].size() != 2)
 296                 {
 297                     if (debug)
 298                     {
 299                         Pout<< "Adding particle from point:" << pointI
 300                             << " coord:" << featureMesh.points()[pointI]
 301                             << " pEdges:" << pointEdges[pointI]
 302                             << endl;
 303                     }
 304
 305                     // Non-manifold point. Create particle.
 306                     cloud.addParticle
 307                     (
 308                         new trackedParticle
 309                         (
 310                             cloud,
 311                             keepPoint,
 312                             cellI,
 313                             featureMesh.points()[pointI],   // endpos
 314                             featureLevels[featI],           // level
 315                             featI,                          // featureMesh
 316                             pointI                          // end point
 317                         )
 318                     );
 319                 }
 320             }
 321         }
 322     }
 323
 324
 325     // Largest refinement level of any feature passed through
 326     labelList maxFeatureLevel(mesh_.nCells(), -1);
 327
 328     // Database to pass into trackedParticle::move
 329     trackedParticle::trackData td(cloud, maxFeatureLevel);
 330
 331     // Track all particles to their end position.
 332     cloud.move(td);
 333
 334     // Reset level
 335     maxFeatureLevel = -1;
 336
 337     // Whether edge has been visited.
 338     List<PackedList<1> > featureEdgeVisited(featureMeshes.size());
 339
 340     forAll(featureMeshes, featI)
 341     {
 342         featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size());
 343         featureEdgeVisited[featI] = 0u;
 344     }
 345
 346     while (true)
 347     {
 348         label nParticles = 0;
 349
 350         // Make particle follow edge.
 351         forAllIter(Cloud<trackedParticle>, cloud, iter)
 352         {
 353             trackedParticle& tp = iter();
 354
 355             label featI = tp.i();
 356             label pointI = tp.j();
 357
 358             const featureEdgeMesh& featureMesh = featureMeshes[featI];
 359             const labelList& pEdges = featureMesh.pointEdges()[pointI];
 360
 361             // Particle now at pointI. Check connected edges to see which one
 362             // we have to visit now.
 363
 364             bool keepParticle = false;
 365
 366             forAll(pEdges, i)
 367             {
 368                 label edgeI = pEdges[i];
 369
 370                 if (featureEdgeVisited[featI].set(edgeI, 1u))
 371                 {
 372                     // Unvisited edge. Make the particle go to the other point
 373                     // on the edge.
 374
 375                     const edge& e = featureMesh.edges()[edgeI];
 376                     label otherPointI = e.otherVertex(pointI);
 377
 378                     tp.end() = featureMesh.points()[otherPointI];
 379                     tp.j() = otherPointI;
 380                     keepParticle = true;
 381                     break;
 382                 }
 383             }
 384
 385             if (!keepParticle)
 386             {
 387                 // Particle at 'knot' where another particle already has been
 388                 // seeded. Delete particle.
 389                 cloud.deleteParticle(tp);
 390             }
 391             else
 392             {
 393                 // Keep particle
 394                 nParticles++;
 395             }
 396         }
 397
 398         reduce(nParticles, sumOp<label>());
 399         if (nParticles == 0)
 400         {
 401             break;
 402         }
 403
 404         // Track all particles to their end position.
 405         cloud.move(td);
 406     }
 407
 408
 409     // See which cells to refine. maxFeatureLevel will hold highest level
 410     // of any feature edge that passed through.
 411
 412     const labelList& cellLevel = meshCutter_.cellLevel();
 413
 414     label oldNRefine = nRefine;
 415
 416     forAll(maxFeatureLevel, cellI)
 417     {
 418         if (maxFeatureLevel[cellI] > cellLevel[cellI])
 419         {
 420             // Mark
 421             if
 422             (
 423                !markForRefine
 424                 (
 425                     0,                      // surface (n/a)
 426                     nAllowRefine,
 427                     refineCell[cellI],
 428                     nRefine
 429                 )
 430             )
 431             {
 432                 // Reached limit
 433                 break;
 434             }
 435         }
 436     }
 437
 438     if
 439     (
 440         returnReduce(nRefine, sumOp<label>())
 441       > returnReduce(nAllowRefine, sumOp<label>())
 442     )
 443     {
 444         Info<< "Reached refinement limit." << endl;
 445     }
 446
 447     return returnReduce(nRefine-oldNRefine,  sumOp<label>());
 448 }
 449
 450
 451 // Mark cells for non-surface intersection based refinement.
 452 Foam::label Foam::meshRefinement::markInternalRefinement
 453 (
 454     const label nAllowRefine,
 455
 456     labelList& refineCell,
 457     label& nRefine
 458 ) const
 459 {
 460     const labelList& cellLevel = meshCutter_.cellLevel();
 461     const pointField& cellCentres = mesh_.cellCentres();
 462
 463     label oldNRefine = nRefine;
 464
 465     // Collect cells to test
 466     pointField testCc(cellLevel.size()-nRefine);
 467     labelList testLevels(cellLevel.size()-nRefine);
 468     label testI = 0;
 469
 470     forAll(cellLevel, cellI)
 471     {
 472         if (refineCell[cellI] == -1)
 473         {
 474             testCc[testI] = cellCentres[cellI];
 475             testLevels[testI] = cellLevel[cellI];
 476             testI++;
 477         }
 478     }
 479
 480     // Do test to see whether cells is inside/outside shell with higher level
 481     labelList maxLevel;
 482     shells_.findHigherLevel(testCc, testLevels, maxLevel);
 483
 484     // Mark for refinement. Note that we didn't store the original cellID so
 485     // now just reloop in same order.
 486     testI = 0;
 487     forAll(cellLevel, cellI)
 488     {
 489         if (refineCell[cellI] == -1)
 490         {
 491             if (maxLevel[testI] > testLevels[testI])
 492             {
 493                 bool reachedLimit = !markForRefine
 494                 (
 495                     maxLevel[testI],    // mark with any positive value
 496                     nAllowRefine,
 497                     refineCell[cellI],
 498                     nRefine
 499                 );
 500
 501                 if (reachedLimit)
 502                 {
 503                     if (debug)
 504                     {
 505                         Pout<< "Stopped refining internal cells"
 506                             << " since reaching my cell limit of "
 507                             << mesh_.nCells()+7*nRefine << endl;
 508                     }
 509                     break;
 510                 }
 511             }
 512             testI++;
 513         }
 514     }
 515
 516     if
 517     (
 518         returnReduce(nRefine, sumOp<label>())
 519       > returnReduce(nAllowRefine, sumOp<label>())
 520     )
 521     {
 522         Info<< "Reached refinement limit." << endl;
 523     }
 524
 525     return returnReduce(nRefine-oldNRefine, sumOp<label>());
 526 }
 527
 528
 529 // Collect faces that are intersected and whose neighbours aren't yet marked
 530 // for refinement.
 531 Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
 532 (
 533     const labelList& refineCell
 534 ) const
 535 {
 536     labelList testFaces(mesh_.nCells());
 537
 538     label nTest = 0;
 539
 540     forAll(surfaceIndex_, faceI)
 541     {
 542         if (surfaceIndex_[faceI] != -1)
 543         {
 544             label own = mesh_.faceOwner()[faceI];
 545
 546             if (mesh_.isInternalFace(faceI))
 547             {
 548                 label nei = mesh_.faceNeighbour()[faceI];
 549
 550                 if (refineCell[own] == -1 || refineCell[nei] == -1)
 551                 {
 552                     testFaces[nTest++] = faceI;
 553                 }
 554             }
 555             else
 556             {
 557                 if (refineCell[own] == -1)
 558                 {
 559                     testFaces[nTest++] = faceI;
 560                 }
 561             }
 562         }
 563     }
 564     testFaces.setSize(nTest);
 565
 566     return testFaces;
 567 }
 568
 569
 570 // Mark cells for surface intersection based refinement.
 571 Foam::label Foam::meshRefinement::markSurfaceRefinement
 572 (
 573     const label nAllowRefine,
 574     const labelList& neiLevel,
 575     const pointField& neiCc,
 576
 577     labelList& refineCell,
 578     label& nRefine
 579 ) const
 580 {
 581     const labelList& cellLevel = meshCutter_.cellLevel();
 582     const pointField& cellCentres = mesh_.cellCentres();
 583
 584     label oldNRefine = nRefine;
 585
 586     // Use cached surfaceIndex_ to detect if any intersection. If so
 587     // re-intersect to determine level wanted.
 588
 589     // Collect candidate faces
 590     // ~~~~~~~~~~~~~~~~~~~~~~~
 591
 592     labelList testFaces(getRefineCandidateFaces(refineCell));
 593
 594     // Collect segments
 595     // ~~~~~~~~~~~~~~~~
 596
 597     pointField start(testFaces.size());
 598     pointField end(testFaces.size());
 599     labelList minLevel(testFaces.size());
 600
 601     forAll(testFaces, i)
 602     {
 603         label faceI = testFaces[i];
 604
 605         label own = mesh_.faceOwner()[faceI];
 606
 607         if (mesh_.isInternalFace(faceI))
 608         {
 609             label nei = mesh_.faceNeighbour()[faceI];
 610
 611             start[i] = cellCentres[own];
 612             end[i] = cellCentres[nei];
 613             minLevel[i] = min(cellLevel[own], cellLevel[nei]);
 614         }
 615         else
 616         {
 617             label bFaceI = faceI - mesh_.nInternalFaces();
 618
 619             start[i] = cellCentres[own];
 620             end[i] = neiCc[bFaceI];
 621             minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
 622         }
 623     }
 624
 625
 626     // Do test for higher intersections
 627     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 628
 629     labelList surfaceHit;
 630     labelList surfaceMinLevel;
 631     surfaces_.findHigherIntersection
 632     (
 633         start,
 634         end,
 635         minLevel,
 636
 637         surfaceHit,
 638         surfaceMinLevel
 639     );
 640
 641
 642     // Mark cells for refinement
 643     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 644
 645     forAll(testFaces, i)
 646     {
 647         label faceI = testFaces[i];
 648
 649         label surfI = surfaceHit[i];
 650
 651         if (surfI != -1)
 652         {
 653             // Found intersection with surface with higher wanted
 654             // refinement. Check if the level field on that surface
 655             // specifies an even higher level. Note:this is weird. Should
 656             // do the check with the surfaceMinLevel whilst intersecting the
 657             // surfaces?
 658
 659             label own = mesh_.faceOwner()[faceI];
 660
 661             if (surfaceMinLevel[i] > cellLevel[own])
 662             {
 663                 // Owner needs refining
 664                 if
 665                 (
 666                    !markForRefine
 667                     (
 668                         surfI,
 669                         nAllowRefine,
 670                         refineCell[own],
 671                         nRefine
 672                     )
 673                 )
 674                 {
 675                     break;
 676                 }
 677             }
 678
 679             if (mesh_.isInternalFace(faceI))
 680             {
 681                 label nei = mesh_.faceNeighbour()[faceI];
 682                 if (surfaceMinLevel[i] > cellLevel[nei])
 683                 {
 684                     // Neighbour needs refining
 685                     if
 686                     (
 687                        !markForRefine
 688                         (
 689                             surfI,
 690                             nAllowRefine,
 691                             refineCell[nei],
 692                             nRefine
 693                         )
 694                     )
 695                     {
 696                         break;
 697                     }
 698                 }
 699             }
 700         }
 701     }
 702
 703     if
 704     (
 705         returnReduce(nRefine, sumOp<label>())
 706       > returnReduce(nAllowRefine, sumOp<label>())
 707     )
 708     {
 709         Info<< "Reached refinement limit." << endl;
 710     }
 711
 712     return returnReduce(nRefine-oldNRefine, sumOp<label>());
 713 }
 714
 715
 716 // Checks if multiple intersections of a cell (by a surface with a higher
 717 // max than the cell level) and if so if the normals at these intersections
 718 // make a large angle.
 719 // Returns false if the nRefine limit has been reached, true otherwise.
 720 bool Foam::meshRefinement::checkCurvature
 721 (
 722     const scalar curvature,
 723     const label nAllowRefine,
 724
 725     const label surfaceLevel,   // current intersection max level
 726     const vector& surfaceNormal,// current intersection normal
 727
 728     const label cellI,
 729
 730     label& cellMaxLevel,        // cached max surface level for this cell
 731     vector& cellMaxNormal,      // cached surface normal for this cell
 732
 733     labelList& refineCell,
 734     label& nRefine
 735 ) const
 736 {
 737     const labelList& cellLevel = meshCutter_.cellLevel();
 738
 739     // Test if surface applicable
 740     if (surfaceLevel > cellLevel[cellI])
 741     {
 742         if (cellMaxLevel == -1)
 743         {
 744             // First visit of cell. Store
 745             cellMaxLevel = surfaceLevel;
 746             cellMaxNormal = surfaceNormal;
 747         }
 748         else
 749         {
 750             // Second or more visit. Check curvature.
 751             if ((cellMaxNormal & surfaceNormal) < curvature)
 752             {
 753                 return markForRefine
 754                 (
 755                     surfaceLevel,   // mark with any non-neg number.
 756                     nAllowRefine,
 757                     refineCell[cellI],
 758                     nRefine
 759                 );
 760             }
 761
 762             // Set normal to that of highest surface. Not really necessary
 763             // over here but we reuse cellMax info when doing coupled faces.
 764             if (surfaceLevel > cellMaxLevel)
 765             {
 766                 cellMaxLevel = surfaceLevel;
 767                 cellMaxNormal = surfaceNormal;
 768             }
 769         }
 770     }
 771
 772     // Did not reach refinement limit.
 773     return true;
 774 }
 775
 776
 777 // Mark cells for surface curvature based refinement.
 778 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 779 (
 780     const scalar curvature,
 781     const label nAllowRefine,
 782     const labelList& neiLevel,
 783     const pointField& neiCc,
 784
 785     labelList& refineCell,
 786     label& nRefine
 787 ) const
 788 {
 789     const labelList& cellLevel = meshCutter_.cellLevel();
 790     const pointField& cellCentres = mesh_.cellCentres();
 791
 792     label oldNRefine = nRefine;
 793
 794     // 1. local test: any cell on more than one surface gets refined
 795     // (if its current level is < max of the surface max level)
 796
 797     // 2. 'global' test: any cell on only one surface with a neighbour
 798     // on a different surface gets refined (if its current level etc.)
 799
 800
 801     // Collect candidate faces (i.e. intersecting any surface and
 802     // owner/neighbour not yet refined.
 803     labelList testFaces(getRefineCandidateFaces(refineCell));
 804
 805     // Collect segments
 806     pointField start(testFaces.size());
 807     pointField end(testFaces.size());
 808     labelList minLevel(testFaces.size());
 809
 810     forAll(testFaces, i)
 811     {
 812         label faceI = testFaces[i];
 813
 814         label own = mesh_.faceOwner()[faceI];
 815
 816         if (mesh_.isInternalFace(faceI))
 817         {
 818             label nei = mesh_.faceNeighbour()[faceI];
 819
 820             start[i] = cellCentres[own];
 821             end[i] = cellCentres[nei];
 822             minLevel[i] = min(cellLevel[own], cellLevel[nei]);
 823         }
 824         else
 825         {
 826             label bFaceI = faceI - mesh_.nInternalFaces();
 827
 828             start[i] = cellCentres[own];
 829             end[i] = neiCc[bFaceI];
 830             minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
 831         }
 832     }
 833
 834     // Test for all intersections (with surfaces of higher max level than
 835     // minLevel) and cache per cell the max surface level and the local normal
 836     // on that surface.
 837     labelList cellMaxLevel(mesh_.nCells(), -1);
 838     vectorField cellMaxNormal(mesh_.nCells());
 839
 840     {
 841         // Per segment the normals of the surfaces
 842         List<vectorList> surfaceNormal;
 843         // Per segment the list of levels of the surfaces
 844         labelListList surfaceLevel;
 845
 846         surfaces_.findAllHigherIntersections
 847         (
 848             start,
 849             end,
 850             minLevel,           // max level of surface has to be bigger
 851                                 // than min level of neighbouring cells
 852             surfaceNormal,
 853             surfaceLevel
 854         );
 855         // Clear out unnecessary data
 856         start.clear();
 857         end.clear();
 858         minLevel.clear();
 859
 860         // Extract per cell information on the surface with the highest max
 861         forAll(testFaces, i)
 862         {
 863             label faceI = testFaces[i];
 864             label own = mesh_.faceOwner()[faceI];
 865
 866             const vectorList& fNormals = surfaceNormal[i];
 867             const labelList& fLevels = surfaceLevel[i];
 868
 869             forAll(fLevels, hitI)
 870             {
 871                 checkCurvature
 872                 (
 873                     curvature,
 874                     nAllowRefine,
 875
 876                     fLevels[hitI],
 877                     fNormals[hitI],
 878
 879                     own,
 880                     cellMaxLevel[own],
 881                     cellMaxNormal[own],
 882
 883                     refineCell,
 884                     nRefine
 885                 );
 886             }
 887
 888             if (mesh_.isInternalFace(faceI))
 889             {
 890                 label nei = mesh_.faceNeighbour()[faceI];
 891
 892                 forAll(fLevels, hitI)
 893                 {
 894                     checkCurvature
 895                     (
 896                         curvature,
 897                         nAllowRefine,
 898
 899                         fLevels[hitI],
 900                         fNormals[hitI],
 901
 902                         nei,
 903                         cellMaxLevel[nei],
 904                         cellMaxNormal[nei],
 905
 906                         refineCell,
 907                         nRefine
 908                     );
 909                 }
 910             }
 911         }
 912     }
 913
 914     // 2. Find out a measure of surface curvature and region edges.
 915     // Send over surface region and surface normal to neighbour cell.
 916
 917     labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 918     vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
 919
 920     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 921     {
 922         label bFaceI = faceI-mesh_.nInternalFaces();
 923         label own = mesh_.faceOwner()[faceI];
 924
 925         neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
 926         neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
 927     }
 928     syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel, false);
 929     syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal, false);
 930
 931     // Loop over all faces. Could only be checkFaces.. except if they're coupled
 932
 933     // Internal faces
 934     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 935     {
 936         label own = mesh_.faceOwner()[faceI];
 937         label nei = mesh_.faceNeighbour()[faceI];
 938
 939         if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
 940         {
 941             // Have valid data on both sides. Check curvature.
 942             if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
 943             {
 944                 // See which side to refine
 945                 if (cellLevel[own] < cellMaxLevel[own])
 946                 {
 947                     if
 948                     (
 949                         !markForRefine
 950                         (
 951                             cellMaxLevel[own],
 952                             nAllowRefine,
 953                             refineCell[own],
 954                             nRefine
 955                         )
 956                     )
 957                     {
 958                         if (debug)
 959                         {
 960                             Pout<< "Stopped refining since reaching my cell"
 961                                 << " limit of " << mesh_.nCells()+7*nRefine
 962                                 << endl;
 963                         }
 964                         break;
 965                     }
 966                 }
 967
 968                 if (cellLevel[nei] < cellMaxLevel[nei])
 969                 {
 970                     if
 971                     (
 972                         !markForRefine
 973                         (
 974                             cellMaxLevel[nei],
 975                             nAllowRefine,
 976                             refineCell[nei],
 977                             nRefine
 978                         )
 979                     )
 980                     {
 981                         if (debug)
 982                         {
 983                             Pout<< "Stopped refining since reaching my cell"
 984                                 << " limit of " << mesh_.nCells()+7*nRefine
 985                                 << endl;
 986                         }
 987                         break;
 988                     }
 989                 }
 990             }
 991         }
 992     }
 993     // Boundary faces
 994     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 995     {
 996         label own = mesh_.faceOwner()[faceI];
 997         label bFaceI = faceI - mesh_.nInternalFaces();
 998
 999         if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
1000         {
1001             // Have valid data on both sides. Check curvature.
1002             if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
1003             {
1004                 if
1005                 (
1006                     !markForRefine
1007                     (
1008                         cellMaxLevel[own],
1009                         nAllowRefine,
1010                         refineCell[own],
1011                         nRefine
1012                     )
1013                 )
1014                 {
1015                     if (debug)
1016                     {
1017                         Pout<< "Stopped refining since reaching my cell"
1018                             << " limit of " << mesh_.nCells()+7*nRefine
1019                             << endl;
1020                     }
1021                     break;
1022                 }
1023             }
1024         }
1025     }
1026
1027     if
1028     (
1029         returnReduce(nRefine, sumOp<label>())
1030       > returnReduce(nAllowRefine, sumOp<label>())
1031     )
1032     {
1033         Info<< "Reached refinement limit." << endl;
1034     }
1035
1036     return returnReduce(nRefine-oldNRefine, sumOp<label>());
1037 }
1038
1039
1040 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1041
1042 // Calculate list of cells to refine. Gets for any edge (start - end)
1043 // whether it intersects the surface. Does more accurate test and checks
1044 // the wanted level on the surface intersected.
1045 // Does approximate precalculation of how many cells can be refined before
1046 // hitting overall limit maxGlobalCells.
1047 Foam::labelList Foam::meshRefinement::refineCandidates
1048 (
1049     const point& keepPoint,
1050     const scalar curvature,
1051
1052     const PtrList<featureEdgeMesh>& featureMeshes,
1053     const labelList& featureLevels,
1054
1055     const bool featureRefinement,
1056     const bool internalRefinement,
1057     const bool surfaceRefinement,
1058     const bool curvatureRefinement,
1059     const label maxGlobalCells,
1060     const label maxLocalCells
1061 ) const
1062 {
1063     label totNCells = mesh_.globalData().nTotalCells();
1064
1065     labelList cellsToRefine;
1066
1067     if (totNCells >= maxGlobalCells)
1068     {
1069         Info<< "No cells marked for refinement since reached limit "
1070             << maxGlobalCells << '.' << endl;
1071     }
1072     else
1073     {
1074         // Every cell I refine adds 7 cells. Estimate the number of cells
1075         // I am allowed to refine.
1076         // Assume perfect distribution so can only refine as much the fraction
1077         // of the mesh I hold. This prediction step prevents us having to do
1078         // lots of reduces to keep count of the total number of cells selected
1079         // for refinement.
1080
1081         //scalar fraction = scalar(mesh_.nCells())/totNCells;
1082         //label myMaxCells = label(maxGlobalCells*fraction);
1083         //label nAllowRefine = (myMaxCells - mesh_.nCells())/7;
1084         ////label nAllowRefine = (maxLocalCells - mesh_.nCells())/7;
1085         //
1086         //Pout<< "refineCandidates:" << nl
1087         //    << "    total cells:" << totNCells << nl
1088         //    << "    local cells:" << mesh_.nCells() << nl
1089         //    << "    local fraction:" << fraction << nl
1090         //    << "    local allowable cells:" << myMaxCells << nl
1091         //    << "    local allowable refinement:" << nAllowRefine << nl
1092         //    << endl;
1093
1094         //- Disable refinement shortcut
1095         label nAllowRefine = labelMax;
1096
1097         // Marked for refinement (>= 0) or not (-1). Actual value is the
1098         // index of the surface it intersects.
1099         labelList refineCell(mesh_.nCells(), -1);
1100         label nRefine = 0;
1101
1102
1103         // Swap neighbouring cell centres and cell level
1104         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1105         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1106         calcNeighbourData(neiLevel, neiCc);
1107
1108
1109
1110         // Cells pierced by feature lines
1111         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1112
1113         if (featureRefinement)
1114         {
1115             label nFeatures = markFeatureRefinement
1116             (
1117                 keepPoint,
1118                 featureMeshes,
1119                 featureLevels,
1120                 nAllowRefine,
1121
1122                 refineCell,
1123                 nRefine
1124             );
1125
1126             Info<< "Marked for refinement due to explicit features    : "
1127                 << nFeatures << " cells."  << endl;
1128         }
1129
1130         // Inside refinement shells
1131         // ~~~~~~~~~~~~~~~~~~~~~~~~
1132
1133         if (internalRefinement)
1134         {
1135             label nShell = markInternalRefinement
1136             (
1137                 nAllowRefine,
1138
1139                 refineCell,
1140                 nRefine
1141             );
1142             Info<< "Marked for refinement due to refinement shells    : "
1143                 << nShell << " cells."  << endl;
1144         }
1145
1146         // Refinement based on intersection of surface
1147         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1148
1149         if (surfaceRefinement)
1150         {
1151             label nSurf = markSurfaceRefinement
1152             (
1153                 nAllowRefine,
1154                 neiLevel,
1155                 neiCc,
1156
1157                 refineCell,
1158                 nRefine
1159             );
1160             Info<< "Marked for refinement due to surface intersection : "
1161                 << nSurf << " cells."  << endl;
1162         }
1163
1164         // Refinement based on curvature of surface
1165         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166
1167         if (curvatureRefinement && (curvature >= -1 && curvature <= 1))
1168         {
1169             label nCurv = markSurfaceCurvatureRefinement
1170             (
1171                 curvature,
1172                 nAllowRefine,
1173                 neiLevel,
1174                 neiCc,
1175
1176                 refineCell,
1177                 nRefine
1178             );
1179             Info<< "Marked for refinement due to curvature/regions    : "
1180                 << nCurv << " cells."  << endl;
1181         }
1182
1183         // Pack cells-to-refine
1184         // ~~~~~~~~~~~~~~~~~~~~
1185
1186         cellsToRefine.setSize(nRefine);
1187         nRefine = 0;
1188
1189         forAll(refineCell, cellI)
1190         {
1191             if (refineCell[cellI] != -1)
1192             {
1193                 cellsToRefine[nRefine++] = cellI;
1194             }
1195         }
1196     }
1197
1198     return cellsToRefine;
1199 }
1200
1201
1202 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
1203 (
1204     const labelList& cellsToRefine
1205 )
1206 {
1207     // Mesh changing engine.
1208     polyTopoChange meshMod(mesh_);
1209
1210     // Play refinement commands into mesh changer.
1211     meshCutter_.setRefinement(cellsToRefine, meshMod);
1212
1213     // Create mesh (no inflation), return map from old to new mesh.
1214     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false);
1215
1216     // Update fields
1217     mesh_.updateMesh(map);
1218
1219     // Optionally inflate mesh
1220     if (map().hasMotionPoints())
1221     {
1222         mesh_.movePoints(map().preMotionPoints());
1223     }
1224
1225     // Update intersection info
1226     updateMesh(map, getChangedFaces(map, cellsToRefine));
1227
1228     return map;
1229 }
1230
1231
1232 // Do refinement of consistent set of cells followed by truncation and
1233 // load balancing.
1234 Foam::autoPtr<Foam::mapDistributePolyMesh>
1235  Foam::meshRefinement::refineAndBalance
1236 (
1237     const string& msg,
1238     decompositionMethod& decomposer,
1239     fvMeshDistribute& distributor,
1240     const labelList& cellsToRefine
1241 )
1242 {
1243     // Do all refinement
1244     refine(cellsToRefine);
1245
1246     if (debug)
1247     {
1248         Pout<< "Writing refined but unbalanced " << msg
1249             << " mesh to time " << mesh_.time().timeName() << endl;
1250         write
1251         (
1252             debug,
1253             mesh_.time().path()
1254            /mesh_.time().timeName()
1255         );
1256         Pout<< "Dumped debug data in = "
1257             << mesh_.time().cpuTimeIncrement() << " s" << endl;
1258
1259         // test all is still synced across proc patches
1260         checkData();
1261     }
1262
1263     Info<< "Refined mesh in = "
1264         << mesh_.time().cpuTimeIncrement() << " s" << endl;
1265     printMeshInfo(debug, "After refinement " + msg);
1266
1267
1268     // Load balancing
1269     // ~~~~~~~~~~~~~~
1270
1271     autoPtr<mapDistributePolyMesh> distMap;
1272
1273     if (Pstream::nProcs() > 1)
1274     {
1275         labelList distribution(decomposer.decompose(mesh_.cellCentres()));
1276         // Get distribution such that baffle faces stay internal to the
1277         // processor.
1278         //labelList distribution(decomposePreserveBaffles(decomposer));
1279
1280         if (debug)
1281         {
1282             Pout<< "Wanted distribution:"
1283                 << distributor.countCells(distribution)
1284                 << endl;
1285         }
1286         // Do actual sending/receiving of mesh
1287         distMap = distributor.distribute(distribution);
1288
1289         // Update numbering of meshRefiner
1290         distribute(distMap);
1291
1292         Info<< "Balanced mesh in = "
1293             << mesh_.time().cpuTimeIncrement() << " s" << endl;
1294
1295         printMeshInfo(debug, "After balancing " + msg);
1296
1297
1298         if (debug)
1299         {
1300             Pout<< "Writing " << msg
1301                 << " mesh to time " << mesh_.time().timeName() << endl;
1302             write
1303             (
1304                 debug,
1305                 mesh_.time().path()
1306                /mesh_.time().timeName()
1307             );
1308             Pout<< "Dumped debug data in = "
1309                 << mesh_.time().cpuTimeIncrement() << " s" << endl;
1310
1311             // test all is still synced across proc patches
1312             checkData();
1313         }
1314     }
1315
1316     return distMap;
1317 }
1318
1319
1320 // ************************************************************************* //