FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

meshRefinement.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*----------------------------------------------------------------------------*/
00025 
00026 #include "meshRefinement.H"
00027 #include <finiteVolume/volMesh.H>
00028 #include <finiteVolume/volFields.H>
00029 #include <finiteVolume/surfaceMesh.H>
00030 #include <OpenFOAM/syncTools.H>
00031 #include <OpenFOAM/Time.H>
00032 #include <dynamicMesh/refinementHistory.H>
00033 #include <autoMesh/refinementSurfaces.H>
00034 #include <decompositionMethods/decompositionMethod.H>
00035 #include <meshTools/regionSplit.H>
00036 #include <dynamicMesh/fvMeshDistribute.H>
00037 #include <OpenFOAM/indirectPrimitivePatch.H>
00038 #include <dynamicMesh/polyTopoChange.H>
00039 #include <dynamicMesh/removeCells.H>
00040 #include <OpenFOAM/mapDistributePolyMesh.H>
00041 #include <dynamicMesh/localPointRegion.H>
00042 #include <OpenFOAM/pointMesh.H>
00043 #include <OpenFOAM/pointFields.H>
00044 #include <OpenFOAM/slipPointPatchFields.H>
00045 #include <OpenFOAM/fixedValuePointPatchFields.H>
00046 #include <OpenFOAM/globalPointPatchFields.H>
00047 #include <OpenFOAM/calculatedPointPatchFields.H>
00048 #include <OpenFOAM/processorPointPatch.H>
00049 #include <OpenFOAM/globalIndex.H>
00050 #include <meshTools/meshTools.H>
00051 #include <OpenFOAM/OFstream.H>
00052 #include <decompositionMethods/geomDecomp.H>
00053 #include <OpenFOAM/Random.H>
00054 #include <meshTools/searchableSurfaces.H>
00055 #include <meshTools/treeBoundBox.H>
00056 #include <finiteVolume/zeroGradientFvPatchFields.H>
00057 
00058 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00059 
00060 namespace Foam
00061 {
00062     defineTypeNameAndDebug(meshRefinement, 0);
00063 }
00064 
00065 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00066 
00067 void Foam::meshRefinement::calcNeighbourData
00068 (
00069     labelList& neiLevel,
00070     pointField& neiCc
00071 )  const
00072 {
00073     const labelList& cellLevel = meshCutter_.cellLevel();
00074     const pointField& cellCentres = mesh_.cellCentres();
00075 
00076     label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
00077 
00078     if (neiLevel.size() != nBoundaryFaces || neiCc.size() != nBoundaryFaces)
00079     {
00080         FatalErrorIn("meshRefinement::calcNeighbour(..)") << "nBoundaries:"
00081             << nBoundaryFaces << " neiLevel:" << neiLevel.size()
00082             << abort(FatalError);
00083     }
00084 
00085     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00086 
00087     labelHashSet addedPatchIDSet = meshedPatches();
00088 
00089     forAll(patches, patchI)
00090     {
00091         const polyPatch& pp = patches[patchI];
00092 
00093         const unallocLabelList& faceCells = pp.faceCells();
00094         const vectorField::subField faceCentres = pp.faceCentres();
00095         const vectorField::subField faceAreas = pp.faceAreas();
00096 
00097         label bFaceI = pp.start()-mesh_.nInternalFaces();
00098 
00099         if (pp.coupled())
00100         {
00101             forAll(faceCells, i)
00102             {
00103                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
00104                 neiCc[bFaceI] = cellCentres[faceCells[i]];
00105                 bFaceI++;
00106             }
00107         }
00108         else if (addedPatchIDSet.found(patchI))
00109         {
00110             // Face was introduced from cell-cell intersection. Try to
00111             // reconstruct other side cell(centre). Three possibilities:
00112             // - cells same size.
00113             // - preserved cell smaller. Not handled.
00114             // - preserved cell larger.
00115             forAll(faceCells, i)
00116             {
00117                 // Extrapolate the face centre.
00118                 vector fn = faceAreas[i];
00119                 fn /= mag(fn)+VSMALL;
00120 
00121                 label own = faceCells[i];
00122                 label ownLevel = cellLevel[own];
00123                 label faceLevel = meshCutter_.getAnchorLevel(pp.start()+i);
00124 
00125                 // Normal distance from face centre to cell centre
00126                 scalar d = ((faceCentres[i] - cellCentres[own]) & fn);
00127                 if (faceLevel > ownLevel)
00128                 {
00129                     // Other cell more refined. Adjust normal distance
00130                     d *= 0.5;
00131                 }
00132                 neiLevel[bFaceI] = cellLevel[ownLevel];
00133                 // Calculate other cell centre by extrapolation
00134                 neiCc[bFaceI] = faceCentres[i] + d*fn;
00135                 bFaceI++;
00136             }
00137         }
00138         else
00139         {
00140             forAll(faceCells, i)
00141             {
00142                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
00143                 neiCc[bFaceI] = faceCentres[i];
00144                 bFaceI++;
00145             }
00146         }
00147     }
00148 
00149     // Swap coupled boundaries. Apply separation to cc since is coordinate.
00150     syncTools::swapBoundaryFaceList(mesh_, neiCc, true);
00151     syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
00152 }
00153 
00154 
00155 // Find intersections of edges (given by their two endpoints) with surfaces.
00156 // Returns first intersection if there are more than one.
00157 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
00158 {
00159     const pointField& cellCentres = mesh_.cellCentres();
00160 
00161     // Stats on edges to test. Count proc faces only once.
00162     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
00163 
00164     {
00165         label nMasterFaces = 0;
00166         forAll(isMasterFace, faceI)
00167         {
00168             if (isMasterFace.get(faceI) == 1)
00169             {
00170                 nMasterFaces++;
00171             }
00172         }
00173         reduce(nMasterFaces, sumOp<label>());
00174 
00175         label nChangedFaces = 0;
00176         forAll(changedFaces, i)
00177         {
00178             if (isMasterFace.get(changedFaces[i]) == 1)
00179             {
00180                 nChangedFaces++;
00181             }
00182         }
00183         reduce(nChangedFaces, sumOp<label>());
00184 
00185         Info<< "Edge intersection testing:" << nl
00186             << "    Number of edges             : " << nMasterFaces << nl
00187             << "    Number of edges to retest   : " << nChangedFaces
00188             << endl;
00189     }
00190 
00191 
00192     // Get boundary face centre and level. Coupled aware.
00193     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
00194     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
00195     calcNeighbourData(neiLevel, neiCc);
00196 
00197     // Collect segments we want to test for
00198     pointField start(changedFaces.size());
00199     pointField end(changedFaces.size());
00200 
00201     forAll(changedFaces, i)
00202     {
00203         label faceI = changedFaces[i];
00204         label own = mesh_.faceOwner()[faceI];
00205 
00206         start[i] = cellCentres[own];
00207         if (mesh_.isInternalFace(faceI))
00208         {
00209             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
00210         }
00211         else
00212         {
00213             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
00214         }
00215     }
00216 
00217     // Do tests in one go
00218     labelList surfaceHit;
00219     {
00220         labelList surfaceLevel;
00221         surfaces_.findHigherIntersection
00222         (
00223             start,
00224             end,
00225             labelList(start.size(), -1),    // accept any intersection
00226             surfaceHit,
00227             surfaceLevel
00228         );
00229     }
00230 
00231     // Keep just surface hit
00232     forAll(surfaceHit, i)
00233     {
00234         surfaceIndex_[changedFaces[i]] = surfaceHit[i];
00235     }
00236 
00237     // Make sure both sides have same information. This should be
00238     // case in general since same vectors but just to make sure.
00239     syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>(), false);
00240 
00241     label nHits = countHits();
00242     label nTotHits = returnReduce(nHits, sumOp<label>());
00243 
00244     Info<< "    Number of intersected edges : " << nTotHits << endl;
00245 
00246     // Set files to same time as mesh
00247     setInstance(mesh_.facesInstance());
00248 }
00249 
00250 
00251 void Foam::meshRefinement::checkData()
00252 {
00253     Pout<< "meshRefinement::checkData() : Checking refinement structure."
00254         << endl;
00255     meshCutter_.checkMesh();
00256 
00257     Pout<< "meshRefinement::checkData() : Checking refinement levels."
00258         << endl;
00259     meshCutter_.checkRefinementLevels(1, labelList(0));
00260 
00261 
00262     label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
00263 
00264     Pout<< "meshRefinement::checkData() : Checking synchronization."
00265         << endl;
00266 
00267     // Check face centres
00268     {
00269         // Boundary face centres
00270         pointField::subList boundaryFc
00271         (
00272             mesh_.faceCentres(),
00273             nBnd,
00274             mesh_.nInternalFaces()
00275         );
00276 
00277         // Get neighbouring face centres
00278         pointField neiBoundaryFc(boundaryFc);
00279         syncTools::swapBoundaryFaceList
00280         (
00281             mesh_,
00282             neiBoundaryFc,
00283             true
00284         );
00285 
00286         // Compare
00287         testSyncBoundaryFaceList
00288         (
00289             mergeDistance_,
00290             "testing faceCentres : ",
00291             boundaryFc,
00292             neiBoundaryFc
00293         );
00294     }
00295     // Check meshRefinement
00296     {
00297         // Get boundary face centre and level. Coupled aware.
00298         labelList neiLevel(nBnd);
00299         pointField neiCc(nBnd);
00300         calcNeighbourData(neiLevel, neiCc);
00301 
00302         // Collect segments we want to test for
00303         pointField start(mesh_.nFaces());
00304         pointField end(mesh_.nFaces());
00305 
00306         forAll(start, faceI)
00307         {
00308             start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
00309 
00310             if (mesh_.isInternalFace(faceI))
00311             {
00312                 end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
00313             }
00314             else
00315             {
00316                 end[faceI] = neiCc[faceI-mesh_.nInternalFaces()];
00317             }
00318         }
00319 
00320         // Do tests in one go
00321         labelList surfaceHit;
00322         {
00323             labelList surfaceLevel;
00324             surfaces_.findHigherIntersection
00325             (
00326                 start,
00327                 end,
00328                 labelList(start.size(), -1),    // accept any intersection
00329                 surfaceHit,
00330                 surfaceLevel
00331             );
00332         }
00333         // Get the coupled hit
00334         labelList neiHit
00335         (
00336             SubList<label>
00337             (
00338                 surfaceHit,
00339                 nBnd,
00340                 mesh_.nInternalFaces()
00341             )
00342         );
00343         syncTools::swapBoundaryFaceList(mesh_, neiHit, false);
00344 
00345         // Check
00346         forAll(surfaceHit, faceI)
00347         {
00348             if (surfaceIndex_[faceI] != surfaceHit[faceI])
00349             {
00350                 if (mesh_.isInternalFace(faceI))
00351                 {
00352                     WarningIn("meshRefinement::checkData()")
00353                         << "Internal face:" << faceI
00354                         << " fc:" << mesh_.faceCentres()[faceI]
00355                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
00356                         << " current:" << surfaceHit[faceI]
00357                         << " ownCc:"
00358                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
00359                         << " neiCc:"
00360                         << mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
00361                         << endl;
00362                 }
00363                 else if
00364                 (
00365                     surfaceIndex_[faceI]
00366                  != neiHit[faceI-mesh_.nInternalFaces()]
00367                 )
00368                 {
00369                     WarningIn("meshRefinement::checkData()")
00370                         << "Boundary face:" << faceI
00371                         << " fc:" << mesh_.faceCentres()[faceI]
00372                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
00373                         << " current:" << surfaceHit[faceI]
00374                         << " ownCc:"
00375                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
00376                         << " end:" << end[faceI]
00377                         << endl;
00378                 }
00379             }
00380         }
00381     }
00382     {
00383         labelList::subList boundarySurface
00384         (
00385             surfaceIndex_,
00386             mesh_.nFaces()-mesh_.nInternalFaces(),
00387             mesh_.nInternalFaces()
00388         );
00389 
00390         labelList neiBoundarySurface(boundarySurface);
00391         syncTools::swapBoundaryFaceList
00392         (
00393             mesh_,
00394             neiBoundarySurface,
00395             false
00396         );
00397 
00398         // Compare
00399         testSyncBoundaryFaceList
00400         (
00401             0,                              // tolerance
00402             "testing surfaceIndex() : ",
00403             boundarySurface,
00404             neiBoundarySurface
00405         );
00406     }
00407 
00408 
00409     // Find duplicate faces
00410     Pout<< "meshRefinement::checkData() : Counting duplicate faces."
00411         << endl;
00412 
00413     labelList duplicateFace
00414     (
00415         localPointRegion::findDuplicateFaces
00416         (
00417             mesh_,
00418             identity(mesh_.nFaces()-mesh_.nInternalFaces())
00419           + mesh_.nInternalFaces()
00420         )
00421     );
00422 
00423     // Count
00424     {
00425         label nDup = 0;
00426 
00427         forAll(duplicateFace, i)
00428         {
00429             if (duplicateFace[i] != -1)
00430             {
00431                 nDup++;
00432             }
00433         }
00434         nDup /= 2;  // will have counted both faces of duplicate
00435         Pout<< "meshRefinement::checkData() : Found " << nDup
00436             << " duplicate pairs of faces." << endl;
00437     }
00438 }
00439 
00440 
00441 void Foam::meshRefinement::setInstance(const fileName& inst)
00442 {
00443     meshCutter_.setInstance(inst);
00444     surfaceIndex_.instance() = inst;
00445 }
00446 
00447 
00448 // Remove cells. Put exposedFaces (output of getExposedFaces(cellsToRemove))
00449 // into exposedPatchIDs.
00450 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
00451 (
00452     const labelList& cellsToRemove,
00453     const labelList& exposedFaces,
00454     const labelList& exposedPatchIDs,
00455     removeCells& cellRemover
00456 )
00457 {
00458     polyTopoChange meshMod(mesh_);
00459 
00460     // Arbitrary: put exposed faces into last patch.
00461     cellRemover.setRefinement
00462     (
00463         cellsToRemove,
00464         exposedFaces,
00465         exposedPatchIDs,
00466         meshMod
00467     );
00468 
00469     // Change the mesh (no inflation)
00470     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00471 
00472     // Update fields
00473     mesh_.updateMesh(map);
00474 
00475     // Move mesh (since morphing might not do this)
00476     if (map().hasMotionPoints())
00477     {
00478         mesh_.movePoints(map().preMotionPoints());
00479     }
00480     else
00481     {
00482         // Delete mesh volumes. No other way to do this?
00483         mesh_.clearOut();
00484     }
00485 
00486     if (overwrite_)
00487     {
00488         mesh_.setInstance(oldInstance_);
00489     }
00490 
00491     // Update local mesh data
00492     cellRemover.updateMesh(map);
00493 
00494     // Update intersections. Recalculate intersections for exposed faces.
00495     labelList newExposedFaces = renumber
00496     (
00497         map().reverseFaceMap(),
00498         exposedFaces
00499     );
00500 
00501     //Pout<< "removeCells : updating intersections for "
00502     //    << newExposedFaces.size() << " newly exposed faces." << endl;
00503 
00504     updateMesh(map, newExposedFaces);
00505 
00506     return map;
00507 }
00508 
00509 
00510 // Determine for multi-processor regions the lowest numbered cell on the lowest
00511 // numbered processor.
00512 void Foam::meshRefinement::getCoupledRegionMaster
00513 (
00514     const globalIndex& globalCells,
00515     const boolList& blockedFace,
00516     const regionSplit& globalRegion,
00517     Map<label>& regionToMaster
00518 ) const
00519 {
00520     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00521 
00522     forAll(patches, patchI)
00523     {
00524         const polyPatch& pp = patches[patchI];
00525 
00526         if (isA<processorPolyPatch>(pp))
00527         {
00528             forAll(pp, i)
00529             {
00530                 label faceI = pp.start()+i;
00531 
00532                 if (!blockedFace[faceI])
00533                 {
00534                     // Only if there is a connection to the neighbour
00535                     // will there be a multi-domain region; if not through
00536                     // this face then through another.
00537 
00538                     label cellI = mesh_.faceOwner()[faceI];
00539                     label globalCellI = globalCells.toGlobal(cellI);
00540 
00541                     Map<label>::iterator iter =
00542                         regionToMaster.find(globalRegion[cellI]);
00543 
00544                     if (iter != regionToMaster.end())
00545                     {
00546                         label master = iter();
00547                         iter() = min(master, globalCellI);
00548                     }
00549                     else
00550                     {
00551                         regionToMaster.insert
00552                         (
00553                             globalRegion[cellI],
00554                             globalCellI
00555                         );
00556                     }
00557                 }
00558             }
00559         }
00560     }
00561 
00562     // Do reduction
00563     Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
00564     Pstream::mapCombineScatter(regionToMaster);
00565 }
00566 
00567 
00568 void Foam::meshRefinement::calcLocalRegions
00569 (
00570     const globalIndex& globalCells,
00571     const labelList& globalRegion,
00572     const Map<label>& coupledRegionToMaster,
00573     const scalarField& cellWeights,
00574 
00575     Map<label>& globalToLocalRegion,
00576     pointField& localPoints,
00577     scalarField& localWeights
00578 ) const
00579 {
00580     globalToLocalRegion.resize(globalRegion.size());
00581     DynamicList<point> localCc(globalRegion.size()/2);
00582     DynamicList<scalar> localWts(globalRegion.size()/2);
00583 
00584     forAll(globalRegion, cellI)
00585     {
00586         Map<label>::const_iterator fndMaster =
00587             coupledRegionToMaster.find(globalRegion[cellI]);
00588 
00589         if (fndMaster != coupledRegionToMaster.end())
00590         {
00591             // Multi-processor region.
00592             if (globalCells.toGlobal(cellI) == fndMaster())
00593             {
00594                 // I am master. Allocate region for me.
00595                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
00596                 localCc.append(mesh_.cellCentres()[cellI]);
00597                 localWts.append(cellWeights[cellI]);
00598             }
00599         }
00600         else
00601         {
00602             // Single processor region.
00603             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
00604             {
00605                 localCc.append(mesh_.cellCentres()[cellI]);
00606                 localWts.append(cellWeights[cellI]);
00607             }
00608         }
00609     }
00610 
00611     localPoints.transfer(localCc);
00612     localWeights.transfer(localWts);
00613 
00614     if (localPoints.size() != globalToLocalRegion.size())
00615     {
00616         FatalErrorIn("calcLocalRegions(..)")
00617             << "localPoints:" << localPoints.size()
00618             << " globalToLocalRegion:" << globalToLocalRegion.size()
00619             << exit(FatalError);
00620     }
00621 }
00622 
00623 
00624 Foam::label Foam::meshRefinement::getShiftedRegion
00625 (
00626     const globalIndex& indexer,
00627     const Map<label>& globalToLocalRegion,
00628     const Map<label>& coupledRegionToShifted,
00629     const label globalRegion
00630 )
00631 {
00632     Map<label>::const_iterator iter =
00633         globalToLocalRegion.find(globalRegion);
00634 
00635     if (iter != globalToLocalRegion.end())
00636     {
00637         // Region is 'owned' locally. Convert local region index into global.
00638         return indexer.toGlobal(iter());
00639     }
00640     else
00641     {
00642         return coupledRegionToShifted[globalRegion];
00643     }
00644 }
00645 
00646 
00647 // Add if not yet present
00648 void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
00649 {
00650     if (findIndex(lst, elem) == -1)
00651     {
00652         label sz = lst.size();
00653         lst.setSize(sz+1);
00654         lst[sz] = elem;
00655     }
00656 }
00657 
00658 
00659 void Foam::meshRefinement::calcRegionRegions
00660 (
00661     const labelList& globalRegion,
00662     const Map<label>& globalToLocalRegion,
00663     const Map<label>& coupledRegionToMaster,
00664     labelListList& regionRegions
00665 ) const
00666 {
00667     // Global region indexing since we now know the shifted regions.
00668     globalIndex shiftIndexer(globalToLocalRegion.size());
00669 
00670     // Redo the coupledRegionToMaster to be in shifted region indexing.
00671     Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
00672     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
00673     {
00674         label region = iter.key();
00675 
00676         Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
00677 
00678         if (fndRegion != globalToLocalRegion.end())
00679         {
00680             // A local cell is master of this region. Get its shifted region.
00681             coupledRegionToShifted.insert
00682             (
00683                 region,
00684                 shiftIndexer.toGlobal(fndRegion())
00685             );
00686         }
00687         Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
00688         Pstream::mapCombineScatter(coupledRegionToShifted);
00689     }
00690 
00691 
00692     // Determine region-region connectivity.
00693     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00694     // This is for all my regions (so my local ones or the ones I am master
00695     // of) the neighbouring region indices.
00696 
00697 
00698     // Transfer lists.
00699     PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
00700     forAll(regionConnectivity, procI)
00701     {
00702         if (procI != Pstream::myProcNo())
00703         {
00704             regionConnectivity.set
00705             (
00706                 procI,
00707                 new HashSet<edge, Hash<edge> >
00708                 (
00709                     coupledRegionToShifted.size()
00710                   / Pstream::nProcs()
00711                 )
00712             );
00713         }
00714     }
00715 
00716 
00717     // Connectivity. For all my local regions the connected regions.
00718     regionRegions.setSize(globalToLocalRegion.size());
00719 
00720     // Add all local connectivity to regionRegions; add all non-local
00721     // to the transferlists.
00722     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
00723     {
00724         label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
00725         label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
00726 
00727         if (ownRegion != neiRegion)
00728         {
00729             label shiftOwnRegion = getShiftedRegion
00730             (
00731                 shiftIndexer,
00732                 globalToLocalRegion,
00733                 coupledRegionToShifted,
00734                 ownRegion
00735             );
00736             label shiftNeiRegion = getShiftedRegion
00737             (
00738                 shiftIndexer,
00739                 globalToLocalRegion,
00740                 coupledRegionToShifted,
00741                 neiRegion
00742             );
00743 
00744 
00745             // Connection between two regions. Send to owner of region.
00746             // - is ownRegion 'owned' by me
00747             // - is neiRegion 'owned' by me
00748 
00749             if (shiftIndexer.isLocal(shiftOwnRegion))
00750             {
00751                 label localI = shiftIndexer.toLocal(shiftOwnRegion);
00752                 addUnique(shiftNeiRegion, regionRegions[localI]);
00753             }
00754             else
00755             {
00756                 label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
00757                 regionConnectivity[masterProc].insert
00758                 (
00759                     edge(shiftOwnRegion, shiftNeiRegion)
00760                 );
00761             }
00762 
00763             if (shiftIndexer.isLocal(shiftNeiRegion))
00764             {
00765                 label localI = shiftIndexer.toLocal(shiftNeiRegion);
00766                 addUnique(shiftOwnRegion, regionRegions[localI]);
00767             }
00768             else
00769             {
00770                 label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
00771                 regionConnectivity[masterProc].insert
00772                 (
00773                     edge(shiftOwnRegion, shiftNeiRegion)
00774                 );
00775             }
00776         }
00777     }
00778 
00779 
00780     // Send
00781     forAll(regionConnectivity, procI)
00782     {
00783         if (procI != Pstream::myProcNo())
00784         {
00785             OPstream str(Pstream::blocking, procI);
00786             str << regionConnectivity[procI];
00787         }
00788     }
00789     // Receive
00790     forAll(regionConnectivity, procI)
00791     {
00792         if (procI != Pstream::myProcNo())
00793         {
00794             IPstream str(Pstream::blocking, procI);
00795             str >> regionConnectivity[procI];
00796         }
00797     }
00798 
00799     // Add to addressing.
00800     forAll(regionConnectivity, procI)
00801     {
00802         if (procI != Pstream::myProcNo())
00803         {
00804             for
00805             (
00806                 HashSet<edge, Hash<edge> >::const_iterator iter =
00807                     regionConnectivity[procI].begin();
00808                 iter != regionConnectivity[procI].end();
00809                 ++iter
00810             )
00811             {
00812                 const edge& e = iter.key();
00813 
00814                 bool someLocal = false;
00815                 if (shiftIndexer.isLocal(e[0]))
00816                 {
00817                     label localI = shiftIndexer.toLocal(e[0]);
00818                     addUnique(e[1], regionRegions[localI]);
00819                     someLocal = true;
00820                 }
00821                 if (shiftIndexer.isLocal(e[1]))
00822                 {
00823                     label localI = shiftIndexer.toLocal(e[1]);
00824                     addUnique(e[0], regionRegions[localI]);
00825                     someLocal = true;
00826                 }
00827 
00828                 if (!someLocal)
00829                 {
00830                     FatalErrorIn("calcRegionRegions(..)")
00831                         << "Received from processor " << procI
00832                         << " connection " << e
00833                         << " where none of the elements is local to me."
00834                         << abort(FatalError);
00835                 }
00836             }
00837         }
00838     }
00839 }
00840 
00841 
00842 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00843 
00844 // Construct from components
00845 Foam::meshRefinement::meshRefinement
00846 (
00847     fvMesh& mesh,
00848     const scalar mergeDistance,
00849     const bool overwrite,
00850     const refinementSurfaces& surfaces,
00851     const shellSurfaces& shells
00852 )
00853 :
00854     mesh_(mesh),
00855     mergeDistance_(mergeDistance),
00856     overwrite_(overwrite),
00857     oldInstance_(mesh.pointsInstance()),
00858     surfaces_(surfaces),
00859     shells_(shells),
00860     meshCutter_
00861     (
00862         mesh,
00863         labelIOList
00864         (
00865             IOobject
00866             (
00867                 "cellLevel",
00868                 mesh_.facesInstance(),
00869                 fvMesh::meshSubDir,
00870                 mesh,
00871                 IOobject::READ_IF_PRESENT,
00872                 IOobject::NO_WRITE,
00873                 false
00874             ),
00875             labelList(mesh_.nCells(), 0)
00876         ),
00877         labelIOList
00878         (
00879             IOobject
00880             (
00881                 "pointLevel",
00882                 mesh_.facesInstance(),
00883                 fvMesh::meshSubDir,
00884                 mesh,
00885                 IOobject::READ_IF_PRESENT,
00886                 IOobject::NO_WRITE,
00887                 false
00888             ),
00889             labelList(mesh_.nPoints(), 0)
00890         ),
00891         refinementHistory
00892         (
00893             IOobject
00894             (
00895                 "refinementHistory",
00896                 mesh_.facesInstance(),
00897                 fvMesh::meshSubDir,
00898                 mesh_,
00899                 IOobject::NO_READ,
00900                 IOobject::NO_WRITE,
00901                 false
00902             ),
00903             List<refinementHistory::splitCell8>(0),
00904             labelList(0)
00905         )   // no unrefinement
00906     ),
00907     surfaceIndex_
00908     (
00909         IOobject
00910         (
00911             "surfaceIndex",
00912             mesh_.facesInstance(),
00913             fvMesh::meshSubDir,
00914             mesh_,
00915             IOobject::NO_READ,
00916             IOobject::NO_WRITE,
00917             false
00918         ),
00919         labelList(mesh_.nFaces(), -1)
00920     ),
00921     userFaceData_(0)
00922 {
00923     // recalculate intersections for all faces
00924     updateIntersections(identity(mesh_.nFaces()));
00925 }
00926 
00927 
00928 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00929 
00930 Foam::label Foam::meshRefinement::countHits() const
00931 {
00932     // Stats on edges to test. Count proc faces only once.
00933     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
00934 
00935     label nHits = 0;
00936 
00937     forAll(surfaceIndex_, faceI)
00938     {
00939         if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
00940         {
00941             nHits++;
00942         }
00943     }
00944     return nHits;
00945 }
00946 
00947 
00948 // Determine distribution to move connected regions onto one processor.
00949 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
00950 (
00951     const scalarField& cellWeights,
00952     const boolList& blockedFace,
00953     const List<labelPair>& explicitConnections,
00954     decompositionMethod& decomposer
00955 ) const
00956 {
00957     // Determine global regions, separated by blockedFaces
00958     regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
00959 
00960     // Now globalRegion has global region per cell. Problem is that
00961     // the region might span multiple domains so we want to get
00962     // a "region master" per domain. Note that multi-processor
00963     // regions can only occur on cells on coupled patches.
00964     // Note: since the number of regions does not change by this the
00965     // process can be seen as just a shift of a region onto a single
00966     // processor.
00967 
00968 
00969     // Global cell numbering engine
00970     globalIndex globalCells(mesh_.nCells());
00971 
00972 
00973     // Determine per coupled region the master cell (lowest numbered cell
00974     // on lowest numbered processor)
00975     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00976     // (does not determine master for non-coupled=fully-local regions)
00977 
00978     Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
00979     getCoupledRegionMaster
00980     (
00981         globalCells,
00982         blockedFace,
00983         globalRegion,
00984         coupledRegionToMaster
00985     );
00986 
00987     // Determine my regions
00988     // ~~~~~~~~~~~~~~~~~~~~
00989     // These are the fully local ones or the coupled ones of which I am master.
00990 
00991     Map<label> globalToLocalRegion;
00992     pointField localPoints;
00993     scalarField localWeights;
00994     calcLocalRegions
00995     (
00996         globalCells,
00997         globalRegion,
00998         coupledRegionToMaster,
00999         cellWeights,
01000 
01001         globalToLocalRegion,
01002         localPoints,
01003         localWeights
01004     );
01005 
01006 
01007 
01008     // Find distribution for regions
01009     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01010 
01011     labelList regionDistribution;
01012 
01013     if (isA<geomDecomp>(decomposer))
01014     {
01015         regionDistribution = decomposer.decompose(localPoints, localWeights);
01016     }
01017     else
01018     {
01019         labelListList regionRegions;
01020         calcRegionRegions
01021         (
01022             globalRegion,
01023             globalToLocalRegion,
01024             coupledRegionToMaster,
01025 
01026             regionRegions
01027         );
01028 
01029         regionDistribution = decomposer.decompose
01030         (
01031             regionRegions,
01032             localPoints,
01033             localWeights
01034         );
01035     }
01036 
01037 
01038 
01039     // Convert region-based decomposition back to cell-based one
01040     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01041 
01042     // Transfer destination processor back to all. Use global reduce for now.
01043     Map<label> regionToDist(coupledRegionToMaster.size());
01044     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
01045     {
01046         label region = iter.key();
01047 
01048         Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
01049 
01050         if (regionFnd != globalToLocalRegion.end())
01051         {
01052             // Master cell is local. Store my distribution.
01053             regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
01054         }
01055         else
01056         {
01057             // Master cell is not on this processor. Make sure it is overridden.
01058             regionToDist.insert(iter.key(), labelMax);
01059         }
01060     }
01061     Pstream::mapCombineGather(regionToDist, minEqOp<label>());
01062     Pstream::mapCombineScatter(regionToDist);
01063 
01064 
01065     // Determine destination for all cells
01066     labelList distribution(mesh_.nCells());
01067 
01068     forAll(globalRegion, cellI)
01069     {
01070         Map<label>::const_iterator fndRegion =
01071             regionToDist.find(globalRegion[cellI]);
01072 
01073         if (fndRegion != regionToDist.end())
01074         {
01075             distribution[cellI] = fndRegion();
01076         }
01077         else
01078         {
01079             // region is local to the processor.
01080             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
01081 
01082             distribution[cellI] = regionDistribution[localRegionI];
01083         }
01084     }
01085 
01086     return distribution;
01087 }
01088 
01089 
01090 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
01091 (
01092     const bool keepZoneFaces,
01093     const bool keepBaffles,
01094     const scalarField& cellWeights,
01095     decompositionMethod& decomposer,
01096     fvMeshDistribute& distributor
01097 )
01098 {
01099     autoPtr<mapDistributePolyMesh> map;
01100 
01101     if (Pstream::parRun())
01102     {
01103         //if (debug_)
01104         //{
01105         //    const_cast<Time&>(mesh_.time())++;
01106         //}
01107 
01108         // Wanted distribution
01109         labelList distribution;
01110 
01111         if (keepZoneFaces || keepBaffles)
01112         {
01113             // Faces where owner and neighbour are not 'connected' so can
01114             // go to different processors.
01115             boolList blockedFace(mesh_.nFaces(), true);
01116             label nUnblocked = 0;
01117             // Pairs of baffles
01118             List<labelPair> couples;
01119 
01120             if (keepZoneFaces)
01121             {
01122                 // Determine decomposition to keep/move surface zones
01123                 // on one processor. The reason is that snapping will make these
01124                 // into baffles, move and convert them back so if they were
01125                 // proc boundaries after baffling&moving the points might be no
01126                 // longer snychronised so recoupling will fail. To prevent this
01127                 // keep owner&neighbour of such a surface zone on the same
01128                 // processor.
01129 
01130                 const wordList& fzNames = surfaces().faceZoneNames();
01131                 const faceZoneMesh& fZones = mesh_.faceZones();
01132 
01133                 // Get faces whose owner and neighbour should stay together,
01134                 // i.e. they are not 'blocked'.
01135 
01136                 forAll(fzNames, surfI)
01137                 {
01138                     if (fzNames[surfI].size())
01139                     {
01140                         // Get zone
01141                         label zoneI = fZones.findZoneID(fzNames[surfI]);
01142 
01143                         const faceZone& fZone = fZones[zoneI];
01144 
01145                         forAll(fZone, i)
01146                         {
01147                             if (blockedFace[fZone[i]])
01148                             {
01149                                 blockedFace[fZone[i]] = false;
01150                                 nUnblocked++;
01151                             }
01152                         }
01153                     }
01154                 }
01155 
01156 
01157                 // If the faceZones are not synchronised the blockedFace
01158                 // might not be synchronised. If you are sure the faceZones
01159                 // are synchronised remove below check.
01160                 syncTools::syncFaceList
01161                 (
01162                     mesh_,
01163                     blockedFace,
01164                     andEqOp<bool>(),    // combine operator
01165                     false               // separation
01166                 );
01167             }
01168             reduce(nUnblocked, sumOp<label>());
01169 
01170             if (keepZoneFaces)
01171             {
01172                 Info<< "Found " << nUnblocked
01173                     << " zoned faces to keep together." << endl;
01174             }
01175 
01176             if (keepBaffles)
01177             {
01178                 // Get boundary baffles that need to stay together.
01179                 couples = getDuplicateFaces   // all baffles
01180                 (
01181                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
01182                    +mesh_.nInternalFaces()
01183                 );
01184             }
01185             label nCouples = returnReduce(couples.size(), sumOp<label>());
01186 
01187             if (keepBaffles)
01188             {
01189                 Info<< "Found " << nCouples << " baffles to keep together."
01190                     << endl;
01191             }
01192 
01193             if (nUnblocked > 0 || nCouples > 0)
01194             {
01195                 Info<< "Applying special decomposition to keep baffles"
01196                     << " and zoned faces together." << endl;
01197 
01198                 distribution = decomposeCombineRegions
01199                 (
01200                     cellWeights,
01201                     blockedFace,
01202                     couples,
01203                     decomposer
01204                 );
01205 
01206                 labelList nProcCells(distributor.countCells(distribution));
01207                 Pstream::listCombineGather(nProcCells, plusEqOp<label>());
01208                 Pstream::listCombineScatter(nProcCells);
01209 
01210                 Info<< "Calculated decomposition:" << endl;
01211                 forAll(nProcCells, procI)
01212                 {
01213                     Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
01214                 }
01215                 Info<< endl;
01216             }
01217             else
01218             {
01219                 // Normal decomposition
01220                 distribution = decomposer.decompose
01221                 (
01222                     mesh_.cellCentres(),
01223                     cellWeights
01224                 );
01225             }
01226         }
01227         else
01228         {
01229             // Normal decomposition
01230             distribution = decomposer.decompose
01231             (
01232                 mesh_.cellCentres(),
01233                 cellWeights
01234             );
01235         }
01236 
01237         if (debug)
01238         {
01239             labelList nProcCells(distributor.countCells(distribution));
01240             Pout<< "Wanted distribution:" << nProcCells << endl;
01241 
01242             Pstream::listCombineGather(nProcCells, plusEqOp<label>());
01243             Pstream::listCombineScatter(nProcCells);
01244 
01245             Pout<< "Wanted resulting decomposition:" << endl;
01246             forAll(nProcCells, procI)
01247             {
01248                 Pout<< "    " << procI << '\t' << nProcCells[procI] << endl;
01249             }
01250             Pout<< endl;
01251         }
01252         // Do actual sending/receiving of mesh
01253         map = distributor.distribute(distribution);
01254 
01255         // Update numbering of meshRefiner
01256         distribute(map);
01257     }
01258     return map;
01259 }
01260 
01261 
01262 // Helper function to get intersected faces
01263 Foam::labelList Foam::meshRefinement::intersectedFaces() const
01264 {
01265     label nBoundaryFaces = 0;
01266 
01267     forAll(surfaceIndex_, faceI)
01268     {
01269         if (surfaceIndex_[faceI] != -1)
01270         {
01271             nBoundaryFaces++;
01272         }
01273     }
01274 
01275     labelList surfaceFaces(nBoundaryFaces);
01276     nBoundaryFaces = 0;
01277 
01278     forAll(surfaceIndex_, faceI)
01279     {
01280         if (surfaceIndex_[faceI] != -1)
01281         {
01282             surfaceFaces[nBoundaryFaces++] = faceI;
01283         }
01284     }
01285     return surfaceFaces;
01286 }
01287 
01288 
01289 // Helper function to get points used by faces
01290 Foam::labelList Foam::meshRefinement::intersectedPoints() const
01291 {
01292     const faceList& faces = mesh_.faces();
01293 
01294     // Mark all points on faces that will become baffles
01295     PackedBoolList isBoundaryPoint(mesh_.nPoints());
01296     label nBoundaryPoints = 0;
01297 
01298     forAll(surfaceIndex_, faceI)
01299     {
01300         if (surfaceIndex_[faceI] != -1)
01301         {
01302             const face& f = faces[faceI];
01303 
01304             forAll(f, fp)
01305             {
01306                 if (isBoundaryPoint.set(f[fp], 1u))
01307                 {
01308                     nBoundaryPoints++;
01309                 }
01310             }
01311         }
01312     }
01313 
01315     //labelList adaptPatchIDs(meshedPatches());
01316     //forAll(adaptPatchIDs, i)
01317     //{
01318     //    label patchI = adaptPatchIDs[i];
01319     //
01320     //    if (patchI != -1)
01321     //    {
01322     //        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
01323     //
01324     //        label faceI = pp.start();
01325     //
01326     //        forAll(pp, i)
01327     //        {
01328     //            const face& f = faces[faceI];
01329     //
01330     //            forAll(f, fp)
01331     //            {
01332     //                if (isBoundaryPoint.set(f[fp], 1u))
01333     //                    nBoundaryPoints++;
01334     //                }
01335     //            }
01336     //            faceI++;
01337     //        }
01338     //    }
01339     //}
01340 
01341 
01342     // Pack
01343     labelList boundaryPoints(nBoundaryPoints);
01344     nBoundaryPoints = 0;
01345     forAll(isBoundaryPoint, pointI)
01346     {
01347         if (isBoundaryPoint.get(pointI) == 1u)
01348         {
01349             boundaryPoints[nBoundaryPoints++] = pointI;
01350         }
01351     }
01352 
01353     return boundaryPoints;
01354 }
01355 
01356 
01357 //- Create patch from set of patches
01358 Foam::autoPtr<Foam::indirectPrimitivePatch> Foam::meshRefinement::makePatch
01359 (
01360     const polyMesh& mesh,
01361     const labelList& patchIDs
01362 )
01363 {
01364     const polyBoundaryMesh& patches = mesh.boundaryMesh();
01365 
01366     // Count faces.
01367     label nFaces = 0;
01368 
01369     forAll(patchIDs, i)
01370     {
01371         const polyPatch& pp = patches[patchIDs[i]];
01372 
01373         nFaces += pp.size();
01374     }
01375 
01376     // Collect faces.
01377     labelList addressing(nFaces);
01378     nFaces = 0;
01379 
01380     forAll(patchIDs, i)
01381     {
01382         const polyPatch& pp = patches[patchIDs[i]];
01383 
01384         label meshFaceI = pp.start();
01385 
01386         forAll(pp, i)
01387         {
01388             addressing[nFaces++] = meshFaceI++;
01389         }
01390     }
01391 
01392     return autoPtr<indirectPrimitivePatch>
01393     (
01394         new indirectPrimitivePatch
01395         (
01396             IndirectList<face>(mesh.faces(), addressing),
01397             mesh.points()
01398         )
01399     );
01400 }
01401 
01402 
01403 // Construct pointVectorField with correct boundary conditions
01404 Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
01405 (
01406     const pointMesh& pMesh,
01407     const labelList& adaptPatchIDs
01408 )
01409 {
01410     const polyMesh& mesh = pMesh();
01411 
01412     // Construct displacement field.
01413     const pointBoundaryMesh& pointPatches = pMesh.boundary();
01414 
01415     wordList patchFieldTypes
01416     (
01417         pointPatches.size(),
01418         slipPointPatchVectorField::typeName
01419     );
01420 
01421     forAll(adaptPatchIDs, i)
01422     {
01423         patchFieldTypes[adaptPatchIDs[i]] =
01424             fixedValuePointPatchVectorField::typeName;
01425     }
01426 
01427     forAll(pointPatches, patchI)
01428     {
01429         if (isA<globalPointPatch>(pointPatches[patchI]))
01430         {
01431             patchFieldTypes[patchI] = globalPointPatchVectorField::typeName;
01432         }
01433         else if (isA<processorPointPatch>(pointPatches[patchI]))
01434         {
01435             patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
01436         }
01437     }
01438 
01439     tmp<pointVectorField> tfld
01440     (
01441         new pointVectorField
01442         (
01443             IOobject
01444             (
01445                 "pointDisplacement",
01446                 mesh.time().timeName(), //timeName(),
01447                 mesh,
01448                 IOobject::NO_READ,
01449                 IOobject::AUTO_WRITE
01450             ),
01451             pMesh,
01452             dimensionedVector("displacement", dimLength, vector::zero),
01453             patchFieldTypes
01454         )
01455     );
01456 
01457     return tfld;
01458 }
01459 
01460 
01461 void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
01462 {
01463     const faceZoneMesh& fZones = mesh.faceZones();
01464 
01465     // Check any zones are present anywhere and in same order
01466 
01467     {
01468         List<wordList> zoneNames(Pstream::nProcs());
01469         zoneNames[Pstream::myProcNo()] = fZones.names();
01470         Pstream::gatherList(zoneNames);
01471         Pstream::scatterList(zoneNames);
01472         // All have same data now. Check.
01473         forAll(zoneNames, procI)
01474         {
01475             if (procI != Pstream::myProcNo())
01476             {
01477                 if (zoneNames[procI] != zoneNames[Pstream::myProcNo()])
01478                 {
01479                     FatalErrorIn
01480                     (
01481                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
01482                     )   << "faceZones are not synchronised on processors." << nl
01483                         << "Processor " << procI << " has faceZones "
01484                         << zoneNames[procI] << nl
01485                         << "Processor " << Pstream::myProcNo()
01486                         << " has faceZones "
01487                         << zoneNames[Pstream::myProcNo()] << nl
01488                         << exit(FatalError);
01489                 }
01490             }
01491         }
01492     }
01493 
01494     // Check that coupled faces are present on both sides.
01495 
01496     labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
01497 
01498     forAll(fZones, zoneI)
01499     {
01500         const faceZone& fZone = fZones[zoneI];
01501 
01502         forAll(fZone, i)
01503         {
01504             label bFaceI = fZone[i]-mesh.nInternalFaces();
01505 
01506             if (bFaceI >= 0)
01507             {
01508                 if (faceToZone[bFaceI] == -1)
01509                 {
01510                     faceToZone[bFaceI] = zoneI;
01511                 }
01512                 else if (faceToZone[bFaceI] == zoneI)
01513                 {
01514                     FatalErrorIn
01515                     (
01516                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
01517                     )   << "Face " << fZone[i] << " in zone "
01518                         << fZone.name()
01519                         << " is twice in zone!"
01520                         << abort(FatalError);
01521                 }
01522                 else
01523                 {
01524                     FatalErrorIn
01525                     (
01526                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
01527                     )   << "Face " << fZone[i] << " in zone "
01528                         << fZone.name()
01529                         << " is also in zone "
01530                         << fZones[faceToZone[bFaceI]].name()
01531                         << abort(FatalError);
01532                 }
01533             }
01534         }
01535     }
01536 
01537     labelList neiFaceToZone(faceToZone);
01538     syncTools::swapBoundaryFaceList(mesh, neiFaceToZone, false);
01539 
01540     forAll(faceToZone, i)
01541     {
01542         if (faceToZone[i] != neiFaceToZone[i])
01543         {
01544             FatalErrorIn
01545             (
01546                 "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
01547             )   << "Face " << mesh.nInternalFaces()+i
01548                 << " is in zone " << faceToZone[i]
01549                 << ", its coupled face is in zone " << neiFaceToZone[i]
01550                 << abort(FatalError);
01551         }
01552     }
01553 }
01554 
01555 
01556 // Adds patch if not yet there. Returns patchID.
01557 Foam::label Foam::meshRefinement::addPatch
01558 (
01559     fvMesh& mesh,
01560     const word& patchName,
01561     const word& patchType
01562 )
01563 {
01564     polyBoundaryMesh& polyPatches =
01565         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
01566 
01567     label patchI = polyPatches.findPatchID(patchName);
01568     if (patchI != -1)
01569     {
01570         if (polyPatches[patchI].type() == patchType)
01571         {
01572             // Already there
01573             return patchI;
01574         }
01575         //else
01576         //{
01577         //    FatalErrorIn
01578         //    (
01579         //        "meshRefinement::addPatch(fvMesh&, const word&, const word&)"
01580         //    )   << "Patch " << patchName << " already exists but with type "
01581         //        << patchType << nl
01582         //        << "Current patch names:" << polyPatches.names()
01583         //        << exit(FatalError);
01584         //}
01585     }
01586 
01587 
01588     label insertPatchI = polyPatches.size();
01589     label startFaceI = mesh.nFaces();
01590 
01591     forAll(polyPatches, patchI)
01592     {
01593         const polyPatch& pp = polyPatches[patchI];
01594 
01595         if (isA<processorPolyPatch>(pp))
01596         {
01597             insertPatchI = patchI;
01598             startFaceI = pp.start();
01599             break;
01600         }
01601     }
01602 
01603 
01604     // Below is all quite a hack. Feel free to change once there is a better
01605     // mechanism to insert and reorder patches.
01606 
01607     // Clear local fields and e.g. polyMesh parallelInfo.
01608     mesh.clearOut();
01609 
01610     label sz = polyPatches.size();
01611 
01612     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
01613 
01614     // Add polyPatch at the end
01615     polyPatches.setSize(sz+1);
01616     polyPatches.set
01617     (
01618         sz,
01619         polyPatch::New
01620         (
01621             patchType,
01622             patchName,
01623             0,              // size
01624             startFaceI,
01625             insertPatchI,
01626             polyPatches
01627         )
01628     );
01629     fvPatches.setSize(sz+1);
01630     fvPatches.set
01631     (
01632         sz,
01633         fvPatch::New
01634         (
01635             polyPatches[sz],  // point to newly added polyPatch
01636             mesh.boundary()
01637         )
01638     );
01639 
01640     addPatchFields<volScalarField>
01641     (
01642         mesh,
01643         calculatedFvPatchField<scalar>::typeName
01644     );
01645     addPatchFields<volVectorField>
01646     (
01647         mesh,
01648         calculatedFvPatchField<vector>::typeName
01649     );
01650     addPatchFields<volSphericalTensorField>
01651     (
01652         mesh,
01653         calculatedFvPatchField<sphericalTensor>::typeName
01654     );
01655     addPatchFields<volSymmTensorField>
01656     (
01657         mesh,
01658         calculatedFvPatchField<symmTensor>::typeName
01659     );
01660     addPatchFields<volTensorField>
01661     (
01662         mesh,
01663         calculatedFvPatchField<tensor>::typeName
01664     );
01665 
01666     // Surface fields
01667 
01668     addPatchFields<surfaceScalarField>
01669     (
01670         mesh,
01671         calculatedFvPatchField<scalar>::typeName
01672     );
01673     addPatchFields<surfaceVectorField>
01674     (
01675         mesh,
01676         calculatedFvPatchField<vector>::typeName
01677     );
01678     addPatchFields<surfaceSphericalTensorField>
01679     (
01680         mesh,
01681         calculatedFvPatchField<sphericalTensor>::typeName
01682     );
01683     addPatchFields<surfaceSymmTensorField>
01684     (
01685         mesh,
01686         calculatedFvPatchField<symmTensor>::typeName
01687     );
01688     addPatchFields<surfaceTensorField>
01689     (
01690         mesh,
01691         calculatedFvPatchField<tensor>::typeName
01692     );
01693 
01694     // Create reordering list
01695     // patches before insert position stay as is
01696     labelList oldToNew(sz+1);
01697     for (label i = 0; i < insertPatchI; i++)
01698     {
01699         oldToNew[i] = i;
01700     }
01701     // patches after insert position move one up
01702     for (label i = insertPatchI; i < sz; i++)
01703     {
01704         oldToNew[i] = i+1;
01705     }
01706     // appended patch gets moved to insert position
01707     oldToNew[sz] = insertPatchI;
01708 
01709     // Shuffle into place
01710     polyPatches.reorder(oldToNew);
01711     fvPatches.reorder(oldToNew);
01712 
01713     reorderPatchFields<volScalarField>(mesh, oldToNew);
01714     reorderPatchFields<volVectorField>(mesh, oldToNew);
01715     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
01716     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
01717     reorderPatchFields<volTensorField>(mesh, oldToNew);
01718     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
01719     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
01720     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
01721     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
01722     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
01723 
01724     return insertPatchI;
01725 }
01726 
01727 
01728 Foam::label Foam::meshRefinement::addMeshedPatch
01729 (
01730     const word& name,
01731     const word& type
01732 )
01733 {
01734     label meshedI = findIndex(meshedPatches_, name);
01735 
01736     if (meshedI != -1)
01737     {
01738         // Already there. Get corresponding polypatch
01739         return mesh_.boundaryMesh().findPatchID(name);
01740     }
01741     else
01742     {
01743         // Add patch
01744         label patchI = addPatch(mesh_, name, type);
01745 
01746         // Store
01747         label sz = meshedPatches_.size();
01748         meshedPatches_.setSize(sz+1);
01749         meshedPatches_[sz] = name;
01750 
01751         return patchI;
01752     }
01753 }
01754 
01755 
01756 Foam::labelList Foam::meshRefinement::meshedPatches() const
01757 {
01758     labelList patchIDs(meshedPatches_.size());
01759     forAll(meshedPatches_, i)
01760     {
01761         patchIDs[i] = mesh_.boundaryMesh().findPatchID(meshedPatches_[i]);
01762 
01763         if (patchIDs[i] == -1)
01764         {
01765             FatalErrorIn("meshRefinement::meshedPatches() const")
01766                 << "Problem : did not find patch " << meshedPatches_[i]
01767                 << abort(FatalError);
01768         }
01769     }
01770 
01771     return patchIDs;
01772 }
01773 
01774 
01775 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
01776 (
01777     const point& keepPoint
01778 )
01779 {
01780     // Determine connected regions. regionSplit is the labelList with the
01781     // region per cell.
01782     regionSplit cellRegion(mesh_);
01783 
01784     label regionI = -1;
01785 
01786     label cellI = mesh_.findCell(keepPoint);
01787 
01788     if (cellI != -1)
01789     {
01790         regionI = cellRegion[cellI];
01791     }
01792 
01793     reduce(regionI, maxOp<label>());
01794 
01795     if (regionI == -1)
01796     {
01797         FatalErrorIn
01798         (
01799             "meshRefinement::splitMeshRegions(const point&)"
01800         )   << "Point " << keepPoint
01801             << " is not inside the mesh." << nl
01802             << "Bounding box of the mesh:" << mesh_.globalData().bb()
01803             << exit(FatalError);
01804     }
01805 
01806     // Subset
01807     // ~~~~~~
01808 
01809     // Get cells to remove
01810     DynamicList<label> cellsToRemove(mesh_.nCells());
01811     forAll(cellRegion, cellI)
01812     {
01813         if (cellRegion[cellI] != regionI)
01814         {
01815             cellsToRemove.append(cellI);
01816         }
01817     }
01818     cellsToRemove.shrink();
01819 
01820     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
01821     reduce(nCellsToKeep, sumOp<label>());
01822 
01823     Info<< "Keeping all cells in region " << regionI
01824         << " containing point " << keepPoint << endl
01825         << "Selected for keeping : "
01826         << nCellsToKeep
01827         << " cells." << endl;
01828 
01829 
01830     // Remove cells
01831     removeCells cellRemover(mesh_);
01832 
01833     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
01834 
01835     if (exposedFaces.size())
01836     {
01837         FatalErrorIn
01838         (
01839             "meshRefinement::splitMeshRegions(const point&)"
01840         )   << "Removing non-reachable cells should only expose boundary faces"
01841             << nl
01842             << "ExposedFaces:" << exposedFaces << abort(FatalError);
01843     }
01844 
01845     return doRemoveCells
01846     (
01847         cellsToRemove,
01848         exposedFaces,
01849         labelList(exposedFaces.size(),-1),  // irrelevant since 0 size.
01850         cellRemover
01851     );
01852 }
01853 
01854 
01855 void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
01856 {
01857     // mesh_ already distributed; distribute my member data
01858 
01859     // surfaceQueries_ ok.
01860 
01861     // refinement
01862     meshCutter_.distribute(map);
01863 
01864     // surfaceIndex is face data.
01865     map.distributeFaceData(surfaceIndex_);
01866 
01867     // maintainedFaces are indices of faces.
01868     forAll(userFaceData_, i)
01869     {
01870         map.distributeFaceData(userFaceData_[i].second());
01871     }
01872 
01873     // Redistribute surface and any fields on it.
01874     {
01875         Random rndGen(653213);
01876 
01877         // Get local mesh bounding box. Single box for now.
01878         List<treeBoundBox> meshBb(1);
01879         treeBoundBox& bb = meshBb[0];
01880         bb = treeBoundBox(mesh_.points());
01881         bb = bb.extend(rndGen, 1E-4);
01882 
01883         // Distribute all geometry (so refinementSurfaces and shellSurfaces)
01884         searchableSurfaces& geometry =
01885             const_cast<searchableSurfaces&>(surfaces_.geometry());
01886 
01887         forAll(geometry, i)
01888         {
01889             autoPtr<mapDistribute> faceMap;
01890             autoPtr<mapDistribute> pointMap;
01891             geometry[i].distribute
01892             (
01893                 meshBb,
01894                 false,          // do not keep outside triangles
01895                 faceMap,
01896                 pointMap
01897             );
01898 
01899             if (faceMap.valid())
01900             {
01901                 // (ab)use the instance() to signal current modification time
01902                 geometry[i].instance() = geometry[i].time().timeName();
01903             }
01904 
01905             faceMap.clear();
01906             pointMap.clear();
01907         }
01908     }
01909 }
01910 
01911 
01912 void Foam::meshRefinement::updateMesh
01913 (
01914     const mapPolyMesh& map,
01915     const labelList& changedFaces
01916 )
01917 {
01918     Map<label> dummyMap(0);
01919 
01920     updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
01921 }
01922 
01923 
01924 void Foam::meshRefinement::storeData
01925 (
01926     const labelList& pointsToStore,
01927     const labelList& facesToStore,
01928     const labelList& cellsToStore
01929 )
01930 {
01931     // For now only meshCutter has storable/retrievable data.
01932     meshCutter_.storeData
01933     (
01934         pointsToStore,
01935         facesToStore,
01936         cellsToStore
01937     );
01938 }
01939 
01940 
01941 void Foam::meshRefinement::updateMesh
01942 (
01943     const mapPolyMesh& map,
01944     const labelList& changedFaces,
01945     const Map<label>& pointsToRestore,
01946     const Map<label>& facesToRestore,
01947     const Map<label>& cellsToRestore
01948 )
01949 {
01950     // For now only meshCutter has storable/retrievable data.
01951 
01952     // Update numbering of cells/vertices.
01953     meshCutter_.updateMesh
01954     (
01955         map,
01956         pointsToRestore,
01957         facesToRestore,
01958         cellsToRestore
01959     );
01960 
01961     // Update surfaceIndex
01962     updateList(map.faceMap(), -1, surfaceIndex_);
01963 
01964     // Update cached intersection information
01965     updateIntersections(changedFaces);
01966 
01967     // Update maintained faces
01968     forAll(userFaceData_, i)
01969     {
01970         labelList& data = userFaceData_[i].second();
01971 
01972         if (userFaceData_[i].first() == KEEPALL)
01973         {
01974             // extend list with face-from-face data
01975             updateList(map.faceMap(), -1, data);
01976         }
01977         else if (userFaceData_[i].first() == MASTERONLY)
01978         {
01979             // keep master only
01980             labelList newFaceData(map.faceMap().size(), -1);
01981 
01982             forAll(newFaceData, faceI)
01983             {
01984                 label oldFaceI = map.faceMap()[faceI];
01985 
01986                 if (oldFaceI >= 0 && map.reverseFaceMap()[oldFaceI] == faceI)
01987                 {
01988                     newFaceData[faceI] = data[oldFaceI];
01989                 }
01990             }
01991             data.transfer(newFaceData);
01992         }
01993         else
01994         {
01995             // remove any face that has been refined i.e. referenced more than
01996             // once.
01997 
01998             // 1. Determine all old faces that get referenced more than once.
01999             // These get marked with -1 in reverseFaceMap
02000             labelList reverseFaceMap(map.reverseFaceMap());
02001 
02002             forAll(map.faceMap(), faceI)
02003             {
02004                 label oldFaceI = map.faceMap()[faceI];
02005 
02006                 if (oldFaceI >= 0)
02007                 {
02008                     if (reverseFaceMap[oldFaceI] != faceI)
02009                     {
02010                         // faceI is slave face. Mark old face.
02011                         reverseFaceMap[oldFaceI] = -1;
02012                     }
02013                 }
02014             }
02015 
02016             // 2. Map only faces with intact reverseFaceMap
02017             labelList newFaceData(map.faceMap().size(), -1);
02018             forAll(newFaceData, faceI)
02019             {
02020                 label oldFaceI = map.faceMap()[faceI];
02021 
02022                 if (oldFaceI >= 0)
02023                 {
02024                     if (reverseFaceMap[oldFaceI] == faceI)
02025                     {
02026                         newFaceData[faceI] = data[oldFaceI];
02027                     }
02028                 }
02029             }
02030             data.transfer(newFaceData);
02031         }
02032     }
02033 }
02034 
02035 
02036 bool Foam::meshRefinement::write() const
02037 {
02038     bool writeOk =
02039         mesh_.write()
02040      && meshCutter_.write()
02041      && surfaceIndex_.write();
02042 
02043 
02044     // Make sure that any distributed surfaces (so ones which probably have
02045     // been changed) get written as well.
02046     // Note: should ideally have some 'modified' flag to say whether it
02047     // has been changed or not.
02048     searchableSurfaces& geometry =
02049         const_cast<searchableSurfaces&>(surfaces_.geometry());
02050 
02051     forAll(geometry, i)
02052     {
02053         searchableSurface& s = geometry[i];
02054 
02055         // Check if instance() of surface is not constant or system.
02056         // Is good hint that surface is distributed.
02057         if
02058         (
02059             s.instance() != s.time().system()
02060          && s.instance() != s.time().caseSystem()
02061          && s.instance() != s.time().constant()
02062          && s.instance() != s.time().caseConstant()
02063         )
02064         {
02065             // Make sure it gets written to current time, not constant.
02066             s.instance() = s.time().timeName();
02067             writeOk = writeOk && s.write();
02068         }
02069     }
02070 
02071     return writeOk;
02072 }
02073 
02074 
02075 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
02076  const
02077 {
02078     const globalMeshData& pData = mesh_.globalData();
02079 
02080     if (debug)
02081     {
02082         Pout<< msg.c_str()
02083             << " : cells(local):" << mesh_.nCells()
02084             << "  faces(local):" << mesh_.nFaces()
02085             << "  points(local):" << mesh_.nPoints()
02086             << endl;
02087     }
02088     else
02089     {
02090         Info<< msg.c_str()
02091             << " : cells:" << pData.nTotalCells()
02092             << "  faces:" << pData.nTotalFaces()
02093             << "  points:" << pData.nTotalPoints()
02094             << endl;
02095     }
02096 
02097 
02098     //if (debug)
02099     {
02100         const labelList& cellLevel = meshCutter_.cellLevel();
02101 
02102         labelList nCells(gMax(cellLevel)+1, 0);
02103 
02104         forAll(cellLevel, cellI)
02105         {
02106             nCells[cellLevel[cellI]]++;
02107         }
02108 
02109         Pstream::listCombineGather(nCells, plusEqOp<label>());
02110         Pstream::listCombineScatter(nCells);
02111 
02112         Info<< "Cells per refinement level:" << endl;
02113         forAll(nCells, levelI)
02114         {
02115             Info<< "    " << levelI << '\t' << nCells[levelI]
02116                 << endl;
02117         }
02118     }
02119 }
02120 
02121 
02122 //- Return either time().constant() or oldInstance
02123 Foam::word Foam::meshRefinement::timeName() const
02124 {
02125     if (overwrite_ && mesh_.time().timeIndex() == 0)
02126     {
02127         return oldInstance_;
02128     }
02129     else
02130     {
02131         return mesh_.time().timeName();
02132     }
02133 }
02134 
02135 
02136 void Foam::meshRefinement::dumpRefinementLevel() const
02137 {
02138     volScalarField volRefLevel
02139     (
02140         IOobject
02141         (
02142             "cellLevel",
02143             timeName(),
02144             mesh_,
02145             IOobject::NO_READ,
02146             IOobject::AUTO_WRITE,
02147             false
02148         ),
02149         mesh_,
02150         dimensionedScalar("zero", dimless, 0),
02151         zeroGradientFvPatchScalarField::typeName
02152     );
02153 
02154     const labelList& cellLevel = meshCutter_.cellLevel();
02155 
02156     forAll(volRefLevel, cellI)
02157     {
02158         volRefLevel[cellI] = cellLevel[cellI];
02159     }
02160 
02161     volRefLevel.write();
02162 
02163 
02164     const pointMesh& pMesh = pointMesh::New(mesh_);
02165 
02166     pointScalarField pointRefLevel
02167     (
02168         IOobject
02169         (
02170             "pointLevel",
02171             timeName(),
02172             mesh_,
02173             IOobject::NO_READ,
02174             IOobject::NO_WRITE,
02175             false
02176         ),
02177         pMesh,
02178         dimensionedScalar("zero", dimless, 0)
02179     );
02180 
02181     const labelList& pointLevel = meshCutter_.pointLevel();
02182 
02183     forAll(pointRefLevel, pointI)
02184     {
02185         pointRefLevel[pointI] = pointLevel[pointI];
02186     }
02187 
02188     pointRefLevel.write();
02189 }
02190 
02191 
02192 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
02193 {
02194     {
02195         const pointField& cellCentres = mesh_.cellCentres();
02196 
02197         OFstream str(prefix + "_edges.obj");
02198         label vertI = 0;
02199         Pout<< "meshRefinement::dumpIntersections :"
02200             << " Writing cellcentre-cellcentre intersections to file "
02201             << str.name() << endl;
02202 
02203 
02204         // Redo all intersections
02205         // ~~~~~~~~~~~~~~~~~~~~~~
02206 
02207         // Get boundary face centre and level. Coupled aware.
02208         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
02209         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
02210         calcNeighbourData(neiLevel, neiCc);
02211 
02212         labelList intersectionFaces(intersectedFaces());
02213 
02214         // Collect segments we want to test for
02215         pointField start(intersectionFaces.size());
02216         pointField end(intersectionFaces.size());
02217 
02218         forAll(intersectionFaces, i)
02219         {
02220             label faceI = intersectionFaces[i];
02221             start[i] = cellCentres[mesh_.faceOwner()[faceI]];
02222 
02223             if (mesh_.isInternalFace(faceI))
02224             {
02225                 end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
02226             }
02227             else
02228             {
02229                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
02230             }
02231         }
02232 
02233         // Do tests in one go
02234         labelList surfaceHit;
02235         List<pointIndexHit> surfaceHitInfo;
02236         surfaces_.findAnyIntersection
02237         (
02238             start,
02239             end,
02240             surfaceHit,
02241             surfaceHitInfo
02242         );
02243 
02244         forAll(intersectionFaces, i)
02245         {
02246             if (surfaceHit[i] != -1)
02247             {
02248                 meshTools::writeOBJ(str, start[i]);
02249                 vertI++;
02250                 meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
02251                 vertI++;
02252                 meshTools::writeOBJ(str, end[i]);
02253                 vertI++;
02254                 str << "l " << vertI-2 << ' ' << vertI-1 << nl
02255                     << "l " << vertI-1 << ' ' << vertI << nl;
02256             }
02257         }
02258     }
02259 
02260     // Convert to vtk format
02261     string cmd
02262     (
02263         "objToVTK " + prefix + "_edges.obj " + prefix + "_edges.vtk > /dev/null"
02264     );
02265     system(cmd.c_str());
02266 
02267     Pout<< endl;
02268 }
02269 
02270 
02271 void Foam::meshRefinement::write
02272 (
02273     const label flag,
02274     const fileName& prefix
02275 ) const
02276 {
02277     if (flag & MESH)
02278     {
02279         write();
02280     }
02281     if (flag & SCALARLEVELS)
02282     {
02283         dumpRefinementLevel();
02284     }
02285     if (flag & OBJINTERSECTIONS && prefix.size())
02286     {
02287         dumpIntersections(prefix);
02288     }
02289 }
02290 
02291 
02292 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines