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