src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*----------------------------------------------------------------------------*/
  25
  26 #include "meshRefinement.H"
  27 #include "volMesh.H"
  28 #include "volFields.H"
  29 #include "surfaceMesh.H"
  30 #include "syncTools.H"
  31 #include "Time.H"
  32 #include "refinementHistory.H"
  33 #include "refinementSurfaces.H"
  34 #include "refinementFeatures.H"
  35 #include "decompositionMethod.H"
  36 #include "regionSplit.H"
  37 #include "fvMeshDistribute.H"
  38 #include "indirectPrimitivePatch.H"
  39 #include "polyTopoChange.H"
  40 #include "removeCells.H"
  41 #include "mapDistributePolyMesh.H"
  42 #include "localPointRegion.H"
  43 #include "pointMesh.H"
  44 #include "pointFields.H"
  45 #include "slipPointPatchFields.H"
  46 #include "fixedValuePointPatchFields.H"
  47 #include "calculatedPointPatchFields.H"
  48 #include "cyclicSlipPointPatchFields.H"
  49 #include "processorPointPatch.H"
  50 #include "globalIndex.H"
  51 #include "meshTools.H"
  52 #include "OFstream.H"
  53 #include "geomDecomp.H"
  54 #include "Random.H"
  55 #include "searchableSurfaces.H"
  56 #include "treeBoundBox.H"
  57 #include "zeroGradientFvPatchFields.H"
  58
  59 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  60
  61 namespace Foam
  62 {
  63     defineTypeNameAndDebug(meshRefinement, 0);
  64 }
  65
  66 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  67
  68 void Foam::meshRefinement::calcNeighbourData
  69 (
  70     labelList& neiLevel,
  71     pointField& neiCc
  72 )  const
  73 {
  74     const labelList& cellLevel = meshCutter_.cellLevel();
  75     const pointField& cellCentres = mesh_.cellCentres();
  76
  77     label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
  78
  79     if (neiLevel.size() != nBoundaryFaces || neiCc.size() != nBoundaryFaces)
  80     {
  81         FatalErrorIn("meshRefinement::calcNeighbour(..)") << "nBoundaries:"
  82             << nBoundaryFaces << " neiLevel:" << neiLevel.size()
  83             << abort(FatalError);
  84     }
  85
  86     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  87
  88     labelHashSet addedPatchIDSet(meshedPatches());
  89
  90     forAll(patches, patchI)
  91     {
  92         const polyPatch& pp = patches[patchI];
  93
  94         const labelUList& faceCells = pp.faceCells();
  95         const vectorField::subField faceCentres = pp.faceCentres();
  96         const vectorField::subField faceAreas = pp.faceAreas();
  97
  98         label bFaceI = pp.start()-mesh_.nInternalFaces();
  99
 100         if (pp.coupled())
 101         {
 102             forAll(faceCells, i)
 103             {
 104                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
 105                 neiCc[bFaceI] = cellCentres[faceCells[i]];
 106                 bFaceI++;
 107             }
 108         }
 109         else if (addedPatchIDSet.found(patchI))
 110         {
 111             // Face was introduced from cell-cell intersection. Try to
 112             // reconstruct other side cell(centre). Three possibilities:
 113             // - cells same size.
 114             // - preserved cell smaller. Not handled.
 115             // - preserved cell larger.
 116             forAll(faceCells, i)
 117             {
 118                 // Extrapolate the face centre.
 119                 vector fn = faceAreas[i];
 120                 fn /= mag(fn)+VSMALL;
 121
 122                 label own = faceCells[i];
 123                 label ownLevel = cellLevel[own];
 124                 label faceLevel = meshCutter_.getAnchorLevel(pp.start()+i);
 125
 126                 // Normal distance from face centre to cell centre
 127                 scalar d = ((faceCentres[i] - cellCentres[own]) & fn);
 128                 if (faceLevel > ownLevel)
 129                 {
 130                     // Other cell more refined. Adjust normal distance
 131                     d *= 0.5;
 132                 }
 133                 neiLevel[bFaceI] = cellLevel[ownLevel];
 134                 // Calculate other cell centre by extrapolation
 135                 neiCc[bFaceI] = faceCentres[i] + d*fn;
 136                 bFaceI++;
 137             }
 138         }
 139         else
 140         {
 141             forAll(faceCells, i)
 142             {
 143                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
 144                 neiCc[bFaceI] = faceCentres[i];
 145                 bFaceI++;
 146             }
 147         }
 148     }
 149
 150     // Swap coupled boundaries. Apply separation to cc since is coordinate.
 151     syncTools::swapBoundaryFacePositions(mesh_, neiCc);
 152     syncTools::swapBoundaryFaceList(mesh_, neiLevel);
 153 }
 154
 155
 156 // Find intersections of edges (given by their two endpoints) with surfaces.
 157 // Returns first intersection if there are more than one.
 158 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 159 {
 160     const pointField& cellCentres = mesh_.cellCentres();
 161
 162     // Stats on edges to test. Count proc faces only once.
 163     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 164
 165     {
 166         label nMasterFaces = 0;
 167         forAll(isMasterFace, faceI)
 168         {
 169             if (isMasterFace.get(faceI) == 1)
 170             {
 171                 nMasterFaces++;
 172             }
 173         }
 174         reduce(nMasterFaces, sumOp<label>());
 175
 176         label nChangedFaces = 0;
 177         forAll(changedFaces, i)
 178         {
 179             if (isMasterFace.get(changedFaces[i]) == 1)
 180             {
 181                 nChangedFaces++;
 182             }
 183         }
 184         reduce(nChangedFaces, sumOp<label>());
 185
 186         Info<< "Edge intersection testing:" << nl
 187             << "    Number of edges             : " << nMasterFaces << nl
 188             << "    Number of edges to retest   : " << nChangedFaces
 189             << endl;
 190     }
 191
 192
 193     // Get boundary face centre and level. Coupled aware.
 194     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 195     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 196     calcNeighbourData(neiLevel, neiCc);
 197
 198     // Collect segments we want to test for
 199     pointField start(changedFaces.size());
 200     pointField end(changedFaces.size());
 201
 202     forAll(changedFaces, i)
 203     {
 204         label faceI = changedFaces[i];
 205         label own = mesh_.faceOwner()[faceI];
 206
 207         start[i] = cellCentres[own];
 208         if (mesh_.isInternalFace(faceI))
 209         {
 210             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
 211         }
 212         else
 213         {
 214             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
 215         }
 216     }
 217
 218     // Extend segments a bit
 219     {
 220         const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
 221         start -= smallVec;
 222         end += smallVec;
 223     }
 224
 225
 226     // Do tests in one go
 227     labelList surfaceHit;
 228     {
 229         labelList surfaceLevel;
 230         surfaces_.findHigherIntersection
 231         (
 232             start,
 233             end,
 234             labelList(start.size(), -1),    // accept any intersection
 235             surfaceHit,
 236             surfaceLevel
 237         );
 238     }
 239
 240     // Keep just surface hit
 241     forAll(surfaceHit, i)
 242     {
 243         surfaceIndex_[changedFaces[i]] = surfaceHit[i];
 244     }
 245
 246     // Make sure both sides have same information. This should be
 247     // case in general since same vectors but just to make sure.
 248     syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>());
 249
 250     label nHits = countHits();
 251     label nTotHits = returnReduce(nHits, sumOp<label>());
 252
 253     Info<< "    Number of intersected edges : " << nTotHits << endl;
 254
 255     // Set files to same time as mesh
 256     setInstance(mesh_.facesInstance());
 257 }
 258
 259
 260 void Foam::meshRefinement::checkData()
 261 {
 262     Pout<< "meshRefinement::checkData() : Checking refinement structure."
 263         << endl;
 264     meshCutter_.checkMesh();
 265
 266     Pout<< "meshRefinement::checkData() : Checking refinement levels."
 267         << endl;
 268     meshCutter_.checkRefinementLevels(1, labelList(0));
 269
 270
 271     label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
 272
 273     Pout<< "meshRefinement::checkData() : Checking synchronization."
 274         << endl;
 275
 276     // Check face centres
 277     {
 278         // Boundary face centres
 279         pointField::subList boundaryFc
 280         (
 281             mesh_.faceCentres(),
 282             nBnd,
 283             mesh_.nInternalFaces()
 284         );
 285
 286         // Get neighbouring face centres
 287         pointField neiBoundaryFc(boundaryFc);
 288         syncTools::syncBoundaryFacePositions
 289         (
 290             mesh_,
 291             neiBoundaryFc,
 292             eqOp<point>()
 293         );
 294
 295         // Compare
 296         testSyncBoundaryFaceList
 297         (
 298             mergeDistance_,
 299             "testing faceCentres : ",
 300             boundaryFc,
 301             neiBoundaryFc
 302         );
 303     }
 304     // Check meshRefinement
 305     {
 306         // Get boundary face centre and level. Coupled aware.
 307         labelList neiLevel(nBnd);
 308         pointField neiCc(nBnd);
 309         calcNeighbourData(neiLevel, neiCc);
 310
 311         // Collect segments we want to test for
 312         pointField start(mesh_.nFaces());
 313         pointField end(mesh_.nFaces());
 314
 315         forAll(start, faceI)
 316         {
 317             start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
 318
 319             if (mesh_.isInternalFace(faceI))
 320             {
 321                 end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
 322             }
 323             else
 324             {
 325                 end[faceI] = neiCc[faceI-mesh_.nInternalFaces()];
 326             }
 327         }
 328
 329         // Extend segments a bit
 330         {
 331             const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
 332             start -= smallVec;
 333             end += smallVec;
 334         }
 335
 336
 337         // Do tests in one go
 338         labelList surfaceHit;
 339         {
 340             labelList surfaceLevel;
 341             surfaces_.findHigherIntersection
 342             (
 343                 start,
 344                 end,
 345                 labelList(start.size(), -1),    // accept any intersection
 346                 surfaceHit,
 347                 surfaceLevel
 348             );
 349         }
 350         // Get the coupled hit
 351         labelList neiHit
 352         (
 353             SubList<label>
 354             (
 355                 surfaceHit,
 356                 nBnd,
 357                 mesh_.nInternalFaces()
 358             )
 359         );
 360         syncTools::swapBoundaryFaceList(mesh_, neiHit);
 361
 362         // Check
 363         forAll(surfaceHit, faceI)
 364         {
 365             if (surfaceIndex_[faceI] != surfaceHit[faceI])
 366             {
 367                 if (mesh_.isInternalFace(faceI))
 368                 {
 369                     WarningIn("meshRefinement::checkData()")
 370                         << "Internal face:" << faceI
 371                         << " fc:" << mesh_.faceCentres()[faceI]
 372                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 373                         << " current:" << surfaceHit[faceI]
 374                         << " ownCc:"
 375                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 376                         << " neiCc:"
 377                         << mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
 378                         << endl;
 379                 }
 380                 else if
 381                 (
 382                     surfaceIndex_[faceI]
 383                  != neiHit[faceI-mesh_.nInternalFaces()]
 384                 )
 385                 {
 386                     WarningIn("meshRefinement::checkData()")
 387                         << "Boundary face:" << faceI
 388                         << " fc:" << mesh_.faceCentres()[faceI]
 389                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 390                         << " current:" << surfaceHit[faceI]
 391                         << " ownCc:"
 392                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 393                         << " end:" << end[faceI]
 394                         << endl;
 395                 }
 396             }
 397         }
 398     }
 399     {
 400         labelList::subList boundarySurface
 401         (
 402             surfaceIndex_,
 403             mesh_.nFaces()-mesh_.nInternalFaces(),
 404             mesh_.nInternalFaces()
 405         );
 406
 407         labelList neiBoundarySurface(boundarySurface);
 408         syncTools::swapBoundaryFaceList
 409         (
 410             mesh_,
 411             neiBoundarySurface
 412         );
 413
 414         // Compare
 415         testSyncBoundaryFaceList
 416         (
 417             0,                              // tolerance
 418             "testing surfaceIndex() : ",
 419             boundarySurface,
 420             neiBoundarySurface
 421         );
 422     }
 423
 424
 425     // Find duplicate faces
 426     Pout<< "meshRefinement::checkData() : Counting duplicate faces."
 427         << endl;
 428
 429     labelList duplicateFace
 430     (
 431         localPointRegion::findDuplicateFaces
 432         (
 433             mesh_,
 434             identity(mesh_.nFaces()-mesh_.nInternalFaces())
 435           + mesh_.nInternalFaces()
 436         )
 437     );
 438
 439     // Count
 440     {
 441         label nDup = 0;
 442
 443         forAll(duplicateFace, i)
 444         {
 445             if (duplicateFace[i] != -1)
 446             {
 447                 nDup++;
 448             }
 449         }
 450         nDup /= 2;  // will have counted both faces of duplicate
 451         Pout<< "meshRefinement::checkData() : Found " << nDup
 452             << " duplicate pairs of faces." << endl;
 453     }
 454 }
 455
 456
 457 void Foam::meshRefinement::setInstance(const fileName& inst)
 458 {
 459     meshCutter_.setInstance(inst);
 460     surfaceIndex_.instance() = inst;
 461 }
 462
 463
 464 // Remove cells. Put exposedFaces (output of getExposedFaces(cellsToRemove))
 465 // into exposedPatchIDs.
 466 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 467 (
 468     const labelList& cellsToRemove,
 469     const labelList& exposedFaces,
 470     const labelList& exposedPatchIDs,
 471     removeCells& cellRemover
 472 )
 473 {
 474     polyTopoChange meshMod(mesh_);
 475
 476     // Arbitrary: put exposed faces into last patch.
 477     cellRemover.setRefinement
 478     (
 479         cellsToRemove,
 480         exposedFaces,
 481         exposedPatchIDs,
 482         meshMod
 483     );
 484
 485     // Change the mesh (no inflation)
 486     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 487
 488     // Update fields
 489     mesh_.updateMesh(map);
 490
 491     // Move mesh (since morphing might not do this)
 492     if (map().hasMotionPoints())
 493     {
 494         mesh_.movePoints(map().preMotionPoints());
 495     }
 496     else
 497     {
 498         // Delete mesh volumes. No other way to do this?
 499         mesh_.clearOut();
 500     }
 501
 502     // Reset the instance for if in overwrite mode
 503     mesh_.setInstance(timeName());
 504     setInstance(mesh_.facesInstance());
 505
 506     // Update local mesh data
 507     cellRemover.updateMesh(map);
 508
 509     // Update intersections. Recalculate intersections for exposed faces.
 510     labelList newExposedFaces = renumber
 511     (
 512         map().reverseFaceMap(),
 513         exposedFaces
 514     );
 515
 516     //Pout<< "removeCells : updating intersections for "
 517     //    << newExposedFaces.size() << " newly exposed faces." << endl;
 518
 519     updateMesh(map, newExposedFaces);
 520
 521     return map;
 522 }
 523
 524
 525 // Determine for multi-processor regions the lowest numbered cell on the lowest
 526 // numbered processor.
 527 void Foam::meshRefinement::getCoupledRegionMaster
 528 (
 529     const globalIndex& globalCells,
 530     const boolList& blockedFace,
 531     const regionSplit& globalRegion,
 532     Map<label>& regionToMaster
 533 ) const
 534 {
 535     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 536
 537     forAll(patches, patchI)
 538     {
 539         const polyPatch& pp = patches[patchI];
 540
 541         if (isA<processorPolyPatch>(pp))
 542         {
 543             forAll(pp, i)
 544             {
 545                 label faceI = pp.start()+i;
 546
 547                 if (!blockedFace[faceI])
 548                 {
 549                     // Only if there is a connection to the neighbour
 550                     // will there be a multi-domain region; if not through
 551                     // this face then through another.
 552
 553                     label cellI = mesh_.faceOwner()[faceI];
 554                     label globalCellI = globalCells.toGlobal(cellI);
 555
 556                     Map<label>::iterator iter =
 557                         regionToMaster.find(globalRegion[cellI]);
 558
 559                     if (iter != regionToMaster.end())
 560                     {
 561                         label master = iter();
 562                         iter() = min(master, globalCellI);
 563                     }
 564                     else
 565                     {
 566                         regionToMaster.insert
 567                         (
 568                             globalRegion[cellI],
 569                             globalCellI
 570                         );
 571                     }
 572                 }
 573             }
 574         }
 575     }
 576
 577     // Do reduction
 578     Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
 579     Pstream::mapCombineScatter(regionToMaster);
 580 }
 581
 582
 583 void Foam::meshRefinement::calcLocalRegions
 584 (
 585     const globalIndex& globalCells,
 586     const labelList& globalRegion,
 587     const Map<label>& coupledRegionToMaster,
 588     const scalarField& cellWeights,
 589
 590     Map<label>& globalToLocalRegion,
 591     pointField& localPoints,
 592     scalarField& localWeights
 593 ) const
 594 {
 595     globalToLocalRegion.resize(globalRegion.size());
 596     DynamicList<point> localCc(globalRegion.size()/2);
 597     DynamicList<scalar> localWts(globalRegion.size()/2);
 598
 599     forAll(globalRegion, cellI)
 600     {
 601         Map<label>::const_iterator fndMaster =
 602             coupledRegionToMaster.find(globalRegion[cellI]);
 603
 604         if (fndMaster != coupledRegionToMaster.end())
 605         {
 606             // Multi-processor region.
 607             if (globalCells.toGlobal(cellI) == fndMaster())
 608             {
 609                 // I am master. Allocate region for me.
 610                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
 611                 localCc.append(mesh_.cellCentres()[cellI]);
 612                 localWts.append(cellWeights[cellI]);
 613             }
 614         }
 615         else
 616         {
 617             // Single processor region.
 618             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
 619             {
 620                 localCc.append(mesh_.cellCentres()[cellI]);
 621                 localWts.append(cellWeights[cellI]);
 622             }
 623         }
 624     }
 625
 626     localPoints.transfer(localCc);
 627     localWeights.transfer(localWts);
 628
 629     if (localPoints.size() != globalToLocalRegion.size())
 630     {
 631         FatalErrorIn("calcLocalRegions(..)")
 632             << "localPoints:" << localPoints.size()
 633             << " globalToLocalRegion:" << globalToLocalRegion.size()
 634             << exit(FatalError);
 635     }
 636 }
 637
 638
 639 Foam::label Foam::meshRefinement::getShiftedRegion
 640 (
 641     const globalIndex& indexer,
 642     const Map<label>& globalToLocalRegion,
 643     const Map<label>& coupledRegionToShifted,
 644     const label globalRegion
 645 )
 646 {
 647     Map<label>::const_iterator iter =
 648         globalToLocalRegion.find(globalRegion);
 649
 650     if (iter != globalToLocalRegion.end())
 651     {
 652         // Region is 'owned' locally. Convert local region index into global.
 653         return indexer.toGlobal(iter());
 654     }
 655     else
 656     {
 657         return coupledRegionToShifted[globalRegion];
 658     }
 659 }
 660
 661
 662 // Add if not yet present
 663 void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
 664 {
 665     if (findIndex(lst, elem) == -1)
 666     {
 667         label sz = lst.size();
 668         lst.setSize(sz+1);
 669         lst[sz] = elem;
 670     }
 671 }
 672
 673
 674 void Foam::meshRefinement::calcRegionRegions
 675 (
 676     const labelList& globalRegion,
 677     const Map<label>& globalToLocalRegion,
 678     const Map<label>& coupledRegionToMaster,
 679     labelListList& regionRegions
 680 ) const
 681 {
 682     // Global region indexing since we now know the shifted regions.
 683     globalIndex shiftIndexer(globalToLocalRegion.size());
 684
 685     // Redo the coupledRegionToMaster to be in shifted region indexing.
 686     Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
 687     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
 688     {
 689         label region = iter.key();
 690
 691         Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
 692
 693         if (fndRegion != globalToLocalRegion.end())
 694         {
 695             // A local cell is master of this region. Get its shifted region.
 696             coupledRegionToShifted.insert
 697             (
 698                 region,
 699                 shiftIndexer.toGlobal(fndRegion())
 700             );
 701         }
 702         Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
 703         Pstream::mapCombineScatter(coupledRegionToShifted);
 704     }
 705
 706
 707     // Determine region-region connectivity.
 708     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 709     // This is for all my regions (so my local ones or the ones I am master
 710     // of) the neighbouring region indices.
 711
 712
 713     // Transfer lists.
 714     PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
 715     forAll(regionConnectivity, procI)
 716     {
 717         if (procI != Pstream::myProcNo())
 718         {
 719             regionConnectivity.set
 720             (
 721                 procI,
 722                 new HashSet<edge, Hash<edge> >
 723                 (
 724                     coupledRegionToShifted.size()
 725                   / Pstream::nProcs()
 726                 )
 727             );
 728         }
 729     }
 730
 731
 732     // Connectivity. For all my local regions the connected regions.
 733     regionRegions.setSize(globalToLocalRegion.size());
 734
 735     // Add all local connectivity to regionRegions; add all non-local
 736     // to the transferlists.
 737     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 738     {
 739         label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
 740         label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
 741
 742         if (ownRegion != neiRegion)
 743         {
 744             label shiftOwnRegion = getShiftedRegion
 745             (
 746                 shiftIndexer,
 747                 globalToLocalRegion,
 748                 coupledRegionToShifted,
 749                 ownRegion
 750             );
 751             label shiftNeiRegion = getShiftedRegion
 752             (
 753                 shiftIndexer,
 754                 globalToLocalRegion,
 755                 coupledRegionToShifted,
 756                 neiRegion
 757             );
 758
 759
 760             // Connection between two regions. Send to owner of region.
 761             // - is ownRegion 'owned' by me
 762             // - is neiRegion 'owned' by me
 763
 764             if (shiftIndexer.isLocal(shiftOwnRegion))
 765             {
 766                 label localI = shiftIndexer.toLocal(shiftOwnRegion);
 767                 addUnique(shiftNeiRegion, regionRegions[localI]);
 768             }
 769             else
 770             {
 771                 label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
 772                 regionConnectivity[masterProc].insert
 773                 (
 774                     edge(shiftOwnRegion, shiftNeiRegion)
 775                 );
 776             }
 777
 778             if (shiftIndexer.isLocal(shiftNeiRegion))
 779             {
 780                 label localI = shiftIndexer.toLocal(shiftNeiRegion);
 781                 addUnique(shiftOwnRegion, regionRegions[localI]);
 782             }
 783             else
 784             {
 785                 label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
 786                 regionConnectivity[masterProc].insert
 787                 (
 788                     edge(shiftOwnRegion, shiftNeiRegion)
 789                 );
 790             }
 791         }
 792     }
 793
 794
 795     // Send
 796     forAll(regionConnectivity, procI)
 797     {
 798         if (procI != Pstream::myProcNo())
 799         {
 800             OPstream str(Pstream::blocking, procI);
 801             str << regionConnectivity[procI];
 802         }
 803     }
 804     // Receive
 805     forAll(regionConnectivity, procI)
 806     {
 807         if (procI != Pstream::myProcNo())
 808         {
 809             IPstream str(Pstream::blocking, procI);
 810             str >> regionConnectivity[procI];
 811         }
 812     }
 813
 814     // Add to addressing.
 815     forAll(regionConnectivity, procI)
 816     {
 817         if (procI != Pstream::myProcNo())
 818         {
 819             for
 820             (
 821                 HashSet<edge, Hash<edge> >::const_iterator iter =
 822                     regionConnectivity[procI].begin();
 823                 iter != regionConnectivity[procI].end();
 824                 ++iter
 825             )
 826             {
 827                 const edge& e = iter.key();
 828
 829                 bool someLocal = false;
 830                 if (shiftIndexer.isLocal(e[0]))
 831                 {
 832                     label localI = shiftIndexer.toLocal(e[0]);
 833                     addUnique(e[1], regionRegions[localI]);
 834                     someLocal = true;
 835                 }
 836                 if (shiftIndexer.isLocal(e[1]))
 837                 {
 838                     label localI = shiftIndexer.toLocal(e[1]);
 839                     addUnique(e[0], regionRegions[localI]);
 840                     someLocal = true;
 841                 }
 842
 843                 if (!someLocal)
 844                 {
 845                     FatalErrorIn("calcRegionRegions(..)")
 846                         << "Received from processor " << procI
 847                         << " connection " << e
 848                         << " where none of the elements is local to me."
 849                         << abort(FatalError);
 850                 }
 851             }
 852         }
 853     }
 854 }
 855
 856
 857 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 858
 859 // Construct from components
 860 Foam::meshRefinement::meshRefinement
 861 (
 862     fvMesh& mesh,
 863     const scalar mergeDistance,
 864     const bool overwrite,
 865     const refinementSurfaces& surfaces,
 866     const refinementFeatures& features,
 867     const shellSurfaces& shells
 868 )
 869 :
 870     mesh_(mesh),
 871     mergeDistance_(mergeDistance),
 872     overwrite_(overwrite),
 873     oldInstance_(mesh.pointsInstance()),
 874     surfaces_(surfaces),
 875     features_(features),
 876     shells_(shells),
 877     meshCutter_
 878     (
 879         mesh,
 880         labelIOList
 881         (
 882             IOobject
 883             (
 884                 "cellLevel",
 885                 mesh_.facesInstance(),
 886                 fvMesh::meshSubDir,
 887                 mesh,
 888                 IOobject::READ_IF_PRESENT,
 889                 IOobject::NO_WRITE,
 890                 false
 891             ),
 892             labelList(mesh_.nCells(), 0)
 893         ),
 894         labelIOList
 895         (
 896             IOobject
 897             (
 898                 "pointLevel",
 899                 mesh_.facesInstance(),
 900                 fvMesh::meshSubDir,
 901                 mesh,
 902                 IOobject::READ_IF_PRESENT,
 903                 IOobject::NO_WRITE,
 904                 false
 905             ),
 906             labelList(mesh_.nPoints(), 0)
 907         ),
 908         refinementHistory
 909         (
 910             IOobject
 911             (
 912                 "refinementHistory",
 913                 mesh_.facesInstance(),
 914                 fvMesh::meshSubDir,
 915                 mesh_,
 916                 IOobject::NO_READ,
 917                 IOobject::NO_WRITE,
 918                 false
 919             ),
 920             List<refinementHistory::splitCell8>(0),
 921             labelList(0)
 922         )   // no unrefinement
 923     ),
 924     surfaceIndex_
 925     (
 926         IOobject
 927         (
 928             "surfaceIndex",
 929             mesh_.facesInstance(),
 930             fvMesh::meshSubDir,
 931             mesh_,
 932             IOobject::NO_READ,
 933             IOobject::NO_WRITE,
 934             false
 935         ),
 936         labelList(mesh_.nFaces(), -1)
 937     ),
 938     userFaceData_(0)
 939 {
 940     // recalculate intersections for all faces
 941     updateIntersections(identity(mesh_.nFaces()));
 942 }
 943
 944
 945 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 946
 947 Foam::label Foam::meshRefinement::countHits() const
 948 {
 949     // Stats on edges to test. Count proc faces only once.
 950     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 951
 952     label nHits = 0;
 953
 954     forAll(surfaceIndex_, faceI)
 955     {
 956         if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
 957         {
 958             nHits++;
 959         }
 960     }
 961     return nHits;
 962 }
 963
 964
 965 // Determine distribution to move connected regions onto one processor.
 966 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 967 (
 968     const scalarField& cellWeights,
 969     const boolList& blockedFace,
 970     const List<labelPair>& explicitConnections,
 971     decompositionMethod& decomposer
 972 ) const
 973 {
 974     // Determine global regions, separated by blockedFaces
 975     regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
 976
 977     // Now globalRegion has global region per cell. Problem is that
 978     // the region might span multiple domains so we want to get
 979     // a "region master" per domain. Note that multi-processor
 980     // regions can only occur on cells on coupled patches.
 981     // Note: since the number of regions does not change by this the
 982     // process can be seen as just a shift of a region onto a single
 983     // processor.
 984
 985
 986     // Global cell numbering engine
 987     globalIndex globalCells(mesh_.nCells());
 988
 989
 990     // Determine per coupled region the master cell (lowest numbered cell
 991     // on lowest numbered processor)
 992     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 993     // (does not determine master for non-coupled=fully-local regions)
 994
 995     Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
 996     getCoupledRegionMaster
 997     (
 998         globalCells,
 999         blockedFace,
1000         globalRegion,
1001         coupledRegionToMaster
1002     );
1003
1004     // Determine my regions
1005     // ~~~~~~~~~~~~~~~~~~~~
1006     // These are the fully local ones or the coupled ones of which I am master.
1007
1008     Map<label> globalToLocalRegion;
1009     pointField localPoints;
1010     scalarField localWeights;
1011     calcLocalRegions
1012     (
1013         globalCells,
1014         globalRegion,
1015         coupledRegionToMaster,
1016         cellWeights,
1017
1018         globalToLocalRegion,
1019         localPoints,
1020         localWeights
1021     );
1022
1023
1024
1025     // Find distribution for regions
1026     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1027
1028     labelList regionDistribution;
1029
1030     if (isA<geomDecomp>(decomposer))
1031     {
1032         regionDistribution = decomposer.decompose(localPoints, localWeights);
1033     }
1034     else
1035     {
1036         labelListList regionRegions;
1037         calcRegionRegions
1038         (
1039             globalRegion,
1040             globalToLocalRegion,
1041             coupledRegionToMaster,
1042
1043             regionRegions
1044         );
1045
1046         regionDistribution = decomposer.decompose
1047         (
1048             regionRegions,
1049             localPoints,
1050             localWeights
1051         );
1052     }
1053
1054
1055
1056     // Convert region-based decomposition back to cell-based one
1057     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1058
1059     // Transfer destination processor back to all. Use global reduce for now.
1060     Map<label> regionToDist(coupledRegionToMaster.size());
1061     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
1062     {
1063         label region = iter.key();
1064
1065         Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
1066
1067         if (regionFnd != globalToLocalRegion.end())
1068         {
1069             // Master cell is local. Store my distribution.
1070             regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
1071         }
1072         else
1073         {
1074             // Master cell is not on this processor. Make sure it is overridden.
1075             regionToDist.insert(iter.key(), labelMax);
1076         }
1077     }
1078     Pstream::mapCombineGather(regionToDist, minEqOp<label>());
1079     Pstream::mapCombineScatter(regionToDist);
1080
1081
1082     // Determine destination for all cells
1083     labelList distribution(mesh_.nCells());
1084
1085     forAll(globalRegion, cellI)
1086     {
1087         Map<label>::const_iterator fndRegion =
1088             regionToDist.find(globalRegion[cellI]);
1089
1090         if (fndRegion != regionToDist.end())
1091         {
1092             distribution[cellI] = fndRegion();
1093         }
1094         else
1095         {
1096             // region is local to the processor.
1097             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
1098
1099             distribution[cellI] = regionDistribution[localRegionI];
1100         }
1101     }
1102
1103     return distribution;
1104 }
1105
1106
1107 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
1108 (
1109     const bool keepZoneFaces,
1110     const bool keepBaffles,
1111     const scalarField& cellWeights,
1112     decompositionMethod& decomposer,
1113     fvMeshDistribute& distributor
1114 )
1115 {
1116     autoPtr<mapDistributePolyMesh> map;
1117
1118     if (Pstream::parRun())
1119     {
1120         //if (debug_)
1121         //{
1122         //    const_cast<Time&>(mesh_.time())++;
1123         //}
1124
1125         // Wanted distribution
1126         labelList distribution;
1127
1128         if (keepZoneFaces || keepBaffles)
1129         {
1130             // Faces where owner and neighbour are not 'connected' so can
1131             // go to different processors.
1132             boolList blockedFace(mesh_.nFaces(), true);
1133             label nUnblocked = 0;
1134             // Pairs of baffles
1135             List<labelPair> couples;
1136
1137             if (keepZoneFaces)
1138             {
1139                 // Determine decomposition to keep/move surface zones
1140                 // on one processor. The reason is that snapping will make these
1141                 // into baffles, move and convert them back so if they were
1142                 // proc boundaries after baffling&moving the points might be no
1143                 // longer snychronised so recoupling will fail. To prevent this
1144                 // keep owner&neighbour of such a surface zone on the same
1145                 // processor.
1146
1147                 const wordList& fzNames = surfaces().faceZoneNames();
1148                 const faceZoneMesh& fZones = mesh_.faceZones();
1149                 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
1150
1151                 // Get faces whose owner and neighbour should stay together,
1152                 // i.e. they are not 'blocked'.
1153
1154                 forAll(fzNames, surfI)
1155                 {
1156                     if (fzNames[surfI].size())
1157                     {
1158                         // Get zone
1159                         const faceZone& fZone = fZones[fzNames[surfI]];
1160
1161                         forAll(fZone, i)
1162                         {
1163                             label faceI = fZone[i];
1164                             if (blockedFace[faceI])
1165                             {
1166                                 if
1167                                 (
1168                                     mesh_.isInternalFace(faceI)
1169                                  || pbm[pbm.whichPatch(faceI)].coupled()
1170                                 )
1171                                 {
1172                                     blockedFace[faceI] = false;
1173                                     nUnblocked++;
1174                                 }
1175                             }
1176                         }
1177                     }
1178                 }
1179
1180
1181                 // If the faceZones are not synchronised the blockedFace
1182                 // might not be synchronised. If you are sure the faceZones
1183                 // are synchronised remove below check.
1184                 syncTools::syncFaceList
1185                 (
1186                     mesh_,
1187                     blockedFace,
1188                     andEqOp<bool>()     // combine operator
1189                 );
1190             }
1191             reduce(nUnblocked, sumOp<label>());
1192
1193             if (keepZoneFaces)
1194             {
1195                 Info<< "Found " << nUnblocked
1196                     << " zoned faces to keep together." << endl;
1197             }
1198
1199             if (keepBaffles)
1200             {
1201                 // Get boundary baffles that need to stay together.
1202                 couples = getDuplicateFaces   // all baffles
1203                 (
1204                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
1205                    +mesh_.nInternalFaces()
1206                 );
1207             }
1208             label nCouples = returnReduce(couples.size(), sumOp<label>());
1209
1210             if (keepBaffles)
1211             {
1212                 Info<< "Found " << nCouples << " baffles to keep together."
1213                     << endl;
1214             }
1215
1216             if (nUnblocked > 0 || nCouples > 0)
1217             {
1218                 Info<< "Applying special decomposition to keep baffles"
1219                     << " and zoned faces together." << endl;
1220
1221                 distribution = decomposeCombineRegions
1222                 (
1223                     cellWeights,
1224                     blockedFace,
1225                     couples,
1226                     decomposer
1227                 );
1228
1229                 labelList nProcCells(distributor.countCells(distribution));
1230                 Pstream::listCombineGather(nProcCells, plusEqOp<label>());
1231                 Pstream::listCombineScatter(nProcCells);
1232
1233                 Info<< "Calculated decomposition:" << endl;
1234                 forAll(nProcCells, procI)
1235                 {
1236                     Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
1237                 }
1238                 Info<< endl;
1239             }
1240             else
1241             {
1242                 // Normal decomposition
1243                 distribution = decomposer.decompose
1244                 (
1245                     mesh_,
1246                     mesh_.cellCentres(),
1247                     cellWeights
1248                 );
1249             }
1250         }
1251         else
1252         {
1253             // Normal decomposition
1254             distribution = decomposer.decompose
1255             (
1256                 mesh_,
1257                 mesh_.cellCentres(),
1258                 cellWeights
1259             );
1260         }
1261
1262         if (debug)
1263         {
1264             labelList nProcCells(distributor.countCells(distribution));
1265             Pout<< "Wanted distribution:" << nProcCells << endl;
1266
1267             Pstream::listCombineGather(nProcCells, plusEqOp<label>());
1268             Pstream::listCombineScatter(nProcCells);
1269
1270             Pout<< "Wanted resulting decomposition:" << endl;
1271             forAll(nProcCells, procI)
1272             {
1273                 Pout<< "    " << procI << '\t' << nProcCells[procI] << endl;
1274             }
1275             Pout<< endl;
1276         }
1277         // Do actual sending/receiving of mesh
1278         map = distributor.distribute(distribution);
1279
1280         // Update numbering of meshRefiner
1281         distribute(map);
1282
1283         // Set correct instance (for if overwrite)
1284         mesh_.setInstance(timeName());
1285         setInstance(mesh_.facesInstance());
1286     }
1287     return map;
1288 }
1289
1290
1291 // Helper function to get intersected faces
1292 Foam::labelList Foam::meshRefinement::intersectedFaces() const
1293 {
1294     label nBoundaryFaces = 0;
1295
1296     forAll(surfaceIndex_, faceI)
1297     {
1298         if (surfaceIndex_[faceI] != -1)
1299         {
1300             nBoundaryFaces++;
1301         }
1302     }
1303
1304     labelList surfaceFaces(nBoundaryFaces);
1305     nBoundaryFaces = 0;
1306
1307     forAll(surfaceIndex_, faceI)
1308     {
1309         if (surfaceIndex_[faceI] != -1)
1310         {
1311             surfaceFaces[nBoundaryFaces++] = faceI;
1312         }
1313     }
1314     return surfaceFaces;
1315 }
1316
1317
1318 // Helper function to get points used by faces
1319 Foam::labelList Foam::meshRefinement::intersectedPoints() const
1320 {
1321     const faceList& faces = mesh_.faces();
1322
1323     // Mark all points on faces that will become baffles
1324     PackedBoolList isBoundaryPoint(mesh_.nPoints());
1325     label nBoundaryPoints = 0;
1326
1327     forAll(surfaceIndex_, faceI)
1328     {
1329         if (surfaceIndex_[faceI] != -1)
1330         {
1331             const face& f = faces[faceI];
1332
1333             forAll(f, fp)
1334             {
1335                 if (isBoundaryPoint.set(f[fp], 1u))
1336                 {
1337                     nBoundaryPoints++;
1338                 }
1339             }
1340         }
1341     }
1342
1343     //// Insert all meshed patches.
1344     //labelList adaptPatchIDs(meshedPatches());
1345     //forAll(adaptPatchIDs, i)
1346     //{
1347     //    label patchI = adaptPatchIDs[i];
1348     //
1349     //    if (patchI != -1)
1350     //    {
1351     //        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
1352     //
1353     //        label faceI = pp.start();
1354     //
1355     //        forAll(pp, i)
1356     //        {
1357     //            const face& f = faces[faceI];
1358     //
1359     //            forAll(f, fp)
1360     //            {
1361     //                if (isBoundaryPoint.set(f[fp], 1u))
1362     //                    nBoundaryPoints++;
1363     //                }
1364     //            }
1365     //            faceI++;
1366     //        }
1367     //    }
1368     //}
1369
1370
1371     // Pack
1372     labelList boundaryPoints(nBoundaryPoints);
1373     nBoundaryPoints = 0;
1374     forAll(isBoundaryPoint, pointI)
1375     {
1376         if (isBoundaryPoint.get(pointI) == 1u)
1377         {
1378             boundaryPoints[nBoundaryPoints++] = pointI;
1379         }
1380     }
1381
1382     return boundaryPoints;
1383 }
1384
1385
1386 //- Create patch from set of patches
1387 Foam::autoPtr<Foam::indirectPrimitivePatch> Foam::meshRefinement::makePatch
1388 (
1389     const polyMesh& mesh,
1390     const labelList& patchIDs
1391 )
1392 {
1393     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1394
1395     // Count faces.
1396     label nFaces = 0;
1397
1398     forAll(patchIDs, i)
1399     {
1400         const polyPatch& pp = patches[patchIDs[i]];
1401
1402         nFaces += pp.size();
1403     }
1404
1405     // Collect faces.
1406     labelList addressing(nFaces);
1407     nFaces = 0;
1408
1409     forAll(patchIDs, i)
1410     {
1411         const polyPatch& pp = patches[patchIDs[i]];
1412
1413         label meshFaceI = pp.start();
1414
1415         forAll(pp, i)
1416         {
1417             addressing[nFaces++] = meshFaceI++;
1418         }
1419     }
1420
1421     return autoPtr<indirectPrimitivePatch>
1422     (
1423         new indirectPrimitivePatch
1424         (
1425             IndirectList<face>(mesh.faces(), addressing),
1426             mesh.points()
1427         )
1428     );
1429 }
1430
1431
1432 // Construct pointVectorField with correct boundary conditions
1433 Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
1434 (
1435     const pointMesh& pMesh,
1436     const labelList& adaptPatchIDs
1437 )
1438 {
1439     const polyMesh& mesh = pMesh();
1440
1441     // Construct displacement field.
1442     const pointBoundaryMesh& pointPatches = pMesh.boundary();
1443
1444     wordList patchFieldTypes
1445     (
1446         pointPatches.size(),
1447         slipPointPatchVectorField::typeName
1448     );
1449
1450     forAll(adaptPatchIDs, i)
1451     {
1452         patchFieldTypes[adaptPatchIDs[i]] =
1453             fixedValuePointPatchVectorField::typeName;
1454     }
1455
1456     forAll(pointPatches, patchI)
1457     {
1458         if (isA<processorPointPatch>(pointPatches[patchI]))
1459         {
1460             patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
1461         }
1462         else if (isA<cyclicPointPatch>(pointPatches[patchI]))
1463         {
1464             patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
1465         }
1466     }
1467
1468     // Note: time().timeName() instead of meshRefinement::timeName() since
1469     // postprocessable field.
1470     tmp<pointVectorField> tfld
1471     (
1472         new pointVectorField
1473         (
1474             IOobject
1475             (
1476                 "pointDisplacement",
1477                 mesh.time().timeName(),
1478                 mesh,
1479                 IOobject::NO_READ,
1480                 IOobject::AUTO_WRITE
1481             ),
1482             pMesh,
1483             dimensionedVector("displacement", dimLength, vector::zero),
1484             patchFieldTypes
1485         )
1486     );
1487     return tfld;
1488 }
1489
1490
1491 void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
1492 {
1493     const faceZoneMesh& fZones = mesh.faceZones();
1494
1495     // Check any zones are present anywhere and in same order
1496
1497     {
1498         List<wordList> zoneNames(Pstream::nProcs());
1499         zoneNames[Pstream::myProcNo()] = fZones.names();
1500         Pstream::gatherList(zoneNames);
1501         Pstream::scatterList(zoneNames);
1502         // All have same data now. Check.
1503         forAll(zoneNames, procI)
1504         {
1505             if (procI != Pstream::myProcNo())
1506             {
1507                 if (zoneNames[procI] != zoneNames[Pstream::myProcNo()])
1508                 {
1509                     FatalErrorIn
1510                     (
1511                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1512                     )   << "faceZones are not synchronised on processors." << nl
1513                         << "Processor " << procI << " has faceZones "
1514                         << zoneNames[procI] << nl
1515                         << "Processor " << Pstream::myProcNo()
1516                         << " has faceZones "
1517                         << zoneNames[Pstream::myProcNo()] << nl
1518                         << exit(FatalError);
1519                 }
1520             }
1521         }
1522     }
1523
1524     // Check that coupled faces are present on both sides.
1525
1526     labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
1527
1528     forAll(fZones, zoneI)
1529     {
1530         const faceZone& fZone = fZones[zoneI];
1531
1532         forAll(fZone, i)
1533         {
1534             label bFaceI = fZone[i]-mesh.nInternalFaces();
1535
1536             if (bFaceI >= 0)
1537             {
1538                 if (faceToZone[bFaceI] == -1)
1539                 {
1540                     faceToZone[bFaceI] = zoneI;
1541                 }
1542                 else if (faceToZone[bFaceI] == zoneI)
1543                 {
1544                     FatalErrorIn
1545                     (
1546                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1547                     )   << "Face " << fZone[i] << " in zone "
1548                         << fZone.name()
1549                         << " is twice in zone!"
1550                         << abort(FatalError);
1551                 }
1552                 else
1553                 {
1554                     FatalErrorIn
1555                     (
1556                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1557                     )   << "Face " << fZone[i] << " in zone "
1558                         << fZone.name()
1559                         << " is also in zone "
1560                         << fZones[faceToZone[bFaceI]].name()
1561                         << abort(FatalError);
1562                 }
1563             }
1564         }
1565     }
1566
1567     labelList neiFaceToZone(faceToZone);
1568     syncTools::swapBoundaryFaceList(mesh, neiFaceToZone);
1569
1570     forAll(faceToZone, i)
1571     {
1572         if (faceToZone[i] != neiFaceToZone[i])
1573         {
1574             FatalErrorIn
1575             (
1576                 "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1577             )   << "Face " << mesh.nInternalFaces()+i
1578                 << " is in zone " << faceToZone[i]
1579                 << ", its coupled face is in zone " << neiFaceToZone[i]
1580                 << abort(FatalError);
1581         }
1582     }
1583 }
1584
1585
1586 Foam::label Foam::meshRefinement::appendPatch
1587 (
1588     fvMesh& mesh,
1589     const label insertPatchI,
1590     const word& patchName,
1591     const dictionary& patchDict
1592 )
1593 {
1594     // Clear local fields and e.g. polyMesh parallelInfo.
1595     mesh.clearOut();
1596
1597     polyBoundaryMesh& polyPatches =
1598         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
1599     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
1600
1601     label patchI = polyPatches.size();
1602
1603     // Add polyPatch at the end
1604     polyPatches.setSize(patchI+1);
1605     polyPatches.set
1606     (
1607         patchI,
1608         polyPatch::New
1609         (
1610             patchName,
1611             patchDict,
1612             insertPatchI,
1613             polyPatches
1614         )
1615     );
1616     fvPatches.setSize(patchI+1);
1617     fvPatches.set
1618     (
1619         patchI,
1620         fvPatch::New
1621         (
1622             polyPatches[patchI],  // point to newly added polyPatch
1623             mesh.boundary()
1624         )
1625     );
1626
1627     addPatchFields<volScalarField>
1628     (
1629         mesh,
1630         calculatedFvPatchField<scalar>::typeName
1631     );
1632     addPatchFields<volVectorField>
1633     (
1634         mesh,
1635         calculatedFvPatchField<vector>::typeName
1636     );
1637     addPatchFields<volSphericalTensorField>
1638     (
1639         mesh,
1640         calculatedFvPatchField<sphericalTensor>::typeName
1641     );
1642     addPatchFields<volSymmTensorField>
1643     (
1644         mesh,
1645         calculatedFvPatchField<symmTensor>::typeName
1646     );
1647     addPatchFields<volTensorField>
1648     (
1649         mesh,
1650         calculatedFvPatchField<tensor>::typeName
1651     );
1652
1653     // Surface fields
1654
1655     addPatchFields<surfaceScalarField>
1656     (
1657         mesh,
1658         calculatedFvPatchField<scalar>::typeName
1659     );
1660     addPatchFields<surfaceVectorField>
1661     (
1662         mesh,
1663         calculatedFvPatchField<vector>::typeName
1664     );
1665     addPatchFields<surfaceSphericalTensorField>
1666     (
1667         mesh,
1668         calculatedFvPatchField<sphericalTensor>::typeName
1669     );
1670     addPatchFields<surfaceSymmTensorField>
1671     (
1672         mesh,
1673         calculatedFvPatchField<symmTensor>::typeName
1674     );
1675     addPatchFields<surfaceTensorField>
1676     (
1677         mesh,
1678         calculatedFvPatchField<tensor>::typeName
1679     );
1680     return patchI;
1681 }
1682
1683
1684 // Adds patch if not yet there. Returns patchID.
1685 Foam::label Foam::meshRefinement::addPatch
1686 (
1687     fvMesh& mesh,
1688     const word& patchName,
1689     const dictionary& patchInfo
1690 )
1691 {
1692     polyBoundaryMesh& polyPatches =
1693         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
1694     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
1695
1696     const label patchI = polyPatches.findPatchID(patchName);
1697     if (patchI != -1)
1698     {
1699         // Already there
1700         return patchI;
1701     }
1702
1703
1704     label insertPatchI = polyPatches.size();
1705     label startFaceI = mesh.nFaces();
1706
1707     forAll(polyPatches, patchI)
1708     {
1709         const polyPatch& pp = polyPatches[patchI];
1710
1711         if (isA<processorPolyPatch>(pp))
1712         {
1713             insertPatchI = patchI;
1714             startFaceI = pp.start();
1715             break;
1716         }
1717     }
1718
1719     dictionary patchDict(patchInfo);
1720     patchDict.set("nFaces", 0);
1721     patchDict.set("startFace", startFaceI);
1722
1723     // Below is all quite a hack. Feel free to change once there is a better
1724     // mechanism to insert and reorder patches.
1725
1726     label addedPatchI = appendPatch(mesh, insertPatchI, patchName, patchDict);
1727
1728
1729     // Create reordering list
1730     // patches before insert position stay as is
1731     labelList oldToNew(addedPatchI+1);
1732     for (label i = 0; i < insertPatchI; i++)
1733     {
1734         oldToNew[i] = i;
1735     }
1736     // patches after insert position move one up
1737     for (label i = insertPatchI; i < addedPatchI; i++)
1738     {
1739         oldToNew[i] = i+1;
1740     }
1741     // appended patch gets moved to insert position
1742     oldToNew[addedPatchI] = insertPatchI;
1743
1744     // Shuffle into place
1745     polyPatches.reorder(oldToNew);
1746     fvPatches.reorder(oldToNew);
1747
1748     reorderPatchFields<volScalarField>(mesh, oldToNew);
1749     reorderPatchFields<volVectorField>(mesh, oldToNew);
1750     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
1751     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
1752     reorderPatchFields<volTensorField>(mesh, oldToNew);
1753     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
1754     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
1755     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
1756     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
1757     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
1758
1759     return insertPatchI;
1760 }
1761
1762
1763 Foam::label Foam::meshRefinement::addMeshedPatch
1764 (
1765     const word& name,
1766     const dictionary& patchInfo
1767 )
1768 {
1769     label meshedI = findIndex(meshedPatches_, name);
1770
1771     if (meshedI != -1)
1772     {
1773         // Already there. Get corresponding polypatch
1774         return mesh_.boundaryMesh().findPatchID(name);
1775     }
1776     else
1777     {
1778         // Add patch
1779         label patchI = addPatch(mesh_, name, patchInfo);
1780
1781         // Store
1782         label sz = meshedPatches_.size();
1783         meshedPatches_.setSize(sz+1);
1784         meshedPatches_[sz] = name;
1785
1786         return patchI;
1787     }
1788 }
1789
1790
1791 Foam::labelList Foam::meshRefinement::meshedPatches() const
1792 {
1793     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1794
1795     DynamicList<label> patchIDs(meshedPatches_.size());
1796     forAll(meshedPatches_, i)
1797     {
1798         label patchI = patches.findPatchID(meshedPatches_[i]);
1799
1800         if (patchI == -1)
1801         {
1802             FatalErrorIn("meshRefinement::meshedPatches() const")
1803                 << "Problem : did not find patch " << meshedPatches_[i]
1804                 << endl << "Valid patches are " << patches.names()
1805                 << abort(FatalError);
1806         }
1807         if (!polyPatch::constraintType(patches[patchI].type()))
1808         {
1809             patchIDs.append(patchI);
1810         }
1811     }
1812
1813     return patchIDs;
1814 }
1815
1816
1817 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
1818 (
1819     const point& keepPoint
1820 )
1821 {
1822     // Determine connected regions. regionSplit is the labelList with the
1823     // region per cell.
1824     regionSplit cellRegion(mesh_);
1825
1826     label regionI = -1;
1827
1828     label cellI = mesh_.findCell(keepPoint);
1829
1830     if (cellI != -1)
1831     {
1832         regionI = cellRegion[cellI];
1833     }
1834
1835     reduce(regionI, maxOp<label>());
1836
1837     if (regionI == -1)
1838     {
1839         FatalErrorIn
1840         (
1841             "meshRefinement::splitMeshRegions(const point&)"
1842         )   << "Point " << keepPoint
1843             << " is not inside the mesh." << nl
1844             << "Bounding box of the mesh:" << mesh_.bounds()
1845             << exit(FatalError);
1846     }
1847
1848     // Subset
1849     // ~~~~~~
1850
1851     // Get cells to remove
1852     DynamicList<label> cellsToRemove(mesh_.nCells());
1853     forAll(cellRegion, cellI)
1854     {
1855         if (cellRegion[cellI] != regionI)
1856         {
1857             cellsToRemove.append(cellI);
1858         }
1859     }
1860     cellsToRemove.shrink();
1861
1862     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
1863     reduce(nCellsToKeep, sumOp<label>());
1864
1865     Info<< "Keeping all cells in region " << regionI
1866         << " containing point " << keepPoint << endl
1867         << "Selected for keeping : "
1868         << nCellsToKeep
1869         << " cells." << endl;
1870
1871
1872     // Remove cells
1873     removeCells cellRemover(mesh_);
1874
1875     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1876
1877     if (exposedFaces.size())
1878     {
1879         FatalErrorIn
1880         (
1881             "meshRefinement::splitMeshRegions(const point&)"
1882         )   << "Removing non-reachable cells should only expose boundary faces"
1883             << nl
1884             << "ExposedFaces:" << exposedFaces << abort(FatalError);
1885     }
1886
1887     return doRemoveCells
1888     (
1889         cellsToRemove,
1890         exposedFaces,
1891         labelList(exposedFaces.size(),-1),  // irrelevant since 0 size.
1892         cellRemover
1893     );
1894 }
1895
1896
1897 void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
1898 {
1899     // mesh_ already distributed; distribute my member data
1900
1901     // surfaceQueries_ ok.
1902
1903     // refinement
1904     meshCutter_.distribute(map);
1905
1906     // surfaceIndex is face data.
1907     map.distributeFaceData(surfaceIndex_);
1908
1909     // maintainedFaces are indices of faces.
1910     forAll(userFaceData_, i)
1911     {
1912         map.distributeFaceData(userFaceData_[i].second());
1913     }
1914
1915     // Redistribute surface and any fields on it.
1916     {
1917         Random rndGen(653213);
1918
1919         // Get local mesh bounding box. Single box for now.
1920         List<treeBoundBox> meshBb(1);
1921         treeBoundBox& bb = meshBb[0];
1922         bb = treeBoundBox(mesh_.points());
1923         bb = bb.extend(rndGen, 1E-4);
1924
1925         // Distribute all geometry (so refinementSurfaces and shellSurfaces)
1926         searchableSurfaces& geometry =
1927             const_cast<searchableSurfaces&>(surfaces_.geometry());
1928
1929         forAll(geometry, i)
1930         {
1931             autoPtr<mapDistribute> faceMap;
1932             autoPtr<mapDistribute> pointMap;
1933             geometry[i].distribute
1934             (
1935                 meshBb,
1936                 false,          // do not keep outside triangles
1937                 faceMap,
1938                 pointMap
1939             );
1940
1941             if (faceMap.valid())
1942             {
1943                 // (ab)use the instance() to signal current modification time
1944                 geometry[i].instance() = geometry[i].time().timeName();
1945             }
1946
1947             faceMap.clear();
1948             pointMap.clear();
1949         }
1950     }
1951 }
1952
1953
1954 // Update local data for a mesh change
1955 void Foam::meshRefinement::updateMesh
1956 (
1957     const mapPolyMesh& map,
1958     const labelList& changedFaces
1959 )
1960 {
1961     Map<label> dummyMap(0);
1962
1963     updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
1964 }
1965
1966
1967 void Foam::meshRefinement::storeData
1968 (
1969     const labelList& pointsToStore,
1970     const labelList& facesToStore,
1971     const labelList& cellsToStore
1972 )
1973 {
1974     // For now only meshCutter has storable/retrievable data.
1975     meshCutter_.storeData
1976     (
1977         pointsToStore,
1978         facesToStore,
1979         cellsToStore
1980     );
1981 }
1982
1983
1984 void Foam::meshRefinement::updateMesh
1985 (
1986     const mapPolyMesh& map,
1987     const labelList& changedFaces,
1988     const Map<label>& pointsToRestore,
1989     const Map<label>& facesToRestore,
1990     const Map<label>& cellsToRestore
1991 )
1992 {
1993     // For now only meshCutter has storable/retrievable data.
1994
1995     // Update numbering of cells/vertices.
1996     meshCutter_.updateMesh
1997     (
1998         map,
1999         pointsToRestore,
2000         facesToRestore,
2001         cellsToRestore
2002     );
2003
2004     // Update surfaceIndex
2005     updateList(map.faceMap(), -1, surfaceIndex_);
2006
2007     // Update cached intersection information
2008     updateIntersections(changedFaces);
2009
2010     // Update maintained faces
2011     forAll(userFaceData_, i)
2012     {
2013         labelList& data = userFaceData_[i].second();
2014
2015         if (userFaceData_[i].first() == KEEPALL)
2016         {
2017             // extend list with face-from-face data
2018             updateList(map.faceMap(), -1, data);
2019         }
2020         else if (userFaceData_[i].first() == MASTERONLY)
2021         {
2022             // keep master only
2023             labelList newFaceData(map.faceMap().size(), -1);
2024
2025             forAll(newFaceData, faceI)
2026             {
2027                 label oldFaceI = map.faceMap()[faceI];
2028
2029                 if (oldFaceI >= 0 && map.reverseFaceMap()[oldFaceI] == faceI)
2030                 {
2031                     newFaceData[faceI] = data[oldFaceI];
2032                 }
2033             }
2034             data.transfer(newFaceData);
2035         }
2036         else
2037         {
2038             // remove any face that has been refined i.e. referenced more than
2039             // once.
2040
2041             // 1. Determine all old faces that get referenced more than once.
2042             // These get marked with -1 in reverseFaceMap
2043             labelList reverseFaceMap(map.reverseFaceMap());
2044
2045             forAll(map.faceMap(), faceI)
2046             {
2047                 label oldFaceI = map.faceMap()[faceI];
2048
2049                 if (oldFaceI >= 0)
2050                 {
2051                     if (reverseFaceMap[oldFaceI] != faceI)
2052                     {
2053                         // faceI is slave face. Mark old face.
2054                         reverseFaceMap[oldFaceI] = -1;
2055                     }
2056                 }
2057             }
2058
2059             // 2. Map only faces with intact reverseFaceMap
2060             labelList newFaceData(map.faceMap().size(), -1);
2061             forAll(newFaceData, faceI)
2062             {
2063                 label oldFaceI = map.faceMap()[faceI];
2064
2065                 if (oldFaceI >= 0)
2066                 {
2067                     if (reverseFaceMap[oldFaceI] == faceI)
2068                     {
2069                         newFaceData[faceI] = data[oldFaceI];
2070                     }
2071                 }
2072             }
2073             data.transfer(newFaceData);
2074         }
2075     }
2076 }
2077
2078
2079 bool Foam::meshRefinement::write() const
2080 {
2081     bool writeOk =
2082         mesh_.write()
2083      && meshCutter_.write()
2084      && surfaceIndex_.write();
2085
2086
2087     // Make sure that any distributed surfaces (so ones which probably have
2088     // been changed) get written as well.
2089     // Note: should ideally have some 'modified' flag to say whether it
2090     // has been changed or not.
2091     searchableSurfaces& geometry =
2092         const_cast<searchableSurfaces&>(surfaces_.geometry());
2093
2094     forAll(geometry, i)
2095     {
2096         searchableSurface& s = geometry[i];
2097
2098         // Check if instance() of surface is not constant or system.
2099         // Is good hint that surface is distributed.
2100         if
2101         (
2102             s.instance() != s.time().system()
2103          && s.instance() != s.time().caseSystem()
2104          && s.instance() != s.time().constant()
2105          && s.instance() != s.time().caseConstant()
2106         )
2107         {
2108             // Make sure it gets written to current time, not constant.
2109             s.instance() = s.time().timeName();
2110             writeOk = writeOk && s.write();
2111         }
2112     }
2113
2114     return writeOk;
2115 }
2116
2117
2118 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
2119  const
2120 {
2121     const globalMeshData& pData = mesh_.globalData();
2122
2123     if (debug)
2124     {
2125         Pout<< msg.c_str()
2126             << " : cells(local):" << mesh_.nCells()
2127             << "  faces(local):" << mesh_.nFaces()
2128             << "  points(local):" << mesh_.nPoints()
2129             << endl;
2130     }
2131
2132     {
2133         PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
2134         label nMasterFaces = 0;
2135         forAll(isMasterFace, i)
2136         {
2137             if (isMasterFace[i])
2138             {
2139                 nMasterFaces++;
2140             }
2141         }
2142
2143         PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh_));
2144         label nMasterPoints = 0;
2145         forAll(isMasterPoint, i)
2146         {
2147             if (isMasterPoint[i])
2148             {
2149                 nMasterPoints++;
2150             }
2151         }
2152
2153         Info<< msg.c_str()
2154             << " : cells:" << pData.nTotalCells()
2155             << "  faces:" << returnReduce(nMasterFaces, sumOp<label>())
2156             << "  points:" << returnReduce(nMasterPoints, sumOp<label>())
2157             << endl;
2158     }
2159
2160
2161     //if (debug)
2162     {
2163         const labelList& cellLevel = meshCutter_.cellLevel();
2164
2165         labelList nCells(gMax(cellLevel)+1, 0);
2166
2167         forAll(cellLevel, cellI)
2168         {
2169             nCells[cellLevel[cellI]]++;
2170         }
2171
2172         Pstream::listCombineGather(nCells, plusEqOp<label>());
2173         Pstream::listCombineScatter(nCells);
2174
2175         Info<< "Cells per refinement level:" << endl;
2176         forAll(nCells, levelI)
2177         {
2178             Info<< "    " << levelI << '\t' << nCells[levelI]
2179                 << endl;
2180         }
2181     }
2182 }
2183
2184
2185 //- Return either time().constant() or oldInstance
2186 Foam::word Foam::meshRefinement::timeName() const
2187 {
2188     if (overwrite_ && mesh_.time().timeIndex() == 0)
2189     {
2190         return oldInstance_;
2191     }
2192     else
2193     {
2194         return mesh_.time().timeName();
2195     }
2196 }
2197
2198
2199 void Foam::meshRefinement::dumpRefinementLevel() const
2200 {
2201     // Note: use time().timeName(), not meshRefinement::timeName()
2202     // so as to dump the fields to 0, not to constant.
2203     volScalarField volRefLevel
2204     (
2205         IOobject
2206         (
2207             "cellLevel",
2208             mesh_.time().timeName(),
2209             mesh_,
2210             IOobject::NO_READ,
2211             IOobject::AUTO_WRITE,
2212             false
2213         ),
2214         mesh_,
2215         dimensionedScalar("zero", dimless, 0),
2216         zeroGradientFvPatchScalarField::typeName
2217     );
2218
2219     const labelList& cellLevel = meshCutter_.cellLevel();
2220
2221     forAll(volRefLevel, cellI)
2222     {
2223         volRefLevel[cellI] = cellLevel[cellI];
2224     }
2225
2226     volRefLevel.write();
2227
2228
2229     const pointMesh& pMesh = pointMesh::New(mesh_);
2230
2231     pointScalarField pointRefLevel
2232     (
2233         IOobject
2234         (
2235             "pointLevel",
2236             mesh_.time().timeName(),
2237             mesh_,
2238             IOobject::NO_READ,
2239             IOobject::NO_WRITE,
2240             false
2241         ),
2242         pMesh,
2243         dimensionedScalar("zero", dimless, 0)
2244     );
2245
2246     const labelList& pointLevel = meshCutter_.pointLevel();
2247
2248     forAll(pointRefLevel, pointI)
2249     {
2250         pointRefLevel[pointI] = pointLevel[pointI];
2251     }
2252
2253     pointRefLevel.write();
2254 }
2255
2256
2257 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
2258 {
2259     {
2260         const pointField& cellCentres = mesh_.cellCentres();
2261
2262         OFstream str(prefix + "_edges.obj");
2263         label vertI = 0;
2264         Pout<< "meshRefinement::dumpIntersections :"
2265             << " Writing cellcentre-cellcentre intersections to file "
2266             << str.name() << endl;
2267
2268
2269         // Redo all intersections
2270         // ~~~~~~~~~~~~~~~~~~~~~~
2271
2272         // Get boundary face centre and level. Coupled aware.
2273         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2274         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2275         calcNeighbourData(neiLevel, neiCc);
2276
2277         labelList intersectionFaces(intersectedFaces());
2278
2279         // Collect segments we want to test for
2280         pointField start(intersectionFaces.size());
2281         pointField end(intersectionFaces.size());
2282
2283         forAll(intersectionFaces, i)
2284         {
2285             label faceI = intersectionFaces[i];
2286             start[i] = cellCentres[mesh_.faceOwner()[faceI]];
2287
2288             if (mesh_.isInternalFace(faceI))
2289             {
2290                 end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
2291             }
2292             else
2293             {
2294                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
2295             }
2296         }
2297
2298         // Extend segments a bit
2299         {
2300             const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
2301             start -= smallVec;
2302             end += smallVec;
2303         }
2304
2305
2306         // Do tests in one go
2307         labelList surfaceHit;
2308         List<pointIndexHit> surfaceHitInfo;
2309         surfaces_.findAnyIntersection
2310         (
2311             start,
2312             end,
2313             surfaceHit,
2314             surfaceHitInfo
2315         );
2316
2317         forAll(intersectionFaces, i)
2318         {
2319             if (surfaceHit[i] != -1)
2320             {
2321                 meshTools::writeOBJ(str, start[i]);
2322                 vertI++;
2323                 meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
2324                 vertI++;
2325                 meshTools::writeOBJ(str, end[i]);
2326                 vertI++;
2327                 str << "l " << vertI-2 << ' ' << vertI-1 << nl
2328                     << "l " << vertI-1 << ' ' << vertI << nl;
2329             }
2330         }
2331     }
2332
2333     // Convert to vtk format
2334     string cmd
2335     (
2336         "objToVTK " + prefix + "_edges.obj " + prefix + "_edges.vtk > /dev/null"
2337     );
2338     system(cmd.c_str());
2339
2340     Pout<< endl;
2341 }
2342
2343
2344 void Foam::meshRefinement::write
2345 (
2346     const label flag,
2347     const fileName& prefix
2348 ) const
2349 {
2350     if (flag & MESH)
2351     {
2352         write();
2353     }
2354     if (flag & SCALARLEVELS)
2355     {
2356         dumpRefinementLevel();
2357     }
2358     if (flag & OBJINTERSECTIONS && prefix.size())
2359     {
2360         dumpIntersections(prefix);
2361     }
2362 }
2363
2364
2365 // ************************************************************************* //