cellDistFuncs.C 10 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9 10 11 12
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

13 14 15 16
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
17 18 19 20 21 22 23

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
24
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
25 26 27 28 29 30 31 32

\*---------------------------------------------------------------------------*/

#include "cellDistFuncs.H"
#include "polyMesh.H"
#include "wallPolyPatch.H"
#include "polyBoundaryMesh.H"

33
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34

35 36 37 38
namespace Foam
{
defineTypeNameAndDebug(cellDistFuncs, 0);
}
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

// Find val in first nElems elements of elems.
Foam::label Foam::cellDistFuncs::findIndex
(
    const label nElems,
    const labelList& elems,
    const label val
)
{
    for (label i = 0; i < nElems; i++)
    {
        if (elems[i] == val)
        {
            return i;
        }
    }
    return -1;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::cellDistFuncs::cellDistFuncs(const polyMesh& mesh)
:
    mesh_(mesh)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::labelHashSet Foam::cellDistFuncs::getPatchIDs
(
74
    const UList<wordRe>& patchNames
75 76
) const
{
77
    return mesh().boundaryMesh().patchSet(patchNames, false);
78 79 80 81 82 83 84 85 86 87 88 89
}


// Return smallest true distance from p to any of wallFaces.
// Note that even if normal hits face we still check other faces.
// Note that wallFaces is untruncated and we explicitly pass in size.
Foam::scalar Foam::cellDistFuncs::smallestDist
(
    const point& p,
    const polyPatch& patch,
    const label nWallFaces,
    const labelList& wallFaces,
90
    label& minFacei
91 92 93 94 95
) const
{
    const pointField& points = patch.points();

    scalar minDist = GREAT;
96
    minFacei = -1;
97

98
    for (label wallFacei = 0; wallFacei < nWallFaces; wallFacei++)
99
    {
100
        label patchFacei = wallFaces[wallFacei];
101

102
        pointHit curHit = patch[patchFacei].nearestPoint(p, points);
103 104 105 106

        if (curHit.distance() < minDist)
        {
            minDist = curHit.distance();
107
            minFacei = patch.start() + patchFacei;
108 109 110 111 112 113 114
        }
    }

    return minDist;
}


115
// Get point neighbours of facei (including facei). Returns number of faces.
116 117 118 119 120
// Note: does not allocate storage but does use linear search to determine
// uniqueness. For polygonal faces this might be quite inefficient.
Foam::label Foam::cellDistFuncs::getPointNeighbours
(
    const primitivePatch& patch,
121
    const label patchFacei,
122 123 124 125 126 127
    labelList& neighbours
) const
{
    label nNeighbours = 0;

    // Add myself
128
    neighbours[nNeighbours++] = patchFacei;
129 130

    // Add all face neighbours
131
    const labelList& faceNeighbours = patch.faceFaces()[patchFacei];
132 133 134 135 136 137 138 139 140 141 142 143 144 145

    forAll(faceNeighbours, faceNeighbourI)
    {
        neighbours[nNeighbours++] = faceNeighbours[faceNeighbourI];
    }

    // Remember part of neighbours that contains edge-connected faces.
    label nEdgeNbs = nNeighbours;


    // Add all point-only neighbours by linear searching in edge neighbours.
    // Assumes that point-only neighbours are not using multiple points on
    // face.

146
    const face& f = patch.localFaces()[patchFacei];
147 148 149

    forAll(f, fp)
    {
150
        label pointi = f[fp];
151

152
        const labelList& pointNbs = patch.pointFaces()[pointi];
153 154 155

        forAll(pointNbs, nbI)
        {
156
            label facei = pointNbs[nbI];
157

158 159
            // Check for facei in edge-neighbours part of neighbours
            if (findIndex(nEdgeNbs, neighbours, facei) == -1)
160
            {
161
                neighbours[nNeighbours++] = facei;
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
            }
        }
    }


    if (debug)
    {
        // Check for duplicates

        // Use hashSet to determine nbs.
        labelHashSet nbs(4*f.size());

        forAll(f, fp)
        {
            const labelList& pointNbs = patch.pointFaces()[f[fp]];
177
            nbs.insert(pointNbs);
178 179 180 181 182 183 184 185 186
        }

        // Subtract ours.
        for (label i = 0; i < nNeighbours; i++)
        {
            label nb = neighbours[i];

            if (!nbs.found(nb))
            {
187
                SeriousErrorInFunction
188
                    << "getPointNeighbours : patchFacei:" << patchFacei
189 190 191 192
                    << " verts:" << f << endl;

                forAll(f, fp)
                {
193
                    SeriousErrorInFunction
194 195 196 197 198 199
                        << "point:" << f[fp] << " pointFaces:"
                        << patch.pointFaces()[f[fp]] << endl;
                }

                for (label i = 0; i < nNeighbours; i++)
                {
200
                    SeriousErrorInFunction
201 202 203 204
                        << "fast nbr:" << neighbours[i]
                        << endl;
                }

205
                FatalErrorInFunction
206
                    << "Problem: fast pointNeighbours routine included " << nb
Andrew Heather's avatar
Andrew Heather committed
207
                    << " which is not in proper neighbour list " << nbs.toc()
208 209 210 211 212
                    << abort(FatalError);
            }
            nbs.erase(nb);
        }

213
        if (nbs.size())
214
        {
215
            FatalErrorInFunction
216 217 218 219 220 221 222 223 224 225
                << "Problem: fast pointNeighbours routine did not find "
                << nbs.toc() << abort(FatalError);
        }
    }

    return nNeighbours;
}


// size of largest patch (out of supplied subset of patches)
226 227 228 229
Foam::label Foam::cellDistFuncs::maxPatchSize
(
    const labelHashSet& patchIDs
) const
230 231 232
{
    label maxSize = 0;

233
    forAll(mesh().boundaryMesh(), patchi)
234
    {
235
        if (patchIDs.found(patchi))
236
        {
237
            const polyPatch& patch = mesh().boundaryMesh()[patchi];
238 239 240 241 242 243 244 245 246

            maxSize = Foam::max(maxSize, patch.size());
        }
    }
    return maxSize;
}


// sum of patch sizes (out of supplied subset of patches)
247 248 249 250 251
Foam::label Foam::cellDistFuncs::sumPatchSize
(
    const labelHashSet& patchIDs
)
const
252 253 254
{
    label sum = 0;

255
    forAll(mesh().boundaryMesh(), patchi)
256
    {
257
        if (patchIDs.found(patchi))
258
        {
259
            const polyPatch& patch = mesh().boundaryMesh()[patchi];
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285

            sum += patch.size();
        }
    }
    return sum;
}


// Gets nearest wall for cells next to wall
void Foam::cellDistFuncs::correctBoundaryFaceCells
(
    const labelHashSet& patchIDs,
    scalarField& wallDistCorrected,
    Map<label>& nearestFace
) const
{
    // Size neighbours array for maximum possible (= size of largest patch)
    label maxPointNeighbours = maxPatchSize(patchIDs);

    labelList neighbours(maxPointNeighbours);


    // Correct all cells with face on wall
    const vectorField& cellCentres = mesh().cellCentres();
    const labelList& faceOwner = mesh().faceOwner();

286
    forAll(mesh().boundaryMesh(), patchi)
287
    {
288
        if (patchIDs.found(patchi))
289
        {
290
            const polyPatch& patch = mesh().boundaryMesh()[patchi];
291 292

            // Check cells with face on wall
293
            forAll(patch, patchFacei)
294 295 296 297
            {
                label nNeighbours = getPointNeighbours
                (
                    patch,
298
                    patchFacei,
299 300 301
                    neighbours
                );

302
                label celli = faceOwner[patch.start() + patchFacei];
303

304
                label minFacei = -1;
305

306
                wallDistCorrected[celli] = smallestDist
307
                (
308
                    cellCentres[celli],
309 310 311
                    patch,
                    nNeighbours,
                    neighbours,
312
                    minFacei
313 314 315
                );

                // Store wallCell and its nearest neighbour
316
                nearestFace.insert(celli, minFacei);
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
            }
        }
    }

}


// Correct all cells connected to wall (via point) and not in nearestFace
void Foam::cellDistFuncs::correctBoundaryPointCells
(
    const labelHashSet& patchIDs,
    scalarField& wallDistCorrected,
    Map<label>& nearestFace
) const
{
    // Correct all (non-visited) cells with point on wall

    const vectorField& cellCentres = mesh().cellCentres();

336
    forAll(mesh().boundaryMesh(), patchi)
337
    {
338
        if (patchIDs.found(patchi))
339
        {
340
            const polyPatch& patch = mesh().boundaryMesh()[patchi];
341 342 343 344

            const labelList& meshPoints = patch.meshPoints();
            const labelListList& pointFaces = patch.pointFaces();

345
            forAll(meshPoints, meshPointi)
346
            {
347
                label vertI = meshPoints[meshPointi];
348

349
                const labelList& neighbours = mesh().pointCells(vertI);
350 351 352

                forAll(neighbours, neighbourI)
                {
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 365 366
                            patch,
                            wallFaces.size(),
                            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
                    }
                }
            }
        }
    }
}


// ************************************************************************* //