src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "primitiveMesh.H"
  27 #include "cell.H"
  28 #include "boundBox.H"
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 bool Foam::primitiveMesh::pointInCellBB
  33 (
  34     const point& p,
  35     label celli,
  36     scalar tol
  37 ) const
  38 {
  39     boundBox bb
  40     (
  41         cells()[celli].points
  42         (
  43             faces(),
  44             points()
  45         ),
  46         false
  47     );
  48
  49     if (tol > SMALL)
  50     {
  51         bb = boundBox
  52         (
  53             bb.min() - tol*bb.span(),
  54             bb.max() + tol*bb.span()
  55         );
  56     }
  57
  58     return bb.contains(p);
  59 }
  60
  61
  62 bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
  63 {
  64     const labelList& f = cells()[celli];
  65     const labelList& owner = this->faceOwner();
  66     const vectorField& cf = faceCentres();
  67     const vectorField& Sf = faceAreas();
  68
  69     bool inCell = true;
  70
  71     forAll(f, facei)
  72     {
  73         label nFace = f[facei];
  74         vector proj = p - cf[nFace];
  75         vector normal = Sf[nFace];
  76         if (owner[nFace] != celli)
  77         {
  78             normal = -normal;
  79         }
  80         inCell = inCell && ((normal & proj) <= 0);
  81     }
  82
  83     return inCell;
  84 }
  85
  86
  87 Foam::label Foam::primitiveMesh::findNearestCell(const point& location) const
  88 {
  89     const vectorField& centres = cellCentres();
  90
  91     label nearestCelli = 0;
  92     scalar minProximity = magSqr(centres[0] - location);
  93
  94     for (label celli = 1; celli < centres.size(); celli++)
  95     {
  96         scalar proximity = magSqr(centres[celli] - location);
  97
  98         if (proximity < minProximity)
  99         {
 100             nearestCelli = celli;
 101             minProximity = proximity;
 102         }
 103     }
 104
 105     return nearestCelli;
 106 }
 107
 108
 109 Foam::label Foam::primitiveMesh::findCell(const point& location) const
 110 {
 111     if (nCells() == 0)
 112     {
 113         return -1;
 114     }
 115
 116     // Find the nearest cell centre to this location
 117     label celli = findNearestCell(location);
 118
 119     // If point is in the nearest cell return
 120     if (pointInCell(location, celli))
 121     {
 122         return celli;
 123     }
 124     else // point is not in the nearest cell so search all cells
 125     {
 126         bool cellFound = false;
 127         label n = 0;
 128
 129         while ((!cellFound) && (n < nCells()))
 130         {
 131             if (pointInCell(location, n))
 132             {
 133                 cellFound = true;
 134                 celli = n;
 135             }
 136             else
 137             {
 138                 n++;
 139             }
 140         }
 141         if (cellFound)
 142         {
 143             return celli;
 144         }
 145         else
 146         {
 147             return -1;
 148         }
 149     }
 150 }
 151
 152
 153 // ************************************************************************* //