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