src/meshTools/cellDist/cellDistFuncs.C

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