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

combineFaces.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 "combineFaces.H"
00027 #include <OpenFOAM/polyMesh.H>
00028 #include "polyTopoChange.H"
00029 #include <dynamicMesh/polyRemoveFace.H>
00030 #include <dynamicMesh/polyAddFace.H>
00031 #include <dynamicMesh/polyModifyFace.H>
00032 #include <dynamicMesh/polyRemovePoint.H>
00033 #include <dynamicMesh/polyAddPoint.H>
00034 #include <OpenFOAM/syncTools.H>
00035 #include <meshTools/meshTools.H>
00036 
00037 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00038 
00039 namespace Foam
00040 {
00041 
00042 defineTypeNameAndDebug(combineFaces, 0);
00043 
00044 }
00045 
00046 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00047 
00048 // Test face for (almost) convexeness. Allows certain convexity before
00049 // complaining.
00050 // For every two consecutive edges calculate the normal. If it differs too
00051 // much from the average face normal complain.
00052 bool Foam::combineFaces::convexFace
00053 (
00054     const scalar minConcaveCos,
00055     const pointField& points,
00056     const face& f
00057 )
00058 {
00059     // Get outwards pointing normal of f.
00060     vector n = f.normal(points);
00061     n /= mag(n);
00062 
00063     // Get edge from f[0] to f[size-1];
00064     vector ePrev(points[f[0]] - points[f[f.size()-1]]);
00065     scalar magEPrev = mag(ePrev);
00066     ePrev /= magEPrev + VSMALL;
00067 
00068     forAll(f, fp0)
00069     {
00070         // Get vertex after fp
00071         label fp1 = f.fcIndex(fp0);
00072 
00073         // Normalized vector between two consecutive points
00074         vector e10(points[f[fp1]] - points[f[fp0]]);
00075         scalar magE10 = mag(e10);
00076         e10 /= magE10 + VSMALL;
00077 
00078         if (magEPrev > SMALL && magE10 > SMALL)
00079         {
00080             vector edgeNormal = ePrev ^ e10;
00081 
00082             if ((edgeNormal & n) < 0)
00083             {
00084                 // Concave. Check angle.
00085                 if ((ePrev & e10) < minConcaveCos)
00086                 {
00087                     return false;
00088                 }
00089             }
00090         }
00091 
00092         ePrev = e10;
00093         magEPrev = magE10;
00094     }
00095 
00096     // Not a single internal angle is concave so face is convex.
00097     return true;
00098 }
00099 
00100 
00101 // Determines if set of faces is valid to collapse into single face.
00102 bool Foam::combineFaces::validFace
00103 (
00104     const scalar minConcaveCos,
00105     const indirectPrimitivePatch& bigFace
00106 )
00107 {
00108     // Get outside vertices (in local vertex numbering)
00109     const labelListList& edgeLoops = bigFace.edgeLoops();
00110 
00111     if (edgeLoops.size() > 1)
00112     {
00113         return false;
00114     }
00115 
00116     // Check for convexness
00117     face f(getOutsideFace(bigFace));
00118 
00119     return convexFace(minConcaveCos, bigFace.points(), f);
00120 }
00121 
00122 
00123 void Foam::combineFaces::regioniseFaces
00124 (
00125     const scalar minCos,
00126     const label cellI,
00127     const labelList& cEdges,
00128     Map<label>& faceRegion
00129 ) const
00130 {
00131     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00132 
00133     forAll(cEdges, i)
00134     {
00135         label edgeI = cEdges[i];
00136 
00137         label f0, f1;
00138         meshTools::getEdgeFaces(mesh_, cellI, edgeI, f0, f1);
00139 
00140         label p0 = patches.whichPatch(f0);
00141         label p1 = patches.whichPatch(f1);
00142 
00143         // Face can be merged if
00144         // - same non-coupled patch
00145         // - small angle
00146         if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
00147         {
00148             vector f0Normal = mesh_.faceAreas()[f0];
00149             f0Normal /= mag(f0Normal);
00150             vector f1Normal = mesh_.faceAreas()[f1];
00151             f1Normal /= mag(f1Normal);
00152 
00153             if ((f0Normal&f1Normal) > minCos)
00154             {
00155                 Map<label>::const_iterator f0Fnd = faceRegion.find(f0);
00156 
00157                 label region0 = -1;
00158                 if (f0Fnd != faceRegion.end())
00159                 {
00160                     region0 = f0Fnd();
00161                 }
00162 
00163                 Map<label>::const_iterator f1Fnd = faceRegion.find(f1);
00164 
00165                 label region1 = -1;
00166                 if (f1Fnd != faceRegion.end())
00167                 {
00168                     region1 = f1Fnd();
00169                 }
00170 
00171                 if (region0 == -1)
00172                 {
00173                     if (region1 == -1)
00174                     {
00175                         label useRegion = faceRegion.size();
00176                         faceRegion.insert(f0, useRegion);
00177                         faceRegion.insert(f1, useRegion);
00178                     }
00179                     else
00180                     {
00181                         faceRegion.insert(f0, region1);
00182                     }
00183                 }
00184                 else
00185                 {
00186                     if (region1 == -1)
00187                     {
00188                         faceRegion.insert(f1, region0);
00189                     }
00190                     else if (region0 != region1)
00191                     {
00192                         // Merge the two regions
00193                         label useRegion = min(region0, region1);
00194                         label freeRegion = max(region0, region1);
00195 
00196                         forAllIter(Map<label>, faceRegion, iter)
00197                         {
00198                             if (iter() == freeRegion)
00199                             {
00200                                 iter() = useRegion;
00201                             }
00202                         }
00203                     }
00204                 }
00205             }
00206         }
00207     }
00208 }
00209 
00210 
00211 bool Foam::combineFaces::faceNeighboursValid
00212 (
00213     const label cellI,
00214     const Map<label>& faceRegion
00215 ) const
00216 {
00217     if (faceRegion.size() <= 1)
00218     {
00219         return true;
00220     }
00221 
00222     const cell& cFaces = mesh_.cells()[cellI];
00223 
00224     DynamicList<label> storage;
00225 
00226     // Test for face collapsing to edge since too many neighbours merged.
00227     forAll(cFaces, cFaceI)
00228     {
00229         label faceI = cFaces[cFaceI];
00230 
00231         if (!faceRegion.found(faceI))
00232         {
00233             const labelList& fEdges = mesh_.faceEdges(faceI, storage);
00234 
00235             // Count number of remaining faces neighbouring faceI. This has
00236             // to be 3 or more.
00237 
00238             // Unregioned neighbouring faces
00239             DynamicList<label> neighbourFaces(cFaces.size());
00240             // Regioned neighbouring faces
00241             labelHashSet neighbourRegions;
00242 
00243             forAll(fEdges, i)
00244             {
00245                 label edgeI = fEdges[i];
00246                 label nbrI = meshTools::otherFace(mesh_, cellI, faceI, edgeI);
00247 
00248                 Map<label>::const_iterator iter = faceRegion.find(nbrI);
00249 
00250                 if (iter == faceRegion.end())
00251                 {
00252                     if (findIndex(neighbourFaces, nbrI) == -1)
00253                     {
00254                         neighbourFaces.append(nbrI);
00255                     }
00256                 }
00257                 else
00258                 {
00259                     neighbourRegions.insert(iter());
00260                 }
00261             }
00262 
00263             if ((neighbourFaces.size()+neighbourRegions.size()) < 3)
00264             {
00265                 return false;
00266             }
00267         }
00268     }
00269     return true;
00270 }
00271 
00272 
00273 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00274 
00275 // Construct from mesh
00276 Foam::combineFaces::combineFaces
00277 (
00278     const polyMesh& mesh,
00279     const bool undoable
00280 )
00281 :
00282     mesh_(mesh),
00283     undoable_(undoable),
00284     masterFace_(0),
00285     faceSetsVertices_(0),
00286     savedPointLabels_(0),
00287     savedPoints_(0)
00288 {}
00289 
00290 
00291 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00292 
00293 Foam::labelListList Foam::combineFaces::getMergeSets
00294 (
00295     const scalar featureCos,
00296     const scalar minConcaveCos,
00297     const labelHashSet& boundaryCells
00298 ) const
00299 {
00300     // Lists of faces that can be merged.
00301     DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
00302     // Storage for on-the-fly cell-edge addressing.
00303     DynamicList<label> storage;
00304 
00305     // On all cells regionise the faces
00306     forAllConstIter(labelHashSet, boundaryCells, iter)
00307     {
00308         label cellI = iter.key();
00309 
00310         const cell& cFaces = mesh_.cells()[cellI];
00311 
00312         const labelList& cEdges = mesh_.cellEdges(cellI, storage);
00313 
00314         // Region per face
00315         Map<label> faceRegion(cFaces.size());
00316         regioniseFaces(featureCos, cellI, cEdges, faceRegion);
00317 
00318         // Now we have in faceRegion for every face the region with planar
00319         // face sharing the same region. We now check whether the resulting
00320         // sets cause a face
00321         // - to become a set of edges since too many faces are merged.
00322         // - to become convex
00323 
00324         if (faceNeighboursValid(cellI, faceRegion))
00325         {
00326             // Create region-to-faces addressing
00327             Map<labelList> regionToFaces(faceRegion.size());
00328 
00329             forAllConstIter(Map<label>, faceRegion, iter)
00330             {
00331                 label faceI = iter.key();
00332                 label region = iter();
00333 
00334                 Map<labelList>::iterator regionFnd = regionToFaces.find(region);
00335 
00336                 if (regionFnd != regionToFaces.end())
00337                 {
00338                     labelList& setFaces = regionFnd();
00339                     label sz = setFaces.size();
00340                     setFaces.setSize(sz+1);
00341                     setFaces[sz] = faceI;
00342                 }
00343                 else
00344                 {
00345                     regionToFaces.insert(region, labelList(1, faceI));
00346                 }
00347             }
00348 
00349             // For every set check if it forms a valid convex face
00350 
00351             forAllConstIter(Map<labelList>, regionToFaces, iter)
00352             {
00353                 // Make face out of setFaces
00354                 indirectPrimitivePatch bigFace
00355                 (
00356                     IndirectList<face>
00357                     (
00358                         mesh_.faces(),
00359                         iter()
00360                     ),
00361                     mesh_.points()
00362                 );
00363 
00364                 // Only store if -only one outside loop -which forms convex face
00365                 if (validFace(minConcaveCos, bigFace))
00366                 {
00367                     allFaceSets.append(iter());
00368                 }
00369             }
00370         }
00371     }
00372 
00373     return allFaceSets.shrink();
00374 }
00375 
00376 
00377 Foam::labelListList Foam::combineFaces::getMergeSets
00378 (
00379     const scalar featureCos,
00380     const scalar minConcaveCos
00381 ) const
00382 {
00383     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00384 
00385     // Pick up all cells on boundary
00386     labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
00387 
00388     forAll(patches, patchI)
00389     {
00390         const polyPatch& patch = patches[patchI];
00391 
00392         if (!patch.coupled())
00393         {
00394             forAll(patch, i)
00395             {
00396                 boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
00397             }
00398         }
00399     }
00400 
00401     return getMergeSets(featureCos, minConcaveCos, boundaryCells);
00402 }
00403 
00404 
00405 // Gets outside edgeloop as a face
00406 // - in same order as faces
00407 // - in mesh vertex labels
00408 Foam::face Foam::combineFaces::getOutsideFace
00409 (
00410     const indirectPrimitivePatch& fp
00411 )
00412 {
00413     if (fp.edgeLoops().size() != 1)
00414     {
00415         FatalErrorIn
00416         (
00417             "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
00418         )   << "Multiple outside loops:" << fp.edgeLoops()
00419             << abort(FatalError);
00420     }
00421 
00422     // Get first boundary edge. Since guaranteed one edgeLoop when in here this
00423     // edge must be on it.
00424     label bEdgeI = fp.nInternalEdges();
00425 
00426     const edge& e = fp.edges()[bEdgeI];
00427 
00428     const labelList& eFaces = fp.edgeFaces()[bEdgeI];
00429 
00430     if (eFaces.size() != 1)
00431     {
00432         FatalErrorIn
00433         (
00434             "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
00435         )   << "boundary edge:" << bEdgeI
00436             << " points:" << fp.meshPoints()[e[0]]
00437             << ' ' << fp.meshPoints()[e[1]]
00438             << " on indirectPrimitivePatch has " << eFaces.size()
00439             << " faces using it" << abort(FatalError);
00440     }
00441 
00442 
00443     // Outside loop
00444     const labelList& outsideLoop = fp.edgeLoops()[0];
00445 
00446 
00447     // Get order of edge e in outside loop
00448     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00449 
00450     // edgeLoopConsistent : true  = edge in same order as outsideloop
00451     //                      false = opposite order
00452     bool edgeLoopConsistent = false;
00453 
00454     {
00455         label index0 = findIndex(outsideLoop, e[0]);
00456         label index1 = findIndex(outsideLoop, e[1]);
00457 
00458         if (index0 == -1 || index1 == -1)
00459         {
00460             FatalErrorIn
00461             (
00462                 "combineFaces::getOutsideFace"
00463                 "(const indirectPrimitivePatch&)"
00464             )   << "Cannot find boundary edge:" << e
00465                 << " points:" << fp.meshPoints()[e[0]]
00466                 << ' ' << fp.meshPoints()[e[1]]
00467                 << " in edgeLoop:" << outsideLoop << abort(FatalError);
00468         }
00469         else if (index1 == outsideLoop.fcIndex(index0))
00470         {
00471             edgeLoopConsistent = true;
00472         }
00473         else if (index0 == outsideLoop.fcIndex(index1))
00474         {
00475             edgeLoopConsistent = false;
00476         }
00477         else
00478         {
00479             FatalErrorIn
00480             (
00481                 "combineFaces::getOutsideFace"
00482                 "(const indirectPrimitivePatch&)"
00483             )   << "Cannot find boundary edge:" << e
00484                 << " points:" << fp.meshPoints()[e[0]]
00485                 << ' ' << fp.meshPoints()[e[1]]
00486                 << " on consecutive points in edgeLoop:"
00487                 << outsideLoop << abort(FatalError);
00488         }
00489     }
00490 
00491 
00492     // Get face in local vertices
00493     const face& localF = fp.localFaces()[eFaces[0]];
00494 
00495     // Get order of edge in localF
00496     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
00497 
00498     // faceEdgeConsistent : true  = edge in same order as localF
00499     //                      false = opposite order
00500     bool faceEdgeConsistent = false;
00501 
00502     {
00503         // Find edge in face.
00504         label index = findIndex(fp.faceEdges()[eFaces[0]], bEdgeI);
00505 
00506         if (index == -1)
00507         {
00508             FatalErrorIn
00509             (
00510                 "combineFaces::getOutsideFace"
00511                 "(const indirectPrimitivePatch&)"
00512             )   << "Cannot find boundary edge:" << e
00513                 << " points:" << fp.meshPoints()[e[0]]
00514                 << ' ' << fp.meshPoints()[e[1]]
00515                 << " in face:" << eFaces[0]
00516                 << " edges:" << fp.faceEdges()[eFaces[0]]
00517                 << abort(FatalError);
00518         }
00519         else if (localF[index] == e[0] && localF.nextLabel(index) == e[1])
00520         {
00521             faceEdgeConsistent = true;
00522         }
00523         else if (localF[index] == e[1] && localF.nextLabel(index) == e[0])
00524         {
00525             faceEdgeConsistent = false;
00526         }
00527         else
00528         {
00529             FatalErrorIn
00530             (
00531                 "combineFaces::getOutsideFace"
00532                 "(const indirectPrimitivePatch&)"
00533             )   << "Cannot find boundary edge:" << e
00534                 << " points:" << fp.meshPoints()[e[0]]
00535                 << ' ' << fp.meshPoints()[e[1]]
00536                 << " in face:" << eFaces[0] << " verts:" << localF
00537                 << abort(FatalError);
00538         }
00539     }
00540 
00541     // Get face in mesh points.
00542     face meshFace(renumber(fp.meshPoints(), outsideLoop));
00543 
00544     if (faceEdgeConsistent != edgeLoopConsistent)
00545     {
00546         reverse(meshFace);
00547     }
00548     return meshFace;
00549 }
00550 
00551 
00552 void Foam::combineFaces::setRefinement
00553 (
00554     const labelListList& faceSets,
00555     polyTopoChange& meshMod
00556 )
00557 {
00558     if (undoable_)
00559     {
00560         masterFace_.setSize(faceSets.size());
00561         faceSetsVertices_.setSize(faceSets.size());
00562         forAll(faceSets, setI)
00563         {
00564             const labelList& setFaces = faceSets[setI];
00565 
00566             masterFace_[setI] = setFaces[0];
00567             faceSetsVertices_[setI] = UIndirectList<face>
00568             (
00569                 mesh_.faces(),
00570                 setFaces
00571             );
00572         }
00573     }
00574 
00575     // Running count of number of faces using a point
00576     labelList nPointFaces(mesh_.nPoints(), 0);
00577 
00578     const labelListList& pointFaces = mesh_.pointFaces();
00579 
00580     forAll(pointFaces, pointI)
00581     {
00582         nPointFaces[pointI] = pointFaces[pointI].size();
00583     }
00584 
00585     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00586 
00587     forAll(faceSets, setI)
00588     {
00589         const labelList& setFaces = faceSets[setI];
00590 
00591         // Check
00592         if (debug)
00593         {
00594             forAll(setFaces, i)
00595             {
00596                 label patchI = patches.whichPatch(setFaces[i]);
00597 
00598                 if (patchI == -1 || patches[patchI].coupled())
00599                 {
00600                     FatalErrorIn
00601                     (
00602                         "combineFaces::setRefinement"
00603                         "(const bool, const labelListList&"
00604                         ", polyTopoChange&)"
00605                     )   << "Can only merge non-coupled boundary faces"
00606                         << " but found internal or coupled face:"
00607                         << setFaces[i] << " in set " << setI
00608                         << abort(FatalError);
00609                 }
00610             }
00611         }
00612 
00613         // Make face out of setFaces
00614         indirectPrimitivePatch bigFace
00615         (
00616             IndirectList<face>
00617             (
00618                 mesh_.faces(),
00619                 setFaces
00620             ),
00621             mesh_.points()
00622         );
00623 
00624         // Get outside vertices (in local vertex numbering)
00625         const labelListList& edgeLoops = bigFace.edgeLoops();
00626 
00627         if (edgeLoops.size() != 1)
00628         {
00629             FatalErrorIn
00630             (
00631                 "combineFaces::setRefinement"
00632                 "(const bool, const labelListList&, polyTopoChange&)"
00633             )   << "Faces to-be-merged " << setFaces
00634                 << " do not form a single big face." << nl
00635                 << abort(FatalError);
00636         }
00637 
00638 
00639         // Delete all faces except master, whose face gets modified.
00640 
00641         // Modify master face
00642         // ~~~~~~~~~~~~~~~~~~
00643 
00644         label masterFaceI = setFaces[0];
00645 
00646         // Get outside face in mesh vertex labels
00647         face outsideFace(getOutsideFace(bigFace));
00648 
00649         label zoneID = mesh_.faceZones().whichZone(masterFaceI);
00650 
00651         bool zoneFlip = false;
00652 
00653         if (zoneID >= 0)
00654         {
00655             const faceZone& fZone = mesh_.faceZones()[zoneID];
00656 
00657             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
00658         }
00659 
00660         label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
00661 
00662         meshMod.setAction
00663         (
00664             polyModifyFace
00665             (
00666                 outsideFace,                    // modified face
00667                 masterFaceI,                    // label of face being modified
00668                 mesh_.faceOwner()[masterFaceI], // owner
00669                 -1,                             // neighbour
00670                 false,                          // face flip
00671                 patchI,                         // patch for face
00672                 false,                          // remove from zone
00673                 zoneID,                         // zone for face
00674                 zoneFlip                        // face flip in zone
00675             )
00676         );
00677 
00678 
00679         // Delete all non-master faces
00680         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
00681 
00682         for (label i = 1; i < setFaces.size(); i++)
00683         {
00684             meshMod.setAction(polyRemoveFace(setFaces[i]));
00685         }
00686 
00687 
00688         // Mark unused points for deletion
00689         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00690 
00691         // Remove count of all faces combined
00692         forAll(setFaces, i)
00693         {
00694             const face& f = mesh_.faces()[setFaces[i]];
00695 
00696             forAll(f, fp)
00697             {
00698                 nPointFaces[f[fp]]--;
00699             }
00700         }
00701         // But keep count on new outside face
00702         forAll(outsideFace, fp)
00703         {
00704             nPointFaces[outsideFace[fp]]++;
00705         }
00706     }
00707 
00708 
00709     // If one or more processors use the point it needs to be kept.
00710     syncTools::syncPointList
00711     (
00712         mesh_,
00713         nPointFaces,
00714         plusEqOp<label>(),
00715         0,                  // null value
00716         false               // no separation
00717     );
00718 
00719     // Remove all unused points. Store position if undoable.
00720     if (!undoable_)
00721     {
00722         forAll(nPointFaces, pointI)
00723         {
00724             if (nPointFaces[pointI] == 0)
00725             {
00726                 meshMod.setAction(polyRemovePoint(pointI));
00727             }
00728         }
00729     }
00730     else
00731     {
00732         // Count removed points
00733         label n = 0;
00734         forAll(nPointFaces, pointI)
00735         {
00736             if (nPointFaces[pointI] == 0)
00737             {
00738                 n++;
00739             }
00740         }
00741 
00742         savedPointLabels_.setSize(n);
00743         savedPoints_.setSize(n);
00744         Map<label> meshToSaved(2*n);
00745 
00746         // Remove points and store position
00747         n = 0;
00748         forAll(nPointFaces, pointI)
00749         {
00750             if (nPointFaces[pointI] == 0)
00751             {
00752                 meshMod.setAction(polyRemovePoint(pointI));
00753 
00754                 savedPointLabels_[n] = pointI;
00755                 savedPoints_[n] = mesh_.points()[pointI];
00756 
00757                 meshToSaved.insert(pointI, n);
00758                 n++;
00759             }
00760         }
00761 
00762         // Update stored vertex labels. Negative indices index into local points
00763         forAll(faceSetsVertices_, setI)
00764         {
00765             faceList& setFaces = faceSetsVertices_[setI];
00766 
00767             forAll(setFaces, i)
00768             {
00769                 face& f = setFaces[i];
00770 
00771                 forAll(f, fp)
00772                 {
00773                     label pointI = f[fp];
00774 
00775                     if (nPointFaces[pointI] == 0)
00776                     {
00777                         f[fp] = -meshToSaved[pointI]-1;
00778                     }
00779                 }
00780             }
00781         }
00782     }
00783 }
00784 
00785 
00786 void Foam::combineFaces::updateMesh(const mapPolyMesh& map)
00787 {
00788     if (undoable_)
00789     {
00790         // Master face just renumbering of point labels
00791         inplaceRenumber(map.reverseFaceMap(), masterFace_);
00792 
00793         // Stored faces refer to backed-up vertices (not changed)
00794         // and normal vertices (need to be renumbered)
00795         forAll(faceSetsVertices_, setI)
00796         {
00797             faceList& faces = faceSetsVertices_[setI];
00798 
00799             forAll(faces, i)
00800             {
00801                 // Note: faces[i] can have negative elements.
00802                 face& f = faces[i];
00803 
00804                 forAll(f, fp)
00805                 {
00806                     label pointI = f[fp];
00807 
00808                     if (pointI >= 0)
00809                     {
00810                         f[fp] = map.reversePointMap()[pointI];
00811 
00812                         if (f[fp] < 0)
00813                         {
00814                             FatalErrorIn
00815                             (
00816                                 "combineFaces::updateMesh"
00817                                 "(const mapPolyMesh&)"
00818                             )   << "In set " << setI << " at position " << i
00819                                 << " with master face "
00820                                 << masterFace_[setI] << nl
00821                                 << "the points of the slave face " << faces[i]
00822                                 << " don't exist anymore."
00823                                 << abort(FatalError);
00824                         }
00825                     }
00826                 }
00827             }
00828         }
00829     }
00830 }
00831 
00832 
00833 // Note that faces on coupled patches are never combined (or at least
00834 // getMergeSets never picks them up. Thus the points on them will never get
00835 // deleted since still used by the coupled faces. So never a risk of undoing
00836 // things (faces/points) on coupled patches.
00837 void Foam::combineFaces::setUnrefinement
00838 (
00839     const labelList& masterFaces,
00840     polyTopoChange& meshMod,
00841     Map<label>& restoredPoints,
00842     Map<label>& restoredFaces,
00843     Map<label>& restoredCells
00844 )
00845 {
00846     if (!undoable_)
00847     {
00848         FatalErrorIn
00849         (
00850             "combineFaces::setUnrefinement"
00851             "(const labelList&, polyTopoChange&"
00852             ", Map<label>&, Map<label>&, Map<label>&)"
00853         )   << "Can only call setUnrefinement if constructed with"
00854             << " unrefinement capability." << exit(FatalError);
00855     }
00856 
00857 
00858     // Restored points
00859     labelList addedPoints(savedPoints_.size(), -1);
00860 
00861     // Invert set-to-master-face
00862     Map<label> masterToSet(masterFace_.size());
00863 
00864     forAll(masterFace_, setI)
00865     {
00866         if (masterFace_[setI] >= 0)
00867         {
00868             masterToSet.insert(masterFace_[setI], setI);
00869         }
00870     }
00871 
00872     forAll(masterFaces, i)
00873     {
00874         label masterFaceI = masterFaces[i];
00875 
00876         Map<label>::const_iterator iter = masterToSet.find(masterFaceI);
00877 
00878         if (iter == masterToSet.end())
00879         {
00880             FatalErrorIn
00881             (
00882                 "combineFaces::setUnrefinement"
00883                 "(const labelList&, polyTopoChange&"
00884                 ", Map<label>&, Map<label>&, Map<label>&)"
00885             )   << "Master face " << masterFaceI
00886                 << " is not the master of one of the merge sets"
00887                 << " or has already been merged"
00888                 << abort(FatalError);
00889         }
00890 
00891         label setI = iter();
00892 
00893 
00894         // Update faces of the merge setI for reintroduced vertices
00895         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00896 
00897         faceList& faces = faceSetsVertices_[setI];
00898 
00899         if (faces.empty())
00900         {
00901             FatalErrorIn
00902             (
00903                 "combineFaces::setUnrefinement"
00904                 "(const labelList&, polyTopoChange&"
00905                 ", Map<label>&, Map<label>&, Map<label>&)"
00906             )   << "Set " << setI << " with master face " << masterFaceI
00907                 << " has already been merged." << abort(FatalError);
00908         }
00909 
00910         forAll(faces, i)
00911         {
00912             face& f = faces[i];
00913 
00914             forAll(f, fp)
00915             {
00916                 label pointI = f[fp];
00917 
00918                 if (pointI < 0)
00919                 {
00920                     label localI = -pointI-1;
00921 
00922                     if (addedPoints[localI] == -1)
00923                     {
00924                         // First occurrence of saved point. Reintroduce point
00925                         addedPoints[localI] = meshMod.setAction
00926                         (
00927                             polyAddPoint
00928                             (
00929                                 savedPoints_[localI],   // point
00930                                 -1,                     // master point
00931                                 -1,                     // zone for point
00932                                 true                    // supports a cell
00933                             )
00934                         );
00935                         restoredPoints.insert
00936                         (
00937                             addedPoints[localI],        // current point label
00938                             savedPointLabels_[localI]   // point label when it
00939                                                         // was stored
00940                         );
00941                     }
00942                     f[fp] = addedPoints[localI];
00943                 }
00944             }
00945         }
00946 
00947 
00948         // Restore
00949         // ~~~~~~~
00950 
00951         label own = mesh_.faceOwner()[masterFaceI];
00952         label zoneID = mesh_.faceZones().whichZone(masterFaceI);
00953         bool zoneFlip = false;
00954         if (zoneID >= 0)
00955         {
00956             const faceZone& fZone = mesh_.faceZones()[zoneID];
00957             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
00958         }
00959         label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
00960 
00961         if (mesh_.boundaryMesh()[patchI].coupled())
00962         {
00963             FatalErrorIn
00964             (
00965                 "combineFaces::setUnrefinement"
00966                 "(const labelList&, polyTopoChange&"
00967                 ", Map<label>&, Map<label>&, Map<label>&)"
00968             )   << "Master face " << masterFaceI << " is on coupled patch "
00969                 << mesh_.boundaryMesh()[patchI].name()
00970                 << abort(FatalError);
00971         }
00972 
00973         //Pout<< "Restoring new master face " << masterFaceI
00974         //    << " to vertices " << faces[0] << endl;
00975 
00976         // Modify the master face.
00977         meshMod.setAction
00978         (
00979             polyModifyFace
00980             (
00981                 faces[0],                       // original face
00982                 masterFaceI,                    // label of face
00983                 own,                            // owner
00984                 -1,                             // neighbour
00985                 false,                          // face flip
00986                 patchI,                         // patch for face
00987                 false,                          // remove from zone
00988                 zoneID,                         // zone for face
00989                 zoneFlip                        // face flip in zone
00990             )
00991         );
00992 
00993         // Add the previously removed faces
00994         for (label i = 1; i < faces.size(); i++)
00995         {
00996             //Pout<< "Restoring removed face with vertices " << faces[i]
00997             //    << endl;
00998 
00999             meshMod.setAction
01000             (
01001                 polyAddFace
01002                 (
01003                     faces[i],        // vertices
01004                     own,                    // owner,
01005                     -1,                     // neighbour,
01006                     -1,                     // masterPointID,
01007                     -1,                     // masterEdgeID,
01008                     masterFaceI,             // masterFaceID,
01009                     false,                  // flipFaceFlux,
01010                     patchI,                 // patchID,
01011                     zoneID,                 // zoneID,
01012                     zoneFlip                // zoneFlip
01013                 )
01014             );
01015         }
01016 
01017         // Clear out restored set
01018         faceSetsVertices_[setI].clear();
01019         masterFace_[setI] = -1;
01020     }
01021 }
01022 
01023 
01024 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines