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