src/meshTools/cellDist/cellDistFuncs.C

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