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

localPointRegion.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 "localPointRegion.H"
00027 #include <OpenFOAM/syncTools.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mapPolyMesh.H>
00030 #include <OpenFOAM/globalIndex.H>
00031 #include <OpenFOAM/indirectPrimitivePatch.H>
00032 
00033 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00034 
00035 namespace Foam
00036 {
00037 
00038 defineTypeNameAndDebug(localPointRegion, 0);
00039 
00040 // Reduction class to get minimum value over face.
00041 class minEqOpFace
00042 {
00043 public:
00044 
00045     void operator()(face& x, const face& y) const
00046     {
00047         if (x.size())
00048         {
00049             label j = 0;
00050             forAll(x, i)
00051             {
00052                 x[i] = min(x[i], y[j]);
00053 
00054                 j = y.rcIndex(j);
00055             }
00056         }
00057     };
00058 };
00059 
00060 
00061 // Dummy transform for faces. Used in synchronisation
00062 void transformList
00063 (
00064     const tensorField& rotTensor,
00065     UList<face>& field
00066 )
00067 {};
00068 
00069 }
00070 
00071 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00072 
00073 // Are two lists identical either in forward or in reverse order.
00074 bool Foam::localPointRegion::isDuplicate
00075 (
00076     const face& f0,
00077     const face& f1,
00078     const bool forward
00079 )
00080 {
00081     label fp1 = findIndex(f1, f0[0]);
00082 
00083     if (fp1 == -1)
00084     {
00085         return false;
00086     }
00087 
00088     forAll(f0, fp0)
00089     {
00090         if (f0[fp0] != f1[fp1])
00091         {
00092             return false;
00093         }
00094 
00095         if (forward)
00096         {
00097             fp1 = f1.fcIndex(fp1);
00098         }
00099         else
00100         {
00101             fp1 = f1.rcIndex(fp1);
00102         }
00103     }
00104     return true;
00105 }
00106 
00107 
00108 // Count regions per point
00109 void Foam::localPointRegion::countPointRegions
00110 (
00111     const polyMesh& mesh,
00112     const boolList& candidatePoint,
00113     const Map<label>& candidateFace,
00114     faceList& minRegion
00115 )
00116 {
00117     // Almost all will have only one so only
00118     // populate Map if more than one.
00119     labelList minPointRegion(mesh.nPoints(), -1);
00120     // From global point to local (multi-region) point numbering
00121     meshPointMap_.resize(candidateFace.size()/100);
00122     // From local (multi-region) point to regions
00123     DynamicList<labelList> pointRegions(meshPointMap_.size());
00124 
00125     // From faces with any duplicated point on it to local face
00126     meshFaceMap_.resize(meshPointMap_.size());
00127 
00128     forAllConstIter(Map<label>, candidateFace, iter)
00129     {
00130         label faceI = iter.key();
00131 
00132         if (!mesh.isInternalFace(faceI))
00133         {
00134             const face& f = mesh.faces()[faceI];
00135 
00136             if (minRegion[faceI].empty())
00137             {
00138                 FatalErrorIn("localPointRegion::countPointRegions(..)")
00139                     << "Face from candidateFace without minRegion set." << endl
00140                     << "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
00141                     << " verts:" << f << abort(FatalError);
00142             }
00143 
00144             forAll(f, fp)
00145             {
00146                 label pointI = f[fp];
00147 
00148                 // Even points which were not candidates for splitting might
00149                 // be on multiple baffles that are being split so check.
00150 
00151                 if (candidatePoint[pointI])
00152                 {
00153                     label region = minRegion[faceI][fp];
00154 
00155                     if (minPointRegion[pointI] == -1)
00156                     {
00157                         minPointRegion[pointI] = region;
00158                     }
00159                     else if (minPointRegion[pointI] != region)
00160                     {
00161                         // Multiple regions for this point. Add.
00162                         Map<label>::iterator iter = meshPointMap_.find(pointI);
00163                         if (iter != meshPointMap_.end())
00164                         {
00165                             labelList& regions = pointRegions[iter()];
00166                             if (findIndex(regions, region) == -1)
00167                             {
00168                                 label sz = regions.size();
00169                                 regions.setSize(sz+1);
00170                                 regions[sz] = region;
00171                             }
00172                         }
00173                         else
00174                         {
00175                             label localPointI = meshPointMap_.size();
00176                             meshPointMap_.insert(pointI, localPointI);
00177                             labelList regions(2);
00178                             regions[0] = minPointRegion[pointI];
00179                             regions[1] = region;
00180                             pointRegions.append(regions);
00181                         }
00182 
00183                         label meshFaceMapI = meshFaceMap_.size();
00184                         meshFaceMap_.insert(faceI, meshFaceMapI);
00185                     }
00186                 }
00187             }
00188         }
00189     }
00190     minPointRegion.clear();
00191 
00192     // Add internal faces that use any duplicated point. Can only have one
00193     // region!
00194     forAllConstIter(Map<label>, candidateFace, iter)
00195     {
00196         label faceI = iter.key();
00197 
00198         if (mesh.isInternalFace(faceI))
00199         {
00200             const face& f = mesh.faces()[faceI];
00201 
00202             forAll(f, fp)
00203             {
00204                 // Note: candidatePoint test not really necessary but
00205                 // speeds up rejection.
00206                 if (candidatePoint[f[fp]] && meshPointMap_.found(f[fp]))
00207                 {
00208                     label meshFaceMapI = meshFaceMap_.size();
00209                     meshFaceMap_.insert(faceI, meshFaceMapI);
00210                 }
00211             }
00212         }
00213     }
00214 
00215 
00216     // Transfer to member data
00217     pointRegions.shrink();
00218     pointRegions_.setSize(pointRegions.size());
00219     forAll(pointRegions, i)
00220     {
00221         pointRegions_[i].transfer(pointRegions[i]);
00222     }
00223 
00224     // Compact minRegion
00225     faceRegions_.setSize(meshFaceMap_.size());
00226     forAllConstIter(Map<label>, meshFaceMap_, iter)
00227     {
00228         faceRegions_[iter()].labelList::transfer(minRegion[iter.key()]);
00229 
00231         //{
00232         //    label faceI = iter.key();
00233         //    const face& f = mesh.faces()[faceI];
00234         //    Pout<< "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
00235         //        << " verts:" << f << endl;
00236         //    forAll(f, fp)
00237         //    {
00238         //        Pout<< "    " << f[fp] << " min:" << faceRegions_[iter()][fp]
00239         //            << endl;
00240         //    }
00241         //    Pout<< endl;
00242         //}
00243     }
00244 
00245     // Compact region numbering
00246     // ? TBD.
00247 }
00248 
00249 
00250 void Foam::localPointRegion::calcPointRegions
00251 (
00252     const polyMesh& mesh,
00253     boolList& candidatePoint
00254 )
00255 {
00256     label nBnd = mesh.nFaces()-mesh.nInternalFaces();
00257     const labelList& faceOwner = mesh.faceOwner();
00258     const labelList& faceNeighbour = mesh.faceNeighbour();
00259 
00260 
00261     syncTools::syncPointList
00262     (
00263         mesh,
00264         candidatePoint,
00265         orEqOp<bool>(),
00266         false,              // nullValue
00267         false               // applySeparation
00268     );
00269 
00270 
00271     // Mark any face/boundaryFace/cell with a point on a candidate point.
00272     // - candidateFace does not necessary have to be a baffle!
00273     // - candidateFace is synchronised (since candidatePoint is)
00274     Map<label> candidateFace(2*nBnd);
00275     label candidateFaceI = 0;
00276 
00277     Map<label> candidateCell(nBnd);
00278     label candidateCellI = 0;
00279 
00280     forAll(mesh.faces(), faceI)
00281     {
00282         const face& f = mesh.faces()[faceI];
00283 
00284         forAll(f, fp)
00285         {
00286             if (candidatePoint[f[fp]])
00287             {
00288                 // Mark face
00289                 if (candidateFace.insert(faceI, candidateFaceI))
00290                 {
00291                     candidateFaceI++;
00292                 }
00293 
00294                 // Mark cells
00295                 if (candidateCell.insert(faceOwner[faceI], candidateCellI))
00296                 {
00297                     candidateCellI++;
00298                 }
00299 
00300                 if (mesh.isInternalFace(faceI))
00301                 {
00302                     label nei = faceNeighbour[faceI];
00303                     if (candidateCell.insert(nei, candidateCellI))
00304                     {
00305                         candidateCellI++;
00306                     }
00307                 }
00308 
00309                 break;
00310             }
00311         }
00312     }
00313 
00314 
00315 
00316     // Get global indices for cells
00317     globalIndex globalCells(mesh.nCells());
00318 
00319 
00320     // Determine for every candidate face per point the minimum region
00321     // (global cell) it is connected to. (candidateFaces are the
00322     // only ones using a
00323     // candidate point so the only ones that can be affected)
00324     faceList minRegion(mesh.nFaces());
00325     forAllConstIter(Map<label>, candidateFace, iter)
00326     {
00327         label faceI = iter.key();
00328         const face& f = mesh.faces()[faceI];
00329 
00330         if (mesh.isInternalFace(faceI))
00331         {
00332             label globOwn = globalCells.toGlobal(faceOwner[faceI]);
00333             label globNei = globalCells.toGlobal(faceNeighbour[faceI]);
00334             minRegion[faceI].setSize(f.size(), min(globOwn, globNei));
00335         }
00336         else
00337         {
00338             label globOwn = globalCells.toGlobal(faceOwner[faceI]);
00339             minRegion[faceI].setSize(f.size(), globOwn);
00340         }
00341     }
00342 
00343     // Now minimize over all faces that are connected through internal
00344     // faces or through cells. This loop iterates over the max number of
00345     // cells connected to a point (=8 for hex mesh)
00346 
00347     while (true)
00348     {
00349         // Transport minimum from face across cell
00350         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00351 
00352         Map<label> minPointValue(100);
00353         label nChanged = 0;
00354         forAllConstIter(Map<label>, candidateCell, iter)
00355         {
00356             minPointValue.clear();
00357 
00358             label cellI = iter.key();
00359             const cell& cFaces = mesh.cells()[cellI];
00360 
00361             // Determine minimum per point
00362             forAll(cFaces, cFaceI)
00363             {
00364                 label faceI = cFaces[cFaceI];
00365 
00366                 if (minRegion[faceI].size())
00367                 {
00368                     const face& f = mesh.faces()[faceI];
00369 
00370                     forAll(f, fp)
00371                     {
00372                         label pointI = f[fp];
00373                         Map<label>::iterator iter = minPointValue.find(pointI);
00374 
00375                         if (iter == minPointValue.end())
00376                         {
00377                             minPointValue.insert(pointI, minRegion[faceI][fp]);
00378                         }
00379                         else
00380                         {
00381                             label currentMin = iter();
00382                             iter() = min(currentMin, minRegion[faceI][fp]);
00383                         }
00384                     }
00385                 }
00386             }
00387 
00388             // Set face minimum from point minimum
00389             forAll(cFaces, cFaceI)
00390             {
00391                 label faceI = cFaces[cFaceI];
00392 
00393                 if (minRegion[faceI].size())
00394                 {
00395                     const face& f = mesh.faces()[faceI];
00396 
00397                     forAll(f, fp)
00398                     {
00399                         label minVal = minPointValue[f[fp]];
00400 
00401                         if (minVal != minRegion[faceI][fp])
00402                         {
00403                             minRegion[faceI][fp] = minVal;
00404                             nChanged++;
00405                         }
00406                     }
00407                 }
00408             }
00409         }
00410 
00411         //Pout<< "nChanged:" << nChanged << endl;
00412 
00413         if (returnReduce(nChanged, sumOp<label>()) == 0)
00414         {
00415             break;
00416         }
00417 
00418 
00419         // Transport minimum across coupled faces
00420         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00421 
00422         syncTools::syncFaceList
00423         (
00424             mesh,
00425             minRegion,
00426             minEqOpFace(),
00427             false               // applySeparation
00428         );
00429     }
00430 
00431 
00432     // Count regions per point
00433     countPointRegions(mesh, candidatePoint, candidateFace, minRegion);
00434     minRegion.clear();
00435 
00436 
00438     //forAllConstIter(Map<label>, meshPointMap_, iter)
00439     //{
00440     //    Pout<< "point:" << iter.key()
00441     //        << " coord:" << mesh.points()[iter.key()]
00442     //        << " regions:" << pointRegions_[iter()] << endl;
00443     //}
00444 }
00445 
00446 
00447 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00448 
00449 Foam::localPointRegion::localPointRegion(const polyMesh& mesh)
00450 :
00451     //nRegions_(0),
00452     meshPointMap_(0),
00453     pointRegions_(0),
00454     meshFaceMap_(0),
00455     faceRegions_(0)
00456 {
00457     const polyBoundaryMesh& patches = mesh.boundaryMesh();
00458 
00459     // Get any point on the outside which is on a non-coupled boundary
00460     boolList candidatePoint(mesh.nPoints(), false);
00461 
00462     forAll(patches, patchI)
00463     {
00464         if (!patches[patchI].coupled())
00465         {
00466             const polyPatch& pp = patches[patchI];
00467 
00468             forAll(pp.meshPoints(), i)
00469             {
00470                 candidatePoint[pp.meshPoints()[i]] = true;
00471             }
00472         }
00473     }
00474 
00475     calcPointRegions(mesh, candidatePoint);
00476 }
00477 
00478 
00479 Foam::localPointRegion::localPointRegion
00480 (
00481     const polyMesh& mesh,
00482     const labelList& candidatePoints
00483 )
00484 :
00485     //nRegions_(0),
00486     meshPointMap_(0),
00487     pointRegions_(0),
00488     meshFaceMap_(0),
00489     faceRegions_(0)
00490 {
00491     boolList candidatePoint(mesh.nPoints(), false);
00492 
00493     forAll(candidatePoints, i)
00494     {
00495         candidatePoint[candidatePoints[i]] = true;
00496     }
00497 
00498     calcPointRegions(mesh, candidatePoint);
00499 }
00500 
00501 
00502 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00503 
00504 // Return a list (in allPatch indices) with either -1 or the face label
00505 // of the face that uses the same vertices.
00506 Foam::labelList Foam::localPointRegion::findDuplicateFaces
00507 (
00508     const primitiveMesh& mesh,
00509     const labelList& boundaryFaces
00510 )
00511 {
00512     // Addressing engine for all boundary faces.
00513     indirectPrimitivePatch allPatch
00514     (
00515         IndirectList<face>(mesh.faces(), boundaryFaces),
00516         mesh.points()
00517     );
00518 
00519     labelList duplicateFace(allPatch.size(), -1);
00520     label nDuplicateFaces = 0;
00521 
00522     // Find all duplicate faces.
00523     forAll(allPatch, bFaceI)
00524     {
00525         const face& f = allPatch.localFaces()[bFaceI];
00526 
00527         // Get faces connected to f[0].
00528         // Check whether share all points with f.
00529         const labelList& pFaces = allPatch.pointFaces()[f[0]];
00530 
00531         forAll(pFaces, i)
00532         {
00533             label otherFaceI = pFaces[i];
00534 
00535             if (otherFaceI > bFaceI)
00536             {
00537                 const face& otherF = allPatch.localFaces()[otherFaceI];
00538 
00539                 if (isDuplicate(f, otherF, true))
00540                 {
00541                     FatalErrorIn
00542                     (
00543                         "findDuplicateFaces(const primitiveMesh&"
00544                         ", const labelList&)"
00545                     )   << "Face:" << bFaceI + mesh.nInternalFaces()
00546                         << " has local points:" << f
00547                         << " which are in same order as face:"
00548                         << otherFaceI + mesh.nInternalFaces()
00549                         << " with local points:" << otherF
00550                         << abort(FatalError);
00551                 }
00552                 else if (isDuplicate(f, otherF, false))
00553                 {
00554                     label meshFace0 = bFaceI + mesh.nInternalFaces();
00555                     label meshFace1 = otherFaceI + mesh.nInternalFaces();
00556 
00557                     if
00558                     (
00559                         duplicateFace[bFaceI] != -1
00560                      || duplicateFace[otherFaceI] != -1
00561                     )
00562                     {
00563                         FatalErrorIn
00564                         (
00565                             "findDuplicateFaces(const primitiveMesh&"
00566                             ", const labelList&)"
00567                         )   << "One of two duplicate faces already marked"
00568                             << " as duplicate." << nl
00569                             << "This means that three or more faces share"
00570                             << " the same points and this is illegal." << nl
00571                             << "Face:" << meshFace0
00572                             << " with local points:" << f
00573                             << " which are in same order as face:"
00574                             << meshFace1
00575                             << " with local points:" << otherF
00576                             << abort(FatalError);
00577                     }
00578 
00579                     duplicateFace[bFaceI] = otherFaceI;
00580                     duplicateFace[otherFaceI] = bFaceI;
00581                     nDuplicateFaces++;
00582                 }
00583             }
00584         }
00585     }
00586 
00587     return duplicateFace;
00588 }
00589 
00590 
00591 void Foam::localPointRegion::updateMesh(const mapPolyMesh& map)
00592 {
00593     {
00594         Map<label> newMap(meshFaceMap_.size());
00595 
00596         forAllConstIter(Map<label>, meshFaceMap_, iter)
00597         {
00598             label newFaceI = map.reverseFaceMap()[iter.key()];
00599 
00600             if (newFaceI >= 0)
00601             {
00602                 newMap.insert(newFaceI, iter());
00603             }
00604         }
00605         meshFaceMap_.transfer(newMap);
00606     }
00607     {
00608         Map<label> newMap(meshPointMap_.size());
00609 
00610         forAllConstIter(Map<label>, meshPointMap_, iter)
00611         {
00612             label newPointI = map.reversePointMap()[iter.key()];
00613 
00614             if (newPointI >= 0)
00615             {
00616                 newMap.insert(newPointI, iter());
00617             }
00618         }
00619 
00620         meshPointMap_.transfer(newMap);
00621     }
00622 }
00623 
00624 
00625 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines