src/foam/meshes/primitiveMesh/primitiveMeshFindCell.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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 // Is the point in the cell bounding box
  33 bool Foam::primitiveMesh::pointInCellBB(const point& p, label celli) const
  34 {
  35     // Make bounding box for check.  All points are local, so no reduce is
  36     // needed
  37     boundBox bb(cells()[celli].points(faces(), points()), false);
  38
  39     return bb.contains(p);
  40 }
  41
  42
  43 // Is the point in the cell
  44 bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const
  45 {
  46     const labelList& f = cells()[celli];
  47     const labelList& owner = this->faceOwner();
  48     const vectorField& cf = faceCentres();
  49     const vectorField& Sf = faceAreas();
  50
  51     bool inCell = true;
  52
  53     forAll(f, facei)
  54     {
  55         label nFace = f[facei];
  56         vector proj = p - cf[nFace];
  57         vector normal = Sf[nFace];
  58         if (owner[nFace] != celli)
  59         {
  60             normal = -normal;
  61         }
  62         inCell = inCell && ((normal & proj) <= 0);
  63     }
  64
  65     return inCell;
  66 }
  67
  68
  69 // Find the cell with the nearest cell centre
  70 Foam::label Foam::primitiveMesh::findNearestCell(const point& location) const
  71 {
  72     const vectorField& centres = cellCentres();
  73
  74     label nearestCelli = 0;
  75     scalar minProximity = magSqr(centres[0] - location);
  76
  77     for (label celli = 1; celli < centres.size(); celli++)
  78     {
  79         scalar proximity = magSqr(centres[celli] - location);
  80
  81         if (proximity < minProximity)
  82         {
  83             nearestCelli = celli;
  84             minProximity = proximity;
  85         }
  86     }
  87
  88     return nearestCelli;
  89 }
  90
  91
  92 // Find cell enclosing this location
  93 Foam::label Foam::primitiveMesh::findCell(const point& location) const
  94 {
  95     if (nCells() == 0)
  96     {
  97         return -1;
  98     }
  99
 100     // Find the nearest cell centre to this location
 101     label celli = findNearestCell(location);
 102
 103     // If point is in the nearest cell return
 104     if (pointInCell(location, celli))
 105     {
 106         return celli;
 107     }
 108     else // point is not in the nearest cell so search all cells
 109     {
 110         bool cellFound = false;
 111         label n = 0;
 112
 113         while ((!cellFound) && (n < nCells()))
 114         {
 115             if (pointInCell(location, n))
 116             {
 117                 cellFound = true;
 118                 celli = n;
 119             }
 120             else
 121             {
 122                 n++;
 123             }
 124         }
 125         if (cellFound)
 126         {
 127             return celli;
 128         }
 129         else
 130         {
 131             return -1;
 132         }
 133     }
 134 }
 135
 136
 137 // ************************************************************************* //