src/meshTools/cellDist/cellDistFuncs.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "cellDistFuncs.H"
  28 #include "polyMesh.H"
  29 #include "wallPolyPatch.H"
  30 #include "polyBoundaryMesh.H"
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(cellDistFuncs, 0);
  38
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 // Find val in first nElems elements of elems.
  45 Foam::label Foam::cellDistFuncs::findIndex
  46 (
  47     const label nElems,
  48     const labelList& elems,
  49     const label val
  50 )
  51 {
  52     for (label i = 0; i < nElems; i++)
  53     {
  54         if (elems[i] == val)
  55         {
  56             return i;
  57         }
  58     }
  59     return -1;
  60 }
  61
  62
  63 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  64
  65 // Construct from mesh
  66 Foam::cellDistFuncs::cellDistFuncs(const polyMesh& mesh)
  67 :
  68     mesh_(mesh)
  69 {}
  70
  71
  72 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  73
  74 // Get patch ids of named patches
  75 Foam::labelHashSet Foam::cellDistFuncs::getPatchIDs
  76 (
  77     const wordList& patchNames
  78 ) const
  79 {
  80     const polyBoundaryMesh& bMesh = mesh().boundaryMesh();
  81
  82     // Construct Set of patchNames for easy checking if included
  83     HashSet<word> patchNamesHash(patchNames.size());
  84
  85     forAll(patchNames, patchI)
  86     {
  87         patchNamesHash.insert(patchNames[patchI]);
  88     }
  89
  90     // Loop over all patches and check if patch name in hashset
  91
  92     labelHashSet patchIDs(bMesh.size());
  93
  94     forAll(bMesh, patchI)
  95     {
  96         const polyPatch& patch = bMesh[patchI];
  97
  98         if (patchNamesHash.found(patch.name()))
  99         {
 100             patchIDs.insert(patchI);
 101         }
 102     }
 103     return patchIDs;
 104 }
 105
 106
 107 // Return smallest true distance from p to any of wallFaces.
 108 // Note that even if normal hits face we still check other faces.
 109 // Note that wallFaces is untruncated and we explicitly pass in size.
 110 Foam::scalar Foam::cellDistFuncs::smallestDist
 111 (
 112     const point& p,
 113     const polyPatch& patch,
 114     const label nWallFaces,
 115     const labelList& wallFaces,
 116     label& minFaceI
 117 ) const
 118 {
 119     const pointField& points = patch.points();
 120
 121     scalar minDist = GREAT;
 122     minFaceI = -1;
 123
 124     for(label wallFaceI = 0; wallFaceI < nWallFaces; wallFaceI++)
 125     {
 126         label patchFaceI = wallFaces[wallFaceI];
 127
 128         pointHit curHit = patch[patchFaceI].nearestPoint(p, points);
 129
 130         if (curHit.distance() < minDist)
 131         {
 132             minDist = curHit.distance();
 133             minFaceI = patch.start() + patchFaceI;
 134         }
 135     }
 136
 137     return minDist;
 138 }
 139
 140
 141 // Get point neighbours of faceI (including faceI). Returns number of faces.
 142 // Note: does not allocate storage but does use linear search to determine
 143 // uniqueness. For polygonal faces this might be quite inefficient.
 144 Foam::label Foam::cellDistFuncs::getPointNeighbours
 145 (
 146     const primitivePatch& patch,
 147     const label patchFaceI,
 148     labelList& neighbours
 149 ) const
 150 {
 151     label nNeighbours = 0;
 152
 153     // Add myself
 154     neighbours[nNeighbours++] = patchFaceI;
 155
 156     // Add all face neighbours
 157     const labelList& faceNeighbours = patch.faceFaces()[patchFaceI];
 158
 159     forAll(faceNeighbours, faceNeighbourI)
 160     {
 161         neighbours[nNeighbours++] = faceNeighbours[faceNeighbourI];
 162     }
 163
 164     // Remember part of neighbours that contains edge-connected faces.
 165     label nEdgeNbs = nNeighbours;
 166
 167
 168     // Add all point-only neighbours by linear searching in edge neighbours.
 169     // Assumes that point-only neighbours are not using multiple points on
 170     // face.
 171
 172     const face& f = patch.localFaces()[patchFaceI];
 173
 174     forAll(f, fp)
 175     {
 176         label pointI = f[fp];
 177
 178         const labelList& pointNbs = patch.pointFaces()[pointI];
 179
 180         forAll(pointNbs, nbI)
 181         {
 182             label faceI = pointNbs[nbI];
 183
 184             // Check for faceI in edge-neighbours part of neighbours
 185             if (findIndex(nEdgeNbs, neighbours, faceI) == -1)
 186             {
 187                 neighbours[nNeighbours++] = faceI;
 188             }
 189         }
 190     }
 191
 192
 193     if (debug)
 194     {
 195         // Check for duplicates
 196
 197         // Use hashSet to determine nbs.
 198         labelHashSet nbs(4*f.size());
 199
 200         forAll(f, fp)
 201         {
 202             const labelList& pointNbs = patch.pointFaces()[f[fp]];
 203
 204             forAll(pointNbs, i)
 205             {
 206                 nbs.insert(pointNbs[i]);
 207             }
 208         }
 209
 210         // Subtract ours.
 211         for (label i = 0; i < nNeighbours; i++)
 212         {
 213             label nb = neighbours[i];
 214
 215             if (!nbs.found(nb))
 216             {
 217                 SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
 218                     << "getPointNeighbours : patchFaceI:" << patchFaceI
 219                     << " verts:" << f << endl;
 220
 221                 forAll(f, fp)
 222                 {
 223                     SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
 224                         << "point:" << f[fp] << " pointFaces:"
 225                         << patch.pointFaces()[f[fp]] << endl;
 226                 }
 227
 228                 for (label i = 0; i < nNeighbours; i++)
 229                 {
 230                     SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
 231                         << "fast nbr:" << neighbours[i]
 232                         << endl;
 233                 }
 234
 235                 FatalErrorIn("getPointNeighbours")
 236                     << "Problem: fast pointNeighbours routine included " << nb
 237                     << " which is not in proper neigbour list " << nbs.toc()
 238                     << abort(FatalError);
 239             }
 240             nbs.erase(nb);
 241         }
 242
 243         if (nbs.size())
 244         {
 245             FatalErrorIn("getPointNeighbours")
 246                 << "Problem: fast pointNeighbours routine did not find "
 247                 << nbs.toc() << abort(FatalError);
 248         }
 249     }
 250
 251     return nNeighbours;
 252 }
 253
 254
 255 // size of largest patch (out of supplied subset of patches)
 256 Foam::label Foam::cellDistFuncs::maxPatchSize(const labelHashSet& patchIDs)
 257  const
 258 {
 259     label maxSize = 0;
 260
 261     forAll(mesh().boundaryMesh(), patchI)
 262     {
 263         if (patchIDs.found(patchI))
 264         {
 265             const polyPatch& patch = mesh().boundaryMesh()[patchI];
 266
 267             maxSize = Foam::max(maxSize, patch.size());
 268         }
 269     }
 270     return maxSize;
 271 }
 272
 273
 274 // sum of patch sizes (out of supplied subset of patches)
 275 Foam::label Foam::cellDistFuncs::sumPatchSize(const labelHashSet& patchIDs)
 276  const
 277 {
 278     label sum = 0;
 279
 280     forAll(mesh().boundaryMesh(), patchI)
 281     {
 282         if (patchIDs.found(patchI))
 283         {
 284             const polyPatch& patch = mesh().boundaryMesh()[patchI];
 285
 286             sum += patch.size();
 287         }
 288     }
 289     return sum;
 290 }
 291
 292
 293 // Gets nearest wall for cells next to wall
 294 void Foam::cellDistFuncs::correctBoundaryFaceCells
 295 (
 296     const labelHashSet& patchIDs,
 297     scalarField& wallDistCorrected,
 298     Map<label>& nearestFace
 299 ) const
 300 {
 301     // Size neighbours array for maximum possible (= size of largest patch)
 302     label maxPointNeighbours = maxPatchSize(patchIDs);
 303
 304     labelList neighbours(maxPointNeighbours);
 305
 306
 307     // Correct all cells with face on wall
 308     const vectorField& cellCentres = mesh().cellCentres();
 309     const labelList& faceOwner = mesh().faceOwner();
 310
 311     forAll(mesh().boundaryMesh(), patchI)
 312     {
 313         if (patchIDs.found(patchI))
 314         {
 315             const polyPatch& patch = mesh().boundaryMesh()[patchI];
 316
 317             // Check cells with face on wall
 318             forAll(patch, patchFaceI)
 319             {
 320                 label nNeighbours = getPointNeighbours
 321                 (
 322                     patch,
 323                     patchFaceI,
 324                     neighbours
 325                 );
 326
 327                 label cellI = faceOwner[patch.start() + patchFaceI];
 328
 329                 label minFaceI = -1;
 330
 331                 wallDistCorrected[cellI] = smallestDist
 332                 (
 333                     cellCentres[cellI],
 334                     patch,
 335                     nNeighbours,
 336                     neighbours,
 337                     minFaceI
 338                 );
 339
 340                 // Store wallCell and its nearest neighbour
 341                 nearestFace.insert(cellI, minFaceI);
 342             }
 343         }
 344     }
 345
 346 }
 347
 348
 349 // Correct all cells connected to wall (via point) and not in nearestFace
 350 void Foam::cellDistFuncs::correctBoundaryPointCells
 351 (
 352     const labelHashSet& patchIDs,
 353     scalarField& wallDistCorrected,
 354     Map<label>& nearestFace
 355 ) const
 356 {
 357     // Correct all (non-visited) cells with point on wall
 358
 359     const labelListList& pointCells = mesh().pointCells();
 360     const vectorField& cellCentres = mesh().cellCentres();
 361
 362     forAll(mesh().boundaryMesh(), patchI)
 363     {
 364         if (patchIDs.found(patchI))
 365         {
 366             const polyPatch& patch = mesh().boundaryMesh()[patchI];
 367
 368             const labelList& meshPoints = patch.meshPoints();
 369             const labelListList& pointFaces = patch.pointFaces();
 370
 371             forAll(meshPoints, meshPointI)
 372             {
 373                 label vertI = meshPoints[meshPointI];
 374
 375                 const labelList& neighbours = pointCells[vertI];
 376
 377                 forAll(neighbours, neighbourI)
 378                 {
 379                     label cellI = neighbours[neighbourI];
 380
 381                     if (!nearestFace.found(cellI))
 382                     {
 383                         const labelList& wallFaces = pointFaces[meshPointI];
 384
 385                         label minFaceI = -1;
 386
 387                         wallDistCorrected[cellI] = smallestDist
 388                         (
 389                             cellCentres[cellI],
 390                             patch,
 391                             wallFaces.size(),
 392                             wallFaces,
 393                             minFaceI
 394                         );
 395
 396                         // Store wallCell and its nearest neighbour
 397                         nearestFace.insert(cellI, minFaceI);
 398                     }
 399                 }
 400             }
 401         }
 402     }
 403 }
 404
 405
 406 // ************************************************************************* //