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

meshRefinementBaffles.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/fvMesh.H>
00028 #include <OpenFOAM/syncTools.H>
00029 #include <OpenFOAM/Time.H>
00030 #include <autoMesh/refinementSurfaces.H>
00031 #include <meshTools/pointSet.H>
00032 #include <meshTools/faceSet.H>
00033 #include <OpenFOAM/indirectPrimitivePatch.H>
00034 #include <dynamicMesh/polyTopoChange.H>
00035 #include <meshTools/meshTools.H>
00036 #include <dynamicMesh/polyModifyFace.H>
00037 #include <dynamicMesh/polyModifyCell.H>
00038 #include <dynamicMesh/polyAddFace.H>
00039 #include <dynamicMesh/polyRemoveFace.H>
00040 #include <dynamicMesh/polyAddPoint.H>
00041 #include <dynamicMesh/localPointRegion.H>
00042 #include <dynamicMesh/duplicatePoints.H>
00043 #include <OpenFOAM/OFstream.H>
00044 #include <meshTools/regionSplit.H>
00045 #include <dynamicMesh/removeCells.H>
00046 
00047 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00048 
00049 // Repatches external face or creates baffle for internal face
00050 // with user specified patches (might be different for both sides).
00051 // Returns label of added face.
00052 Foam::label Foam::meshRefinement::createBaffle
00053 (
00054     const label faceI,
00055     const label ownPatch,
00056     const label neiPatch,
00057     polyTopoChange& meshMod
00058 ) const
00059 {
00060     const face& f = mesh_.faces()[faceI];
00061     label zoneID = mesh_.faceZones().whichZone(faceI);
00062     bool zoneFlip = false;
00063 
00064     if (zoneID >= 0)
00065     {
00066         const faceZone& fZone = mesh_.faceZones()[zoneID];
00067         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
00068     }
00069 
00070     meshMod.setAction
00071     (
00072         polyModifyFace
00073         (
00074             f,                          // modified face
00075             faceI,                      // label of face
00076             mesh_.faceOwner()[faceI],   // owner
00077             -1,                         // neighbour
00078             false,                      // face flip
00079             ownPatch,                   // patch for face
00080             false,                      // remove from zone
00081             zoneID,                     // zone for face
00082             zoneFlip                    // face flip in zone
00083         )
00084     );
00085 
00086 
00087     label dupFaceI = -1;
00088 
00089     if (mesh_.isInternalFace(faceI))
00090     {
00091         if (neiPatch == -1)
00092         {
00093             FatalErrorIn
00094             (
00095                 "meshRefinement::createBaffle"
00096                 "(const label, const label, const label, polyTopoChange&)"
00097             )   << "No neighbour patch for internal face " << faceI
00098                 << " fc:" << mesh_.faceCentres()[faceI]
00099                 << " ownPatch:" << ownPatch << abort(FatalError);
00100         }
00101 
00102         bool reverseFlip = false;
00103         if (zoneID >= 0)
00104         {
00105             reverseFlip = !zoneFlip;
00106         }
00107 
00108         dupFaceI = meshMod.setAction
00109         (
00110             polyAddFace
00111             (
00112                 f.reverseFace(),            // modified face
00113                 mesh_.faceNeighbour()[faceI],// owner
00114                 -1,                         // neighbour
00115                 -1,                         // masterPointID
00116                 -1,                         // masterEdgeID
00117                 faceI,                      // masterFaceID,
00118                 true,                       // face flip
00119                 neiPatch,                   // patch for face
00120                 zoneID,                     // zone for face
00121                 reverseFlip                 // face flip in zone
00122             )
00123         );
00124     }
00125     return dupFaceI;
00126 }
00127 
00128 
00129 // Get an estimate for the patch the internal face should get. Bit heuristic.
00130 Foam::label Foam::meshRefinement::getBafflePatch
00131 (
00132     const labelList& facePatch,
00133     const label faceI
00134 ) const
00135 {
00136     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00137 
00138     // Loop over face points
00139     // for each point check all faces patch IDs
00140     // as soon as an ID >= 0 is found, break and assign that ID
00141     // to the current face.
00142     // Check first for real patch (so proper surface intersection and then
00143     // in facePatch array for patches to block off faces
00144 
00145     forAll(mesh_.faces()[faceI], fp)
00146     {
00147         label pointI = mesh_.faces()[faceI][fp];
00148 
00149         forAll(mesh_.pointFaces()[pointI], pf)
00150         {
00151             label pFaceI = mesh_.pointFaces()[pointI][pf];
00152 
00153             label patchI = patches.whichPatch(pFaceI);
00154 
00155             if (patchI != -1 && !patches[patchI].coupled())
00156             {
00157                 return patchI;
00158             }
00159             else if (facePatch[pFaceI] != -1)
00160             {
00161                 return facePatch[pFaceI];
00162             }
00163         }
00164     }
00165 
00166     // Loop over owner and neighbour cells, looking for the first face with a
00167     // valid patch number
00168     const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00169 
00170     forAll(ownFaces, i)
00171     {
00172         label cFaceI = ownFaces[i];
00173 
00174         label patchI = patches.whichPatch(cFaceI);
00175 
00176         if (patchI != -1 && !patches[patchI].coupled())
00177         {
00178             return patchI;
00179         }
00180         else if (facePatch[cFaceI] != -1)
00181         {
00182             return facePatch[cFaceI];
00183         }
00184     }
00185 
00186     if (mesh_.isInternalFace(faceI))
00187     {
00188         const cell& neiFaces = mesh_.cells()[mesh_.faceNeighbour()[faceI]];
00189 
00190         forAll(neiFaces, i)
00191         {
00192             label cFaceI = neiFaces[i];
00193 
00194             label patchI = patches.whichPatch(cFaceI);
00195 
00196             if (patchI != -1 && !patches[patchI].coupled())
00197             {
00198                 return patchI;
00199             }
00200             else if (facePatch[cFaceI] != -1)
00201             {
00202                 return facePatch[cFaceI];
00203             }
00204         }
00205     }
00206 
00207     WarningIn
00208     (
00209         "meshRefinement::getBafflePatch(const labelList&, const label)"
00210     )   << "Could not find boundary face neighbouring internal face "
00211         << faceI << " with face centre " << mesh_.faceCentres()[faceI]
00212         << nl
00213         << "Using arbitrary patch " << 0 << " instead." << endl;
00214 
00215     return 0;
00216 }
00217 
00218 
00219 // Determine patches for baffles.
00220 void Foam::meshRefinement::getBafflePatches
00221 (
00222     const labelList& globalToPatch,
00223     const labelList& neiLevel,
00224     const pointField& neiCc,
00225 
00226     labelList& ownPatch,
00227     labelList& neiPatch
00228 ) const
00229 {
00230     autoPtr<OFstream> str;
00231     label vertI = 0;
00232     if (debug&OBJINTERSECTIONS)
00233     {
00234         str.reset
00235         (
00236             new OFstream
00237             (
00238                 mesh_.time().path()/timeName()/"intersections.obj"
00239             )
00240         );
00241 
00242         Pout<< "getBafflePatches : Writing surface intersections to file "
00243             << str().name() << nl << endl;
00244     }
00245 
00246     const pointField& cellCentres = mesh_.cellCentres();
00247 
00248     // Build list of surfaces that are not to be baffled.
00249     const wordList& cellZoneNames = surfaces_.cellZoneNames();
00250 
00251     labelList surfacesToBaffle(cellZoneNames.size());
00252     label baffleI = 0;
00253     forAll(cellZoneNames, surfI)
00254     {
00255         if (cellZoneNames[surfI].size())
00256         {
00257             if (debug)
00258             {
00259                 Pout<< "getBafflePatches : Not baffling surface "
00260                     << surfaces_.names()[surfI] << endl;
00261             }
00262         }
00263         else
00264         {
00265             surfacesToBaffle[baffleI++] = surfI;
00266         }
00267     }
00268     surfacesToBaffle.setSize(baffleI);
00269 
00270     ownPatch.setSize(mesh_.nFaces());
00271     ownPatch = -1;
00272     neiPatch.setSize(mesh_.nFaces());
00273     neiPatch = -1;
00274 
00275 
00276     // Collect candidate faces
00277     // ~~~~~~~~~~~~~~~~~~~~~~~
00278 
00279     labelList testFaces(intersectedFaces());
00280 
00281     // Collect segments
00282     // ~~~~~~~~~~~~~~~~
00283 
00284     pointField start(testFaces.size());
00285     pointField end(testFaces.size());
00286 
00287     forAll(testFaces, i)
00288     {
00289         label faceI = testFaces[i];
00290 
00291         label own = mesh_.faceOwner()[faceI];
00292 
00293         if (mesh_.isInternalFace(faceI))
00294         {
00295             start[i] = cellCentres[own];
00296             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
00297         }
00298         else
00299         {
00300             start[i] = cellCentres[own];
00301             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
00302         }
00303     }
00304 
00305 
00306     // Do test for intersections
00307     // ~~~~~~~~~~~~~~~~~~~~~~~~~
00308     labelList surface1;
00309     List<pointIndexHit> hit1;
00310     labelList region1;
00311     labelList surface2;
00312     List<pointIndexHit> hit2;
00313     labelList region2;
00314     surfaces_.findNearestIntersection
00315     (
00316         surfacesToBaffle,
00317         start,
00318         end,
00319 
00320         surface1,
00321         hit1,
00322         region1,
00323 
00324         surface2,
00325         hit2,
00326         region2
00327     );
00328 
00329     forAll(testFaces, i)
00330     {
00331         label faceI = testFaces[i];
00332 
00333         if (hit1[i].hit() && hit2[i].hit())
00334         {
00335             if (str.valid())
00336             {
00337                 meshTools::writeOBJ(str(), start[i]);
00338                 vertI++;
00339                 meshTools::writeOBJ(str(), hit1[i].rawPoint());
00340                 vertI++;
00341                 meshTools::writeOBJ(str(), hit2[i].rawPoint());
00342                 vertI++;
00343                 meshTools::writeOBJ(str(), end[i]);
00344                 vertI++;
00345                 str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
00346                 str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
00347                 str()<< "l " << vertI-1 << ' ' << vertI << nl;
00348             }
00349 
00350             // Pick up the patches
00351             ownPatch[faceI] = globalToPatch
00352             [
00353                 surfaces_.globalRegion(surface1[i], region1[i])
00354             ];
00355             neiPatch[faceI] = globalToPatch
00356             [
00357                 surfaces_.globalRegion(surface2[i], region2[i])
00358             ];
00359 
00360             if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
00361             {
00362                 FatalErrorIn("getBafflePatches(..)")
00363                     << "problem." << abort(FatalError);
00364             }
00365         }
00366     }
00367 
00368     // No need to parallel sync since intersection data (surfaceIndex_ etc.)
00369     // already guaranteed to be synced...
00370     // However:
00371     // - owncc and neicc are reversed on different procs so might pick
00372     //   up different regions reversed? No problem. Neighbour on one processor
00373     //   might not be owner on the other processor but the neighbour is
00374     //   not used when creating baffles from proc faces.
00375     // - tolerances issues occasionally crop up.
00376     syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
00377     syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>(), false);
00378 }
00379 
00380 
00381 // Create baffle for every face where ownPatch != -1
00382 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
00383 (
00384     const labelList& ownPatch,
00385     const labelList& neiPatch
00386 )
00387 {
00388     if
00389     (
00390         ownPatch.size() != mesh_.nFaces()
00391      || neiPatch.size() != mesh_.nFaces()
00392     )
00393     {
00394         FatalErrorIn
00395         (
00396             "meshRefinement::createBaffles"
00397             "(const labelList&, const labelList&)"
00398         )   << "Illegal size :"
00399             << " ownPatch:" << ownPatch.size()
00400             << " neiPatch:" << neiPatch.size()
00401             << ". Should be number of faces:" << mesh_.nFaces()
00402             << abort(FatalError);
00403     }
00404 
00405     if (debug)
00406     {
00407         labelList syncedOwnPatch(ownPatch);
00408         syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>(), false);
00409         labelList syncedNeiPatch(neiPatch);
00410         syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>(), false);
00411 
00412         forAll(syncedOwnPatch, faceI)
00413         {
00414             if
00415             (
00416                 (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
00417              || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
00418             )
00419             {
00420                 FatalErrorIn
00421                 (
00422                     "meshRefinement::createBaffles"
00423                     "(const labelList&, const labelList&)"
00424                 )   << "Non synchronised at face:" << faceI
00425                     << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
00426                     << " fc:" << mesh_.faceCentres()[faceI] << endl
00427                     << "ownPatch:" << ownPatch[faceI]
00428                     << " syncedOwnPatch:" << syncedOwnPatch[faceI]
00429                     << " neiPatch:" << neiPatch[faceI]
00430                     << " syncedNeiPatch:" << syncedNeiPatch[faceI]
00431                     << abort(FatalError);
00432             }
00433         }
00434     }
00435 
00436     // Topochange container
00437     polyTopoChange meshMod(mesh_);
00438 
00439     label nBaffles = 0;
00440 
00441     forAll(ownPatch, faceI)
00442     {
00443         if (ownPatch[faceI] != -1)
00444         {
00445             // Create baffle or repatch face. Return label of inserted baffle
00446             // face.
00447             createBaffle
00448             (
00449                 faceI,
00450                 ownPatch[faceI],   // owner side patch
00451                 neiPatch[faceI],   // neighbour side patch
00452                 meshMod
00453             );
00454             nBaffles++;
00455         }
00456     }
00457 
00458     // Change the mesh (no inflation, parallel sync)
00459     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00460 
00461     // Update fields
00462     mesh_.updateMesh(map);
00463 
00464     // Move mesh if in inflation mode
00465     if (map().hasMotionPoints())
00466     {
00467         mesh_.movePoints(map().preMotionPoints());
00468     }
00469     else
00470     {
00471         // Delete mesh volumes.
00472         mesh_.clearOut();
00473     }
00474 
00475     if (overwrite())
00476     {
00477         mesh_.setInstance(oldInstance());
00478     }
00479 
00480     //- Redo the intersections on the newly create baffle faces. Note that
00481     //  this changes also the cell centre positions.
00482     faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
00483 
00484     const labelList& reverseFaceMap = map().reverseFaceMap();
00485     const labelList& faceMap = map().faceMap();
00486 
00487     // Pick up owner side of baffle
00488     forAll(ownPatch, oldFaceI)
00489     {
00490         label faceI = reverseFaceMap[oldFaceI];
00491 
00492         if (ownPatch[oldFaceI] != -1 && faceI >= 0)
00493         {
00494             const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00495 
00496             forAll(ownFaces, i)
00497             {
00498                 baffledFacesSet.insert(ownFaces[i]);
00499             }
00500         }
00501     }
00502     // Pick up neighbour side of baffle (added faces)
00503     forAll(faceMap, faceI)
00504     {
00505         label oldFaceI = faceMap[faceI];
00506 
00507         if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
00508         {
00509             const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00510 
00511             forAll(ownFaces, i)
00512             {
00513                 baffledFacesSet.insert(ownFaces[i]);
00514             }
00515         }
00516     }
00517     baffledFacesSet.sync(mesh_);
00518 
00519     updateMesh(map, baffledFacesSet.toc());
00520 
00521     return map;
00522 }
00523 
00524 
00525 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
00526 // (this information is recalculated instead of maintained since would be too
00527 // hard across splitMeshRegions).
00528 Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
00529 (
00530     const labelList& testFaces
00531 ) const
00532 {
00533     labelList duplicateFace
00534     (
00535         localPointRegion::findDuplicateFaces
00536         (
00537             mesh_,
00538             testFaces
00539         )
00540     );
00541 
00542     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00543 
00544     // Convert into list of coupled face pairs (mesh face labels).
00545     List<labelPair> duplicateFaces(testFaces.size());
00546     label dupI = 0;
00547 
00548     forAll(duplicateFace, i)
00549     {
00550         label otherFaceI = duplicateFace[i];
00551 
00552         if (otherFaceI != -1 && i < otherFaceI)
00553         {
00554             label meshFace0 = testFaces[i];
00555             label patch0 = patches.whichPatch(meshFace0);
00556             label meshFace1 = testFaces[otherFaceI];
00557             label patch1 = patches.whichPatch(meshFace1);
00558 
00559             if
00560             (
00561                 (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
00562              || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
00563             )
00564             {
00565                 FatalErrorIn
00566                 (
00567                     "meshRefinement::getDuplicateFaces"
00568                     "(const bool, const labelList&)"
00569                 )   << "One of two duplicate faces is on"
00570                     << " processorPolyPatch."
00571                     << "This is not allowed." << nl
00572                     << "Face:" << meshFace0
00573                     << " is on patch:" << patches[patch0].name()
00574                     << nl
00575                     << "Face:" << meshFace1
00576                     << " is on patch:" << patches[patch1].name()
00577                     << abort(FatalError);
00578             }
00579 
00580             duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
00581         }
00582     }
00583     duplicateFaces.setSize(dupI);
00584 
00585     Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
00586         << " pairs of duplicate faces." << nl << endl;
00587 
00588 
00589     if (debug)
00590     {
00591         faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
00592 
00593         forAll(duplicateFaces, i)
00594         {
00595             duplicateFaceSet.insert(duplicateFaces[i][0]);
00596             duplicateFaceSet.insert(duplicateFaces[i][1]);
00597         }
00598         Pout<< "Writing duplicate faces (baffles) to faceSet "
00599             << duplicateFaceSet.name() << nl << endl;
00600         duplicateFaceSet.write();
00601     }
00602 
00603     return duplicateFaces;
00604 }
00605 
00606 
00607 // Extract those baffles (duplicate) faces that are on the edge of a baffle
00608 // region. These are candidates for merging.
00609 // Done by counting the number of baffles faces per mesh edge. If edge
00610 // has 2 boundary faces and both are baffle faces it is the edge of a baffle
00611 // region.
00612 Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
00613 (
00614     const List<labelPair>& couples
00615 ) const
00616 {
00617     // Construct addressing engine for all duplicate faces (only one
00618     // for each pair)
00619 
00620     // Duplicate faces in mesh labels (first face of each pair only)
00621     // (reused later on to mark off filtered couples. see below)
00622     labelList duplicateFaces(couples.size());
00623 
00624 
00625     // All duplicate faces on edge of the patch are to be merged.
00626     // So we count for all edges of duplicate faces how many duplicate
00627     // faces use them.
00628     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
00629 
00630 
00631 
00632     // Count number of boundary faces per edge
00633     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00634 
00635     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00636 
00637     forAll(patches, patchI)
00638     {
00639         const polyPatch& pp = patches[patchI];
00640 
00641         // Count number of boundary faces. Discard coupled boundary faces.
00642         if (!pp.coupled())
00643         {
00644             label faceI = pp.start();
00645 
00646             forAll(pp, i)
00647             {
00648                 const labelList& fEdges = mesh_.faceEdges(faceI);
00649 
00650                 forAll(fEdges, fEdgeI)
00651                 {
00652                     nBafflesPerEdge[fEdges[fEdgeI]]++;
00653                 }
00654                 faceI++;
00655             }
00656         }
00657     }
00658 
00659 
00660     DynamicList<label> fe0;
00661     DynamicList<label> fe1;
00662 
00663 
00664     // Count number of duplicate boundary faces per edge
00665     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00666 
00667     forAll(couples, i)
00668     {
00669         const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
00670 
00671         forAll(fEdges0, fEdgeI)
00672         {
00673             nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
00674         }
00675 
00676         const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
00677 
00678         forAll(fEdges1, fEdgeI)
00679         {
00680             nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000;
00681         }
00682     }
00683 
00684     // Add nBaffles on shared edges
00685     syncTools::syncEdgeList
00686     (
00687         mesh_,
00688         nBafflesPerEdge,
00689         plusEqOp<label>(),  // in-place add
00690         0,                  // initial value
00691         false               // no separation
00692     );
00693 
00694 
00695     // Baffles which are not next to other boundaries and baffles will have
00696     // value 2*1000000+2*1
00697 
00698     List<labelPair> filteredCouples(couples.size());
00699     label filterI = 0;
00700 
00701     forAll(couples, i)
00702     {
00703         const labelPair& couple = couples[i];
00704 
00705         if
00706         (
00707             patches.whichPatch(couple.first())
00708          == patches.whichPatch(couple.second())
00709         )
00710         {
00711             const labelList& fEdges = mesh_.faceEdges(couples[i].first());
00712 
00713             forAll(fEdges, fEdgeI)
00714             {
00715                 label edgeI = fEdges[fEdgeI];
00716 
00717                 if (nBafflesPerEdge[edgeI] == 2*1000000+2*1)
00718                 {
00719                     filteredCouples[filterI++] = couples[i];
00720                     break;
00721                 }
00722             }
00723         }
00724     }
00725     filteredCouples.setSize(filterI);
00726 
00727     //Info<< "filterDuplicateFaces : from "
00728     //    << returnReduce(couples.size(), sumOp<label>())
00729     //    << " down to "
00730     //    << returnReduce(filteredCouples.size(), sumOp<label>())
00731     //    << " baffles." << nl << endl;
00732 
00733     return filteredCouples;
00734 }
00735 
00736 
00737 // Merge baffles. Gets pairs of faces.
00738 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
00739 (
00740     const List<labelPair>& couples
00741 )
00742 {
00743     // Mesh change engine
00744     polyTopoChange meshMod(mesh_);
00745 
00746     const faceList& faces = mesh_.faces();
00747     const labelList& faceOwner = mesh_.faceOwner();
00748     const faceZoneMesh& faceZones = mesh_.faceZones();
00749 
00750     forAll(couples, i)
00751     {
00752         label face0 = couples[i].first();
00753         label face1 = couples[i].second();
00754 
00755         // face1 < 0 signals a coupled face that has been converted to baffle.
00756 
00757         label own0 = faceOwner[face0];
00758         label own1 = faceOwner[face1];
00759 
00760         if (face1 < 0 || own0 < own1)
00761         {
00762             // Use face0 as the new internal face.
00763             label zoneID = faceZones.whichZone(face0);
00764             bool zoneFlip = false;
00765 
00766             if (zoneID >= 0)
00767             {
00768                 const faceZone& fZone = faceZones[zoneID];
00769                 zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
00770             }
00771 
00772             label nei = (face1 < 0 ? -1 : own1);
00773 
00774             meshMod.setAction(polyRemoveFace(face1));
00775             meshMod.setAction
00776             (
00777                 polyModifyFace
00778                 (
00779                     faces[face0],           // modified face
00780                     face0,                  // label of face being modified
00781                     own0,                   // owner
00782                     nei,                    // neighbour
00783                     false,                  // face flip
00784                     -1,                     // patch for face
00785                     false,                  // remove from zone
00786                     zoneID,                 // zone for face
00787                     zoneFlip                // face flip in zone
00788                 )
00789             );
00790         }
00791         else
00792         {
00793             // Use face1 as the new internal face.
00794             label zoneID = faceZones.whichZone(face1);
00795             bool zoneFlip = false;
00796 
00797             if (zoneID >= 0)
00798             {
00799                 const faceZone& fZone = faceZones[zoneID];
00800                 zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
00801             }
00802 
00803             meshMod.setAction(polyRemoveFace(face0));
00804             meshMod.setAction
00805             (
00806                 polyModifyFace
00807                 (
00808                     faces[face1],           // modified face
00809                     face1,                  // label of face being modified
00810                     own1,                   // owner
00811                     own0,                   // neighbour
00812                     false,                  // face flip
00813                     -1,                     // patch for face
00814                     false,                  // remove from zone
00815                     zoneID,                 // zone for face
00816                     zoneFlip                // face flip in zone
00817                 )
00818             );
00819         }
00820     }
00821 
00822     // Change the mesh (no inflation)
00823     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00824 
00825     // Update fields
00826     mesh_.updateMesh(map);
00827 
00828     // Move mesh (since morphing does not do this)
00829     if (map().hasMotionPoints())
00830     {
00831         mesh_.movePoints(map().preMotionPoints());
00832     }
00833     else
00834     {
00835         // Delete mesh volumes.
00836         mesh_.clearOut();
00837     }
00838 
00839     if (overwrite())
00840     {
00841         mesh_.setInstance(oldInstance());
00842     }
00843 
00844     // Update intersections. Recalculate intersections on merged faces since
00845     // this seems to give problems? Note: should not be nessecary since
00846     // baffles preserve intersections from when they were created.
00847     labelList newExposedFaces(2*couples.size());
00848     label newI = 0;
00849 
00850     forAll(couples, i)
00851     {
00852         label newFace0 = map().reverseFaceMap()[couples[i].first()];
00853         if (newFace0 != -1)
00854         {
00855             newExposedFaces[newI++] = newFace0;
00856         }
00857 
00858         label newFace1 = map().reverseFaceMap()[couples[i].second()];
00859         if (newFace1 != -1)
00860         {
00861             newExposedFaces[newI++] = newFace1;
00862         }
00863     }
00864     newExposedFaces.setSize(newI);
00865     updateMesh(map, newExposedFaces);
00866 
00867     return map;
00868 }
00869 
00870 
00871 // Finds region per cell for cells inside closed named surfaces
00872 void Foam::meshRefinement::findCellZoneGeometric
00873 (
00874     const labelList& closedNamedSurfaces,   // indices of closed surfaces
00875     labelList& namedSurfaceIndex,           // per face index of named surface
00876     const labelList& surfaceToCellZone,     // cell zone index per surface
00877 
00878     labelList& cellToZone
00879 ) const
00880 {
00881     const pointField& cellCentres = mesh_.cellCentres();
00882     const labelList& faceOwner = mesh_.faceOwner();
00883     const labelList& faceNeighbour = mesh_.faceNeighbour();
00884 
00885     // Check if cell centre is inside
00886     labelList insideSurfaces;
00887     surfaces_.findInside
00888     (
00889         closedNamedSurfaces,
00890         cellCentres,
00891         insideSurfaces
00892     );
00893 
00894     forAll(insideSurfaces, cellI)
00895     {
00896         if (cellToZone[cellI] == -2)
00897         {
00898             label surfI = insideSurfaces[cellI];
00899 
00900             if (surfI != -1)
00901             {
00902                 cellToZone[cellI] = surfaceToCellZone[surfI];
00903             }
00904         }
00905     }
00906 
00907 
00908     // Some cells with cell centres close to surface might have
00909     // had been put into wrong surface. Recheck with perturbed cell centre.
00910 
00911 
00912     // 1. Collect points
00913 
00914     // Count points to test.
00915     label nCandidates = 0;
00916     forAll(namedSurfaceIndex, faceI)
00917     {
00918         label surfI = namedSurfaceIndex[faceI];
00919 
00920         if (surfI != -1)
00921         {
00922             if (mesh_.isInternalFace(faceI))
00923             {
00924                 nCandidates += 2;
00925             }
00926             else
00927             {
00928                 nCandidates += 1;
00929             }
00930         }
00931     }
00932 
00933     // Collect points.
00934     pointField candidatePoints(nCandidates);
00935     nCandidates = 0;
00936     forAll(namedSurfaceIndex, faceI)
00937     {
00938         label surfI = namedSurfaceIndex[faceI];
00939 
00940         if (surfI != -1)
00941         {
00942             label own = faceOwner[faceI];
00943             const point& ownCc = cellCentres[own];
00944 
00945             if (mesh_.isInternalFace(faceI))
00946             {
00947                 label nei = faceNeighbour[faceI];
00948                 const point& neiCc = cellCentres[nei];
00949                 // Perturbed cc
00950                 const vector d = 1E-4*(neiCc - ownCc);
00951                 candidatePoints[nCandidates++] = ownCc-d;
00952                 candidatePoints[nCandidates++] = neiCc+d;
00953             }
00954             else
00955             {
00956                 const point& neiFc = mesh_.faceCentres()[faceI];
00957                 // Perturbed cc
00958                 const vector d = 1E-4*(neiFc - ownCc);
00959                 candidatePoints[nCandidates++] = ownCc-d;
00960             }
00961         }
00962     }
00963 
00964 
00965     // 2. Test points for inside
00966 
00967     surfaces_.findInside
00968     (
00969         closedNamedSurfaces,
00970         candidatePoints,
00971         insideSurfaces
00972     );
00973 
00974 
00975     // 3. Update zone information
00976 
00977     nCandidates = 0;
00978     forAll(namedSurfaceIndex, faceI)
00979     {
00980         label surfI = namedSurfaceIndex[faceI];
00981 
00982         if (surfI != -1)
00983         {
00984             label own = faceOwner[faceI];
00985 
00986             if (mesh_.isInternalFace(faceI))
00987             {
00988                 label ownSurfI = insideSurfaces[nCandidates++];
00989                 if (ownSurfI != -1)
00990                 {
00991                     cellToZone[own] = surfaceToCellZone[ownSurfI];
00992                 }
00993 
00994                 label neiSurfI = insideSurfaces[nCandidates++];
00995                 if (neiSurfI != -1)
00996                 {
00997                     label nei = faceNeighbour[faceI];
00998 
00999                     cellToZone[nei] = surfaceToCellZone[neiSurfI];
01000                 }
01001             }
01002             else
01003             {
01004                 label ownSurfI = insideSurfaces[nCandidates++];
01005                 if (ownSurfI != -1)
01006                 {
01007                     cellToZone[own] = surfaceToCellZone[ownSurfI];
01008                 }
01009             }
01010         }
01011     }
01012 
01013 
01014     // Adapt the namedSurfaceIndex
01015     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
01016     // for if any cells were not completely covered.
01017 
01018     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01019     {
01020         label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01021         label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
01022 
01023         if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
01024         {
01025             // Give face the zone of max cell zone
01026             namedSurfaceIndex[faceI] = findIndex
01027             (
01028                 surfaceToCellZone,
01029                 max(ownZone, neiZone)
01030             );
01031         }
01032     }
01033 
01034     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
01035     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01036 
01037     forAll(patches, patchI)
01038     {
01039         const polyPatch& pp = patches[patchI];
01040 
01041         if (pp.coupled())
01042         {
01043             forAll(pp, i)
01044             {
01045                 label faceI = pp.start()+i;
01046                 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01047                 neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
01048             }
01049         }
01050     }
01051     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
01052 
01053     forAll(patches, patchI)
01054     {
01055         const polyPatch& pp = patches[patchI];
01056 
01057         if (pp.coupled())
01058         {
01059             forAll(pp, i)
01060             {
01061                 label faceI = pp.start()+i;
01062                 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01063                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
01064 
01065                 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
01066                 {
01067                     // Give face the max cell zone
01068                     namedSurfaceIndex[faceI] = findIndex
01069                     (
01070                         surfaceToCellZone,
01071                         max(ownZone, neiZone)
01072                     );
01073                 }
01074             }
01075         }
01076     }
01077 
01078     // Sync
01079     syncTools::syncFaceList
01080     (
01081         mesh_,
01082         namedSurfaceIndex,
01083         maxEqOp<label>(),
01084         false
01085     );
01086 }
01087 
01088 
01089 bool Foam::meshRefinement::calcRegionToZone
01090 (
01091     const label surfZoneI,
01092     const label ownRegion,
01093     const label neiRegion,
01094 
01095     labelList& regionToCellZone
01096 ) const
01097 {
01098     bool changed = false;
01099 
01100     // Check whether inbetween different regions
01101     if (ownRegion != neiRegion)
01102     {
01103         // Jump. Change one of the sides to my type.
01104 
01105         // 1. Interface between my type and unset region.
01106         // Set region to keepRegion
01107 
01108         if (regionToCellZone[ownRegion] == -2)
01109         {
01110             if (regionToCellZone[neiRegion] == surfZoneI)
01111             {
01112                 // Face between unset and my region. Put unset
01113                 // region into keepRegion
01114                 regionToCellZone[ownRegion] = -1;
01115                 changed = true;
01116             }
01117             else if (regionToCellZone[neiRegion] != -2)
01118             {
01119                 // Face between unset and other region.
01120                 // Put unset region into my region
01121                 regionToCellZone[ownRegion] = surfZoneI;
01122                 changed = true;
01123             }
01124         }
01125         else if (regionToCellZone[neiRegion] == -2)
01126         {
01127             if (regionToCellZone[ownRegion] == surfZoneI)
01128             {
01129                 // Face between unset and my region. Put unset
01130                 // region into keepRegion
01131                 regionToCellZone[neiRegion] = -1;
01132                 changed = true;
01133             }
01134             else if (regionToCellZone[ownRegion] != -2)
01135             {
01136                 // Face between unset and other region.
01137                 // Put unset region into my region
01138                 regionToCellZone[neiRegion] = surfZoneI;
01139                 changed = true;
01140             }
01141         }
01142     }
01143     return changed;
01144 }
01145 
01146 
01147 // Finds region per cell. Assumes:
01148 // - region containing keepPoint does not go into a cellZone
01149 // - all other regions can be found by crossing faces marked in
01150 //   namedSurfaceIndex.
01151 void Foam::meshRefinement::findCellZoneTopo
01152 (
01153     const point& keepPoint,
01154     const labelList& namedSurfaceIndex,
01155     const labelList& surfaceToCellZone,
01156     labelList& cellToZone
01157 ) const
01158 {
01159     // Analyse regions. Reuse regionsplit
01160     boolList blockedFace(mesh_.nFaces());
01161 
01162     forAll(namedSurfaceIndex, faceI)
01163     {
01164         if (namedSurfaceIndex[faceI] == -1)
01165         {
01166             blockedFace[faceI] = false;
01167         }
01168         else
01169         {
01170             blockedFace[faceI] = true;
01171         }
01172     }
01173     // No need to sync since namedSurfaceIndex already is synced
01174 
01175     // Set region per cell based on walking
01176     regionSplit cellRegion(mesh_, blockedFace);
01177     blockedFace.clear();
01178 
01179     // Per mesh region the zone the cell should be put in.
01180     // -2   : not analysed yet
01181     // -1   : keepPoint region. Not put into any cellzone.
01182     // >= 0 : index of cellZone
01183     labelList regionToCellZone(cellRegion.nRegions(), -2);
01184 
01185     // See which cells already are set in the cellToZone (from geometric
01186     // searching) and use these to take over their zones.
01187     // Note: could be improved to count number of cells per region.
01188     forAll(cellToZone, cellI)
01189     {
01190         if (cellToZone[cellI] != -2)
01191         {
01192             regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
01193         }
01194     }
01195 
01196 
01197 
01198     // Find the region containing the keepPoint
01199     label keepRegionI = -1;
01200 
01201     label cellI = mesh_.findCell(keepPoint);
01202 
01203     if (cellI != -1)
01204     {
01205         keepRegionI = cellRegion[cellI];
01206     }
01207     reduce(keepRegionI, maxOp<label>());
01208 
01209     Info<< "Found point " << keepPoint << " in cell " << cellI
01210         << " in global region " << keepRegionI
01211         << " out of " << cellRegion.nRegions() << " regions." << endl;
01212 
01213     if (keepRegionI == -1)
01214     {
01215         FatalErrorIn
01216         (
01217             "meshRefinement::findCellZoneTopo"
01218             "(const point&, const labelList&, const labelList&, labelList&)"
01219         )   << "Point " << keepPoint
01220             << " is not inside the mesh." << nl
01221             << "Bounding box of the mesh:" << mesh_.globalData().bb()
01222             << exit(FatalError);
01223     }
01224 
01225     // Mark default region with zone -1.
01226     if (regionToCellZone[keepRegionI] == -2)
01227     {
01228         regionToCellZone[keepRegionI] = -1;
01229     }
01230 
01231 
01232     // Find correspondence between cell zone and surface
01233     // by changing cell zone every time we cross a surface.
01234     while (true)
01235     {
01236         // Synchronise regionToCellZone.
01237         // Note:
01238         // - region numbers are identical on all processors
01239         // - keepRegion is identical ,,
01240         // - cellZones are identical ,,
01241         // This done at top of loop to account for geometric matching
01242         // not being synchronised.
01243         Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
01244         Pstream::listCombineScatter(regionToCellZone);
01245 
01246 
01247         bool changed = false;
01248 
01249         // Internal faces
01250 
01251         for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01252         {
01253             label surfI = namedSurfaceIndex[faceI];
01254 
01255             if (surfI != -1)
01256             {
01257                 // Calculate region to zone from cellRegions on either side
01258                 // of internal face.
01259                 bool changedCell = calcRegionToZone
01260                 (
01261                     surfaceToCellZone[surfI],
01262                     cellRegion[mesh_.faceOwner()[faceI]],
01263                     cellRegion[mesh_.faceNeighbour()[faceI]],
01264                     regionToCellZone
01265                 );
01266 
01267                 changed = changed | changedCell;
01268             }
01269         }
01270 
01271         // Boundary faces
01272 
01273         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01274 
01275         // Get coupled neighbour cellRegion
01276         labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
01277         forAll(patches, patchI)
01278         {
01279             const polyPatch& pp = patches[patchI];
01280 
01281             if (pp.coupled())
01282             {
01283                 forAll(pp, i)
01284                 {
01285                     label faceI = pp.start()+i;
01286                     neiCellRegion[faceI-mesh_.nInternalFaces()] =
01287                         cellRegion[mesh_.faceOwner()[faceI]];
01288                 }
01289             }
01290         }
01291         syncTools::swapBoundaryFaceList(mesh_, neiCellRegion, false);
01292 
01293         // Calculate region to zone from cellRegions on either side of coupled
01294         // face.
01295         forAll(patches, patchI)
01296         {
01297             const polyPatch& pp = patches[patchI];
01298 
01299             if (pp.coupled())
01300             {
01301                 forAll(pp, i)
01302                 {
01303                     label faceI = pp.start()+i;
01304 
01305                     label surfI = namedSurfaceIndex[faceI];
01306 
01307                     if (surfI != -1)
01308                     {
01309                         bool changedCell = calcRegionToZone
01310                         (
01311                             surfaceToCellZone[surfI],
01312                             cellRegion[mesh_.faceOwner()[faceI]],
01313                             neiCellRegion[faceI-mesh_.nInternalFaces()],
01314                             regionToCellZone
01315                         );
01316 
01317                         changed = changed | changedCell;
01318                     }
01319                 }
01320             }
01321         }
01322 
01323         if (!returnReduce(changed, orOp<bool>()))
01324         {
01325             break;
01326         }
01327     }
01328 
01329 
01330     forAll(regionToCellZone, regionI)
01331     {
01332         label zoneI = regionToCellZone[regionI];
01333 
01334         if (zoneI ==  -2)
01335         {
01336             FatalErrorIn
01337             (
01338                 "meshRefinement::findCellZoneTopo"
01339                 "(const point&, const labelList&, const labelList&, labelList&)"
01340             )   << "For region " << regionI << " haven't set cell zone."
01341                 << exit(FatalError);
01342         }
01343     }
01344 
01345     if (debug)
01346     {
01347         forAll(regionToCellZone, regionI)
01348         {
01349             Pout<< "Region " << regionI
01350                 << " becomes cellZone:" << regionToCellZone[regionI]
01351                 << endl;
01352         }
01353     }
01354 
01355     // Rework into cellToZone
01356     forAll(cellToZone, cellI)
01357     {
01358         cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
01359     }
01360 }
01361 
01362 
01363 // Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
01364 // faces inbetween same cell zone.
01365 void Foam::meshRefinement::makeConsistentFaceIndex
01366 (
01367     const labelList& cellToZone,
01368     labelList& namedSurfaceIndex
01369 ) const
01370 {
01371     const labelList& faceOwner = mesh_.faceOwner();
01372     const labelList& faceNeighbour = mesh_.faceNeighbour();
01373 
01374     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01375     {
01376         label ownZone = cellToZone[faceOwner[faceI]];
01377         label neiZone = cellToZone[faceNeighbour[faceI]];
01378 
01379         if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
01380         {
01381             namedSurfaceIndex[faceI] = -1;
01382         }
01383         else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
01384         {
01385             FatalErrorIn("meshRefinement::zonify()")
01386                 << "Different cell zones on either side of face " << faceI
01387                 << " at " << mesh_.faceCentres()[faceI]
01388                 << " but face not marked with a surface."
01389                 << abort(FatalError);
01390         }
01391     }
01392 
01393     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01394 
01395     // Get coupled neighbour cellZone
01396     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
01397     forAll(patches, patchI)
01398     {
01399         const polyPatch& pp = patches[patchI];
01400 
01401         if (pp.coupled())
01402         {
01403             forAll(pp, i)
01404             {
01405                 label faceI = pp.start()+i;
01406                 neiCellZone[faceI-mesh_.nInternalFaces()] =
01407                     cellToZone[mesh_.faceOwner()[faceI]];
01408             }
01409         }
01410     }
01411     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
01412 
01413     // Use coupled cellZone to do check
01414     forAll(patches, patchI)
01415     {
01416         const polyPatch& pp = patches[patchI];
01417 
01418         if (pp.coupled())
01419         {
01420             forAll(pp, i)
01421             {
01422                 label faceI = pp.start()+i;
01423 
01424                 label ownZone = cellToZone[faceOwner[faceI]];
01425                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
01426 
01427                 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
01428                 {
01429                     namedSurfaceIndex[faceI] = -1;
01430                 }
01431                 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
01432                 {
01433                     FatalErrorIn("meshRefinement::zonify()")
01434                         << "Different cell zones on either side of face "
01435                         << faceI << " at " << mesh_.faceCentres()[faceI]
01436                         << " but face not marked with a surface."
01437                         << abort(FatalError);
01438                 }
01439             }
01440         }
01441     }
01442 }
01443 
01444 
01445 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
01446 
01447 // Split off unreachable areas of mesh.
01448 void Foam::meshRefinement::baffleAndSplitMesh
01449 (
01450     const bool handleSnapProblems,
01451     const bool removeEdgeConnectedCells,
01452     const scalarField& perpendicularAngle,
01453     const bool mergeFreeStanding,
01454     const dictionary& motionDict,
01455     Time& runTime,
01456     const labelList& globalToPatch,
01457     const point& keepPoint
01458 )
01459 {
01460     // Introduce baffles
01461     // ~~~~~~~~~~~~~~~~~
01462 
01463     // Split the mesh along internal faces wherever there is a pierce between
01464     // two cell centres.
01465 
01466     Info<< "Introducing baffles for "
01467         << returnReduce(countHits(), sumOp<label>())
01468         << " faces that are intersected by the surface." << nl << endl;
01469 
01470     // Swap neighbouring cell centres and cell level
01471     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
01472     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
01473     calcNeighbourData(neiLevel, neiCc);
01474 
01475     labelList ownPatch, neiPatch;
01476     getBafflePatches
01477     (
01478         globalToPatch,
01479         neiLevel,
01480         neiCc,
01481 
01482         ownPatch,
01483         neiPatch
01484     );
01485 
01486     if (debug)
01487     {
01488         runTime++;
01489     }
01490 
01491     createBaffles(ownPatch, neiPatch);
01492 
01493     if (debug)
01494     {
01495         // Debug:test all is still synced across proc patches
01496         checkData();
01497     }
01498 
01499     Info<< "Created baffles in = "
01500         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01501 
01502     printMeshInfo(debug, "After introducing baffles");
01503 
01504     if (debug)
01505     {
01506         Pout<< "Writing baffled mesh to time " << timeName()
01507             << endl;
01508         write(debug, runTime.path()/"baffles");
01509         Pout<< "Dumped debug data in = "
01510             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01511     }
01512 
01513 
01514     // Introduce baffles to delete problem cells
01515     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01516 
01517     // Create some additional baffles where we want surface cells removed.
01518 
01519     if (handleSnapProblems)
01520     {
01521         Info<< nl
01522             << "Introducing baffles to block off problem cells" << nl
01523             << "----------------------------------------------" << nl
01524             << endl;
01525 
01526         labelList facePatch
01527         (
01528             markFacesOnProblemCells
01529             (
01530                 motionDict,
01531                 removeEdgeConnectedCells,
01532                 perpendicularAngle,
01533                 globalToPatch
01534             )
01535             //markFacesOnProblemCellsGeometric(motionDict)
01536         );
01537         Info<< "Analyzed problem cells in = "
01538             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01539 
01540         if (debug)
01541         {
01542             faceSet problemTopo(mesh_, "problemFacesTopo", 100);
01543 
01544             const labelList facePatchTopo
01545             (
01546                 markFacesOnProblemCells
01547                 (
01548                     motionDict,
01549                     removeEdgeConnectedCells,
01550                     perpendicularAngle,
01551                     globalToPatch
01552                 )
01553             );
01554             forAll(facePatchTopo, faceI)
01555             {
01556                 if (facePatchTopo[faceI] != -1)
01557                 {
01558                     problemTopo.insert(faceI);
01559                 }
01560             }
01561             Pout<< "Dumping " << problemTopo.size()
01562                 << " problem faces to " << problemTopo.objectPath() << endl;
01563             problemTopo.write();
01564         }
01565 
01566         Info<< "Introducing baffles to delete problem cells." << nl << endl;
01567 
01568         if (debug)
01569         {
01570             runTime++;
01571         }
01572 
01573         // Create baffles with same owner and neighbour for now.
01574         createBaffles(facePatch, facePatch);
01575 
01576         if (debug)
01577         {
01578             // Debug:test all is still synced across proc patches
01579             checkData();
01580         }
01581         Info<< "Created baffles in = "
01582             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01583 
01584         printMeshInfo(debug, "After introducing baffles");
01585 
01586         if (debug)
01587         {
01588             Pout<< "Writing extra baffled mesh to time "
01589                 << timeName() << endl;
01590             write(debug, runTime.path()/"extraBaffles");
01591             Pout<< "Dumped debug data in = "
01592                 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01593         }
01594     }
01595 
01596 
01597     // Select part of mesh
01598     // ~~~~~~~~~~~~~~~~~~~
01599 
01600     Info<< nl
01601         << "Remove unreachable sections of mesh" << nl
01602         << "-----------------------------------" << nl
01603         << endl;
01604 
01605     if (debug)
01606     {
01607         runTime++;
01608     }
01609 
01610     splitMeshRegions(keepPoint);
01611 
01612     if (debug)
01613     {
01614         // Debug:test all is still synced across proc patches
01615         checkData();
01616     }
01617     Info<< "Split mesh in = "
01618         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01619 
01620     printMeshInfo(debug, "After subsetting");
01621 
01622     if (debug)
01623     {
01624         Pout<< "Writing subsetted mesh to time " << timeName()
01625             << endl;
01626         write(debug, runTime.path()/timeName());
01627         Pout<< "Dumped debug data in = "
01628             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01629     }
01630 
01631 
01632     // Merge baffles
01633     // ~~~~~~~~~~~~~
01634 
01635     if (mergeFreeStanding)
01636     {
01637         Info<< nl
01638             << "Merge free-standing baffles" << nl
01639             << "---------------------------" << nl
01640             << endl;
01641 
01642         if (debug)
01643         {
01644             runTime++;
01645         }
01646 
01647         // List of pairs of freestanding baffle faces.
01648         List<labelPair> couples
01649         (
01650             filterDuplicateFaces    // filter out freestanding baffles
01651             (
01652                 getDuplicateFaces   // get all baffles
01653                 (
01654                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
01655                    +mesh_.nInternalFaces()
01656                 )
01657             )
01658         );
01659 
01660         label nCouples = couples.size();
01661         reduce(nCouples, sumOp<label>());
01662 
01663         Info<< "Detected free-standing baffles : " << nCouples << endl;
01664 
01665         if (nCouples > 0)
01666         {
01667             // Actually merge baffles. Note: not exactly parallellized. Should
01668             // convert baffle faces into processor faces if they resulted
01669             // from them.
01670             mergeBaffles(couples);
01671 
01672             if (debug)
01673             {
01674                 // Debug:test all is still synced across proc patches
01675                 checkData();
01676             }
01677         }
01678         Info<< "Merged free-standing baffles in = "
01679             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01680     }
01681 }
01682 
01683 
01684 // Split off (with optional buffer layers) unreachable areas of mesh.
01685 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
01686 (
01687     const label nBufferLayers,
01688     const labelList& globalToPatch,
01689     const point& keepPoint
01690 )
01691 {
01692     // Determine patches to put intersections into
01693     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01694 
01695     // Swap neighbouring cell centres and cell level
01696     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
01697     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
01698     calcNeighbourData(neiLevel, neiCc);
01699 
01700     labelList ownPatch, neiPatch;
01701     getBafflePatches
01702     (
01703         globalToPatch,
01704         neiLevel,
01705         neiCc,
01706 
01707         ownPatch,
01708         neiPatch
01709     );
01710 
01711     // Analyse regions. Reuse regionsplit
01712     boolList blockedFace(mesh_.nFaces(), false);
01713 
01714     forAll(ownPatch, faceI)
01715     {
01716         if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
01717         {
01718             blockedFace[faceI] = true;
01719         }
01720     }
01721     syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>(), false);
01722 
01723     // Set region per cell based on walking
01724     regionSplit cellRegion(mesh_, blockedFace);
01725     blockedFace.clear();
01726 
01727     // Find the region containing the keepPoint
01728     label keepRegionI = -1;
01729 
01730     label cellI = mesh_.findCell(keepPoint);
01731 
01732     if (cellI != -1)
01733     {
01734         keepRegionI = cellRegion[cellI];
01735     }
01736     reduce(keepRegionI, maxOp<label>());
01737 
01738     Info<< "Found point " << keepPoint << " in cell " << cellI
01739         << " in global region " << keepRegionI
01740         << " out of " << cellRegion.nRegions() << " regions." << endl;
01741 
01742     if (keepRegionI == -1)
01743     {
01744         FatalErrorIn
01745         (
01746             "meshRefinement::splitMesh"
01747             "(const label, const labelList&, const point&)"
01748         )   << "Point " << keepPoint
01749             << " is not inside the mesh." << nl
01750             << "Bounding box of the mesh:" << mesh_.globalData().bb()
01751             << exit(FatalError);
01752     }
01753 
01754 
01755     // Walk out nBufferlayers from region split
01756     // (modifies cellRegion, ownPatch)
01757     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01758     // Takes over face patch onto points and then back to faces and cells
01759     // (so cell-face-point walk)
01760 
01761     const labelList& faceOwner = mesh_.faceOwner();
01762     const labelList& faceNeighbour = mesh_.faceNeighbour();
01763 
01764     // Patch for exposed faces for lack of anything sensible.
01765     label defaultPatch = 0;
01766     if (globalToPatch.size())
01767     {
01768         defaultPatch = globalToPatch[0];
01769     }
01770 
01771     for (label i = 0; i < nBufferLayers; i++)
01772     {
01773         // 1. From cells (via faces) to points
01774 
01775         labelList pointBaffle(mesh_.nPoints(), -1);
01776 
01777         forAll(faceNeighbour, faceI)
01778         {
01779             const face& f = mesh_.faces()[faceI];
01780 
01781             label ownRegion = cellRegion[faceOwner[faceI]];
01782             label neiRegion = cellRegion[faceNeighbour[faceI]];
01783 
01784             if (ownRegion == keepRegionI && neiRegion != keepRegionI)
01785             {
01786                 // Note max(..) since possibly regionSplit might have split
01787                 // off extra unreachable parts of mesh. Note: or can this only
01788                 // happen for boundary faces?
01789                 forAll(f, fp)
01790                 {
01791                     pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
01792                 }
01793             }
01794             else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
01795             {
01796                 label newPatchI = neiPatch[faceI];
01797                 if (newPatchI == -1)
01798                 {
01799                     newPatchI = max(defaultPatch, ownPatch[faceI]);
01800                 }
01801                 forAll(f, fp)
01802                 {
01803                     pointBaffle[f[fp]] = newPatchI;
01804                 }
01805             }
01806         }
01807         for
01808         (
01809             label faceI = mesh_.nInternalFaces();
01810             faceI < mesh_.nFaces();
01811             faceI++
01812         )
01813         {
01814             const face& f = mesh_.faces()[faceI];
01815 
01816             label ownRegion = cellRegion[faceOwner[faceI]];
01817 
01818             if (ownRegion == keepRegionI)
01819             {
01820                 forAll(f, fp)
01821                 {
01822                     pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
01823                 }
01824             }
01825         }
01826 
01827         // Sync
01828         syncTools::syncPointList
01829         (
01830             mesh_,
01831             pointBaffle,
01832             maxEqOp<label>(),
01833             -1,                 // null value
01834             false               // no separation
01835         );
01836 
01837 
01838         // 2. From points back to faces
01839 
01840         const labelListList& pointFaces = mesh_.pointFaces();
01841 
01842         forAll(pointFaces, pointI)
01843         {
01844             if (pointBaffle[pointI] != -1)
01845             {
01846                 const labelList& pFaces = pointFaces[pointI];
01847 
01848                 forAll(pFaces, pFaceI)
01849                 {
01850                     label faceI = pFaces[pFaceI];
01851 
01852                     if (ownPatch[faceI] == -1)
01853                     {
01854                         ownPatch[faceI] = pointBaffle[pointI];
01855                     }
01856                 }
01857             }
01858         }
01859         syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
01860 
01861 
01862         // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
01863 
01864         labelList newOwnPatch(ownPatch);
01865 
01866         forAll(ownPatch, faceI)
01867         {
01868             if (ownPatch[faceI] != -1)
01869             {
01870                 label own = faceOwner[faceI];
01871 
01872                 if (cellRegion[own] != keepRegionI)
01873                 {
01874                     cellRegion[own] = keepRegionI;
01875 
01876                     const cell& ownFaces = mesh_.cells()[own];
01877                     forAll(ownFaces, j)
01878                     {
01879                         if (ownPatch[ownFaces[j]] == -1)
01880                         {
01881                             newOwnPatch[ownFaces[j]] = ownPatch[faceI];
01882                         }
01883                     }
01884                 }
01885                 if (mesh_.isInternalFace(faceI))
01886                 {
01887                     label nei = faceNeighbour[faceI];
01888 
01889                     if (cellRegion[nei] != keepRegionI)
01890                     {
01891                         cellRegion[nei] = keepRegionI;
01892 
01893                         const cell& neiFaces = mesh_.cells()[nei];
01894                         forAll(neiFaces, j)
01895                         {
01896                             if (ownPatch[neiFaces[j]] == -1)
01897                             {
01898                                 newOwnPatch[neiFaces[j]] = ownPatch[faceI];
01899                             }
01900                         }
01901                     }
01902                 }
01903             }
01904         }
01905 
01906         ownPatch.transfer(newOwnPatch);
01907 
01908         syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
01909     }
01910 
01911 
01912 
01913     // Subset
01914     // ~~~~~~
01915 
01916     // Get cells to remove
01917     DynamicList<label> cellsToRemove(mesh_.nCells());
01918     forAll(cellRegion, cellI)
01919     {
01920         if (cellRegion[cellI] != keepRegionI)
01921         {
01922             cellsToRemove.append(cellI);
01923         }
01924     }
01925     cellsToRemove.shrink();
01926 
01927     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
01928     reduce(nCellsToKeep, sumOp<label>());
01929 
01930     Info<< "Keeping all cells in region " << keepRegionI
01931         << " containing point " << keepPoint << endl
01932         << "Selected for keeping : " << nCellsToKeep
01933         << " cells." << endl;
01934 
01935 
01936     // Remove cells
01937     removeCells cellRemover(mesh_);
01938 
01939     // Pick up patches for exposed faces
01940     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
01941     labelList exposedPatches(exposedFaces.size());
01942 
01943     forAll(exposedFaces, i)
01944     {
01945         label faceI = exposedFaces[i];
01946 
01947         if (ownPatch[faceI] != -1)
01948         {
01949             exposedPatches[i] = ownPatch[faceI];
01950         }
01951         else
01952         {
01953             WarningIn("meshRefinement::splitMesh(..)")
01954                 << "For exposed face " << faceI
01955                 << " fc:" << mesh_.faceCentres()[faceI]
01956                 << " found no patch." << endl
01957                 << "    Taking patch " << defaultPatch
01958                 << " instead." << endl;
01959             exposedPatches[i] = defaultPatch;
01960         }
01961     }
01962 
01963     return doRemoveCells
01964     (
01965         cellsToRemove,
01966         exposedFaces,
01967         exposedPatches,
01968         cellRemover
01969     );
01970 }
01971 
01972 
01973 // Find boundary points that connect to more than one cell region and
01974 // split them.
01975 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
01976 {
01977     // Topochange container
01978     polyTopoChange meshMod(mesh_);
01979 
01980 
01981     // Analyse which points need to be duplicated
01982     localPointRegion regionSide(mesh_);
01983 
01984     label nNonManifPoints = returnReduce
01985     (
01986         regionSide.meshPointMap().size(),
01987         sumOp<label>()
01988     );
01989 
01990     Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
01991         << " non-manifold points (out of "
01992         << mesh_.globalData().nTotalPoints()
01993         << ')' << endl;
01994 
01995     // Topo change engine
01996     duplicatePoints pointDuplicator(mesh_);
01997 
01998     // Insert changes into meshMod
01999     pointDuplicator.setRefinement(regionSide, meshMod);
02000 
02001     // Change the mesh (no inflation, parallel sync)
02002     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
02003 
02004     // Update fields
02005     mesh_.updateMesh(map);
02006 
02007     // Move mesh if in inflation mode
02008     if (map().hasMotionPoints())
02009     {
02010         mesh_.movePoints(map().preMotionPoints());
02011     }
02012     else
02013     {
02014         // Delete mesh volumes.
02015         mesh_.clearOut();
02016     }
02017 
02018     if (overwrite())
02019     {
02020         mesh_.setInstance(oldInstance());
02021     }
02022 
02023     // Update intersections. Is mapping only (no faces created, positions stay
02024     // same) so no need to recalculate intersections.
02025     updateMesh(map, labelList(0));
02026 
02027     return map;
02028 }
02029 
02030 
02031 // Zoning
02032 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
02033 (
02034     const point& keepPoint,
02035     const bool allowFreeStandingZoneFaces
02036 )
02037 {
02038     const wordList& cellZoneNames = surfaces_.cellZoneNames();
02039     const wordList& faceZoneNames = surfaces_.faceZoneNames();
02040 
02041     labelList namedSurfaces(surfaces_.getNamedSurfaces());
02042 
02043     boolList isNamedSurface(cellZoneNames.size(), false);
02044 
02045     forAll(namedSurfaces, i)
02046     {
02047         label surfI = namedSurfaces[i];
02048 
02049         isNamedSurface[surfI] = true;
02050 
02051         Info<< "Surface : " << surfaces_.names()[surfI] << nl
02052             << "    faceZone : " << faceZoneNames[surfI] << nl
02053             << "    cellZone : " << cellZoneNames[surfI] << endl;
02054     }
02055 
02056 
02057     // Add zones to mesh
02058 
02059     labelList surfaceToFaceZone(faceZoneNames.size(), -1);
02060     {
02061         faceZoneMesh& faceZones = mesh_.faceZones();
02062 
02063         forAll(namedSurfaces, i)
02064         {
02065             label surfI = namedSurfaces[i];
02066 
02067             label zoneI = faceZones.findZoneID(faceZoneNames[surfI]);
02068 
02069             if (zoneI == -1)
02070             {
02071                 zoneI = faceZones.size();
02072                 faceZones.setSize(zoneI+1);
02073                 faceZones.set
02074                 (
02075                     zoneI,
02076                     new faceZone
02077                     (
02078                         faceZoneNames[surfI],   //name
02079                         labelList(0),           //addressing
02080                         boolList(0),            //flipmap
02081                         zoneI,                  //index
02082                         faceZones               //faceZoneMesh
02083                     )
02084                 );
02085             }
02086 
02087             if (debug)
02088             {
02089                 Pout<< "Faces on " << surfaces_.names()[surfI]
02090                     << " will go into faceZone " << zoneI << endl;
02091             }
02092             surfaceToFaceZone[surfI] = zoneI;
02093         }
02094 
02095         // Check they are synced
02096         List<wordList> allFaceZones(Pstream::nProcs());
02097         allFaceZones[Pstream::myProcNo()] = faceZones.names();
02098         Pstream::gatherList(allFaceZones);
02099         Pstream::scatterList(allFaceZones);
02100 
02101         for (label procI = 1; procI < allFaceZones.size(); procI++)
02102         {
02103             if (allFaceZones[procI] != allFaceZones[0])
02104             {
02105                 FatalErrorIn
02106                 (
02107                     "meshRefinement::zonify"
02108                     "(const label, const point&)"
02109                 )   << "Zones not synchronised among processors." << nl
02110                     << " Processor0 has faceZones:" << allFaceZones[0]
02111                     << " , processor" << procI
02112                     << " has faceZones:" << allFaceZones[procI]
02113                     << exit(FatalError);
02114             }
02115         }
02116     }
02117 
02118     labelList surfaceToCellZone(cellZoneNames.size(), -1);
02119     {
02120         cellZoneMesh& cellZones = mesh_.cellZones();
02121 
02122         forAll(namedSurfaces, i)
02123         {
02124             label surfI = namedSurfaces[i];
02125 
02126             label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
02127 
02128             if (zoneI == -1)
02129             {
02130                 zoneI = cellZones.size();
02131                 cellZones.setSize(zoneI+1);
02132                 cellZones.set
02133                 (
02134                     zoneI,
02135                     new cellZone
02136                     (
02137                         cellZoneNames[surfI],   //name
02138                         labelList(0),           //addressing
02139                         zoneI,                  //index
02140                         cellZones               //cellZoneMesh
02141                     )
02142                 );
02143             }
02144 
02145             if (debug)
02146             {
02147                 Pout<< "Cells inside " << surfaces_.names()[surfI]
02148                     << " will go into cellZone " << zoneI << endl;
02149             }
02150             surfaceToCellZone[surfI] = zoneI;
02151         }
02152 
02153         // Check they are synced
02154         List<wordList> allCellZones(Pstream::nProcs());
02155         allCellZones[Pstream::myProcNo()] = cellZones.names();
02156         Pstream::gatherList(allCellZones);
02157         Pstream::scatterList(allCellZones);
02158 
02159         for (label procI = 1; procI < allCellZones.size(); procI++)
02160         {
02161             if (allCellZones[procI] != allCellZones[0])
02162             {
02163                 FatalErrorIn
02164                 (
02165                     "meshRefinement::zonify"
02166                     "(const label, const point&)"
02167                 )   << "Zones not synchronised among processors." << nl
02168                     << " Processor0 has cellZones:" << allCellZones[0]
02169                     << " , processor" << procI
02170                     << " has cellZones:" << allCellZones[procI]
02171                     << exit(FatalError);
02172             }
02173         }
02174     }
02175 
02176 
02177 
02178     const pointField& cellCentres = mesh_.cellCentres();
02179     const labelList& faceOwner = mesh_.faceOwner();
02180     const labelList& faceNeighbour = mesh_.faceNeighbour();
02181     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
02182 
02183 
02184     // Mark faces intersecting zoned surfaces
02185     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02186 
02187 
02188     // Like surfaceIndex_ but only for named surfaces.
02189     labelList namedSurfaceIndex(mesh_.nFaces(), -1);
02190 
02191     {
02192         // Statistics: number of faces per faceZone
02193         labelList nSurfFaces(faceZoneNames.size(), 0);
02194 
02195         // Swap neighbouring cell centres and cell level
02196         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
02197         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
02198         calcNeighbourData(neiLevel, neiCc);
02199 
02200         // Note: for all internal faces? internal + coupled?
02201         // Since zonify is run after baffling the surfaceIndex_ on baffles is
02202         // not synchronised across both baffle faces. Fortunately we don't
02203         // do zonify baffle faces anyway (they are normal boundary faces).
02204 
02205         // Collect candidate faces
02206         // ~~~~~~~~~~~~~~~~~~~~~~~
02207 
02208         labelList testFaces(intersectedFaces());
02209 
02210         // Collect segments
02211         // ~~~~~~~~~~~~~~~~
02212 
02213         pointField start(testFaces.size());
02214         pointField end(testFaces.size());
02215 
02216         forAll(testFaces, i)
02217         {
02218             label faceI = testFaces[i];
02219 
02220             if (mesh_.isInternalFace(faceI))
02221             {
02222                 start[i] = cellCentres[faceOwner[faceI]];
02223                 end[i] = cellCentres[faceNeighbour[faceI]];
02224             }
02225             else
02226             {
02227                 start[i] = cellCentres[faceOwner[faceI]];
02228                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
02229             }
02230         }
02231 
02232 
02233         // Do test for intersections
02234         // ~~~~~~~~~~~~~~~~~~~~~~~~~
02235         // Note that we intersect all intersected faces again. Could reuse
02236         // the information already in surfaceIndex_.
02237 
02238         labelList surface1;
02239         labelList surface2;
02240         {
02241             List<pointIndexHit> hit1;
02242             labelList region1;
02243             List<pointIndexHit> hit2;
02244             labelList region2;
02245             surfaces_.findNearestIntersection
02246             (
02247                 namedSurfaces,
02248                 start,
02249                 end,
02250 
02251                 surface1,
02252                 hit1,
02253                 region1,
02254                 surface2,
02255                 hit2,
02256                 region2
02257             );
02258         }
02259 
02260         forAll(testFaces, i)
02261         {
02262             label faceI = testFaces[i];
02263 
02264             if (surface1[i] != -1)
02265             {
02266                 // If both hit should probably choose nearest. For later.
02267                 namedSurfaceIndex[faceI] = surface1[i];
02268                 nSurfFaces[surface1[i]]++;
02269             }
02270             else if (surface2[i] != -1)
02271             {
02272                 namedSurfaceIndex[faceI] = surface2[i];
02273                 nSurfFaces[surface2[i]]++;
02274             }
02275         }
02276 
02277 
02278         // surfaceIndex migh have different surfaces on both sides if
02279         // there happen to be a (obviously thin) surface with different
02280         // regions between the cell centres. If one is on a named surface
02281         // and the other is not this might give problems so sync.
02282         syncTools::syncFaceList
02283         (
02284             mesh_,
02285             namedSurfaceIndex,
02286             maxEqOp<label>(),
02287             false
02288         );
02289 
02290         // Print a bit
02291         if (debug)
02292         {
02293             forAll(nSurfFaces, surfI)
02294             {
02295                 Pout<< "Surface:"
02296                     << surfaces_.names()[surfI]
02297                     << "  nZoneFaces:" << nSurfFaces[surfI] << nl;
02298             }
02299             Pout<< endl;
02300         }
02301     }
02302 
02303 
02304     // Put the cells into the correct zone
02305     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02306 
02307     // Closed surfaces with cellZone specified.
02308     labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
02309 
02310     // Zone per cell:
02311     // -2 : unset
02312     // -1 : not in any zone
02313     // >=0: zoneID
02314     labelList cellToZone(mesh_.nCells(), -2);
02315 
02316 
02317     // Set using geometric test
02318     // ~~~~~~~~~~~~~~~~~~~~~~~~
02319 
02320     if (closedNamedSurfaces.size())
02321     {
02322         findCellZoneGeometric
02323         (
02324             closedNamedSurfaces,   // indices of closed surfaces
02325             namedSurfaceIndex,     // per face index of named surface
02326             surfaceToCellZone,     // cell zone index per surface
02327             cellToZone
02328         );
02329     }
02330 
02331     // Set using walking
02332     // ~~~~~~~~~~~~~~~~~
02333 
02334     //if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
02335     {
02336         // Topological walk
02337         findCellZoneTopo
02338         (
02339             keepPoint,
02340             namedSurfaceIndex,
02341             surfaceToCellZone,
02342             cellToZone
02343         );
02344     }
02345 
02346 
02347     // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
02348     if (!allowFreeStandingZoneFaces)
02349     {
02350         makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
02351     }
02352 
02353     // Topochange container
02354     polyTopoChange meshMod(mesh_);
02355 
02356 
02357     // Put the faces into the correct zone
02358     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02359 
02360     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
02361     {
02362         label surfI = namedSurfaceIndex[faceI];
02363 
02364         if (surfI != -1)
02365         {
02366             // Orient face zone to have slave cells in max cell zone.
02367             label ownZone = cellToZone[faceOwner[faceI]];
02368             label neiZone = cellToZone[faceNeighbour[faceI]];
02369 
02370             bool flip;
02371             if (ownZone == max(ownZone, neiZone))
02372             {
02373                 flip = false;
02374             }
02375             else
02376             {
02377                 flip = true;
02378             }
02379 
02380             meshMod.setAction
02381             (
02382                 polyModifyFace
02383                 (
02384                     mesh_.faces()[faceI],           // modified face
02385                     faceI,                          // label of face
02386                     faceOwner[faceI],               // owner
02387                     faceNeighbour[faceI],           // neighbour
02388                     false,                          // face flip
02389                     -1,                             // patch for face
02390                     false,                          // remove from zone
02391                     surfaceToFaceZone[surfI],       // zone for face
02392                     flip                            // face flip in zone
02393                 )
02394             );
02395         }
02396     }
02397 
02398     // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
02399     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
02400     forAll(patches, patchI)
02401     {
02402         const polyPatch& pp = patches[patchI];
02403 
02404         if (pp.coupled())
02405         {
02406             forAll(pp, i)
02407             {
02408                 label faceI = pp.start()+i;
02409                 neiCellZone[faceI-mesh_.nInternalFaces()] =
02410                     cellToZone[mesh_.faceOwner()[faceI]];
02411             }
02412         }
02413     }
02414     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
02415 
02416     // Get per face whether is it master (of a coupled set of faces)
02417     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
02418 
02419     // Set owner as no-flip
02420     forAll(patches, patchI)
02421     {
02422         const polyPatch& pp = patches[patchI];
02423 
02424         label faceI = pp.start();
02425 
02426         forAll(pp, i)
02427         {
02428             label surfI = namedSurfaceIndex[faceI];
02429 
02430             if (surfI != -1)
02431             {
02432                 label ownZone = cellToZone[faceOwner[faceI]];
02433                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
02434 
02435                 bool flip;
02436 
02437                 label maxZone = max(ownZone, neiZone);
02438 
02439                 if (maxZone == -1)
02440                 {
02441                     flip = false;
02442                 }
02443                 else if (ownZone == neiZone)
02444                 {
02445                     // Free-standing zone face or coupled boundary. Keep master
02446                     // face unflipped.
02447                     flip = !isMasterFace[faceI];
02448                 }
02449                 else if (neiZone == maxZone)
02450                 {
02451                     flip = true;
02452                 }
02453                 else
02454                 {
02455                     flip = false;
02456                 }
02457 
02458                 meshMod.setAction
02459                 (
02460                     polyModifyFace
02461                     (
02462                         mesh_.faces()[faceI],           // modified face
02463                         faceI,                          // label of face
02464                         faceOwner[faceI],               // owner
02465                         -1,                             // neighbour
02466                         false,                          // face flip
02467                         patchI,                         // patch for face
02468                         false,                          // remove from zone
02469                         surfaceToFaceZone[surfI],       // zone for face
02470                         flip                            // face flip in zone
02471                     )
02472                 );
02473             }
02474             faceI++;
02475         }
02476     }
02477 
02478 
02479     // Put the cells into the correct zone
02480     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02481 
02482     forAll(cellToZone, cellI)
02483     {
02484         label zoneI = cellToZone[cellI];
02485 
02486         if (zoneI >= 0)
02487         {
02488             meshMod.setAction
02489             (
02490                 polyModifyCell
02491                 (
02492                     cellI,
02493                     false,          // removeFromZone
02494                     zoneI
02495                 )
02496             );
02497         }
02498     }
02499 
02500     // Change the mesh (no inflation, parallel sync)
02501     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
02502 
02503     // Update fields
02504     mesh_.updateMesh(map);
02505 
02506     // Move mesh if in inflation mode
02507     if (map().hasMotionPoints())
02508     {
02509         mesh_.movePoints(map().preMotionPoints());
02510     }
02511     else
02512     {
02513         // Delete mesh volumes.
02514         mesh_.clearOut();
02515     }
02516 
02517     if (overwrite())
02518     {
02519         mesh_.setInstance(oldInstance());
02520     }
02521 
02522     // None of the faces has changed, only the zones. Still...
02523     updateMesh(map, labelList());
02524 
02525     return map;
02526 }
02527 
02528 
02529 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines