cellDistFuncs.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) 2011-2012 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 "cellDistFuncs.H"
27 #include "polyMesh.H"
28 #include "wallPolyPatch.H"
29 #include "polyBoundaryMesh.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 defineTypeNameAndDebug(cellDistFuncs, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 // Find val in first nElems elements of elems.
42 Foam::label Foam::cellDistFuncs::findIndex
43 (
44  const label nElems,
45  const labelList& elems,
46  const label val
47 )
48 {
49  for (label i = 0; i < nElems; i++)
50  {
51  if (elems[i] == val)
52  {
53  return i;
54  }
55  }
56  return -1;
57 }
58 
59 
60 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
61 
62 Foam::cellDistFuncs::cellDistFuncs(const polyMesh& mesh)
63 :
64  mesh_(mesh)
65 {}
66 
67 
68 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
69 
71 (
72  const wordReList& patchNames
73 ) const
74 {
75  return mesh().boundaryMesh().patchSet(patchNames, false);
76 }
77 
78 
79 // Return smallest true distance from p to any of wallFaces.
80 // Note that even if normal hits face we still check other faces.
81 // Note that wallFaces is untruncated and we explicitly pass in size.
83 (
84  const point& p,
85  const polyPatch& patch,
86  const label nWallFaces,
87  const labelList& wallFaces,
88  label& minFaceI
89 ) const
90 {
91  const pointField& points = patch.points();
92 
93  scalar minDist = GREAT;
94  minFaceI = -1;
95 
96  for (label wallFaceI = 0; wallFaceI < nWallFaces; wallFaceI++)
97  {
98  label patchFaceI = wallFaces[wallFaceI];
99 
100  pointHit curHit = patch[patchFaceI].nearestPoint(p, points);
101 
102  if (curHit.distance() < minDist)
103  {
104  minDist = curHit.distance();
105  minFaceI = patch.start() + patchFaceI;
106  }
107  }
108 
109  return minDist;
110 }
111 
112 
113 // Get point neighbours of faceI (including faceI). Returns number of faces.
114 // Note: does not allocate storage but does use linear search to determine
115 // uniqueness. For polygonal faces this might be quite inefficient.
117 (
118  const primitivePatch& patch,
119  const label patchFaceI,
120  labelList& neighbours
121 ) const
122 {
123  label nNeighbours = 0;
124 
125  // Add myself
126  neighbours[nNeighbours++] = patchFaceI;
127 
128  // Add all face neighbours
129  const labelList& faceNeighbours = patch.faceFaces()[patchFaceI];
130 
131  forAll(faceNeighbours, faceNeighbourI)
132  {
133  neighbours[nNeighbours++] = faceNeighbours[faceNeighbourI];
134  }
135 
136  // Remember part of neighbours that contains edge-connected faces.
137  label nEdgeNbs = nNeighbours;
138 
139 
140  // Add all point-only neighbours by linear searching in edge neighbours.
141  // Assumes that point-only neighbours are not using multiple points on
142  // face.
143 
144  const face& f = patch.localFaces()[patchFaceI];
145 
146  forAll(f, fp)
147  {
148  label pointI = f[fp];
149 
150  const labelList& pointNbs = patch.pointFaces()[pointI];
151 
152  forAll(pointNbs, nbI)
153  {
154  label faceI = pointNbs[nbI];
155 
156  // Check for faceI in edge-neighbours part of neighbours
157  if (findIndex(nEdgeNbs, neighbours, faceI) == -1)
158  {
159  neighbours[nNeighbours++] = faceI;
160  }
161  }
162  }
163 
164 
165  if (debug)
166  {
167  // Check for duplicates
168 
169  // Use hashSet to determine nbs.
170  labelHashSet nbs(4*f.size());
171 
172  forAll(f, fp)
173  {
174  const labelList& pointNbs = patch.pointFaces()[f[fp]];
175 
176  forAll(pointNbs, i)
177  {
178  nbs.insert(pointNbs[i]);
179  }
180  }
181 
182  // Subtract ours.
183  for (label i = 0; i < nNeighbours; i++)
184  {
185  label nb = neighbours[i];
186 
187  if (!nbs.found(nb))
188  {
189  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
190  << "getPointNeighbours : patchFaceI:" << patchFaceI
191  << " verts:" << f << endl;
192 
193  forAll(f, fp)
194  {
195  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
196  << "point:" << f[fp] << " pointFaces:"
197  << patch.pointFaces()[f[fp]] << endl;
198  }
199 
200  for (label i = 0; i < nNeighbours; i++)
201  {
202  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
203  << "fast nbr:" << neighbours[i]
204  << endl;
205  }
206 
207  FatalErrorIn("getPointNeighbours")
208  << "Problem: fast pointNeighbours routine included " << nb
209  << " which is not in proper neigbour list " << nbs.toc()
210  << abort(FatalError);
211  }
212  nbs.erase(nb);
213  }
214 
215  if (nbs.size())
216  {
217  FatalErrorIn("getPointNeighbours")
218  << "Problem: fast pointNeighbours routine did not find "
219  << nbs.toc() << abort(FatalError);
220  }
221  }
222 
223  return nNeighbours;
224 }
225 
226 
227 // size of largest patch (out of supplied subset of patches)
229 (
230  const labelHashSet& patchIDs
231 ) const
232 {
233  label maxSize = 0;
234 
235  forAll(mesh().boundaryMesh(), patchI)
236  {
237  if (patchIDs.found(patchI))
238  {
239  const polyPatch& patch = mesh().boundaryMesh()[patchI];
240 
241  maxSize = Foam::max(maxSize, patch.size());
242  }
243  }
244  return maxSize;
245 }
246 
247 
248 // sum of patch sizes (out of supplied subset of patches)
250 (
251  const labelHashSet& patchIDs
252 )
253 const
254 {
255  label sum = 0;
256 
257  forAll(mesh().boundaryMesh(), patchI)
258  {
259  if (patchIDs.found(patchI))
260  {
261  const polyPatch& patch = mesh().boundaryMesh()[patchI];
262 
263  sum += patch.size();
264  }
265  }
266  return sum;
267 }
268 
269 
270 // Gets nearest wall for cells next to wall
272 (
273  const labelHashSet& patchIDs,
274  scalarField& wallDistCorrected,
275  Map<label>& nearestFace
276 ) const
277 {
278  // Size neighbours array for maximum possible (= size of largest patch)
279  label maxPointNeighbours = maxPatchSize(patchIDs);
280 
281  labelList neighbours(maxPointNeighbours);
282 
283 
284  // Correct all cells with face on wall
285  const vectorField& cellCentres = mesh().cellCentres();
286  const labelList& faceOwner = mesh().faceOwner();
287 
288  forAll(mesh().boundaryMesh(), patchI)
289  {
290  if (patchIDs.found(patchI))
291  {
292  const polyPatch& patch = mesh().boundaryMesh()[patchI];
293 
294  // Check cells with face on wall
295  forAll(patch, patchFaceI)
296  {
297  label nNeighbours = getPointNeighbours
298  (
299  patch,
300  patchFaceI,
301  neighbours
302  );
303 
304  label cellI = faceOwner[patch.start() + patchFaceI];
305 
306  label minFaceI = -1;
307 
308  wallDistCorrected[cellI] = smallestDist
309  (
310  cellCentres[cellI],
311  patch,
312  nNeighbours,
313  neighbours,
314  minFaceI
315  );
316 
317  // Store wallCell and its nearest neighbour
318  nearestFace.insert(cellI, minFaceI);
319  }
320  }
321  }
322 
323 }
324 
325 
326 // Correct all cells connected to wall (via point) and not in nearestFace
328 (
329  const labelHashSet& patchIDs,
330  scalarField& wallDistCorrected,
331  Map<label>& nearestFace
332 ) const
333 {
334  // Correct all (non-visited) cells with point on wall
335 
336  const vectorField& cellCentres = mesh().cellCentres();
337 
338  forAll(mesh().boundaryMesh(), patchI)
339  {
340  if (patchIDs.found(patchI))
341  {
342  const polyPatch& patch = mesh().boundaryMesh()[patchI];
343 
344  const labelList& meshPoints = patch.meshPoints();
345  const labelListList& pointFaces = patch.pointFaces();
346 
347  forAll(meshPoints, meshPointI)
348  {
349  label vertI = meshPoints[meshPointI];
350 
351  const labelList& neighbours = mesh().pointCells(vertI);
352 
353  forAll(neighbours, neighbourI)
354  {
355  label cellI = neighbours[neighbourI];
356 
357  if (!nearestFace.found(cellI))
358  {
359  const labelList& wallFaces = pointFaces[meshPointI];
360 
361  label minFaceI = -1;
362 
363  wallDistCorrected[cellI] = smallestDist
364  (
365  cellCentres[cellI],
366  patch,
367  wallFaces.size(),
368  wallFaces,
369  minFaceI
370  );
371 
372  // Store wallCell and its nearest neighbour
373  nearestFace.insert(cellI, minFaceI);
374  }
375  }
376  }
377  }
378  }
379 }
380 
381 
382 // ************************************************************************* //
#define SeriousErrorIn(functionName)
Report an error message using Foam::SeriousError.
const pointField & points
const labelListList & faceFaces() const
Return face-face addressing.
void correctBoundaryPointCells(const labelHashSet &patchIDs, scalarField &wallDistCorrected, Map< label > &nearestFace) const
Correct all cells connected to wall (via point). Sets values in.
label nWallFaces(0)
labelList f(nPoints)
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
void correctBoundaryFaceCells(const labelHashSet &patchIDs, scalarField &wallDistCorrected, Map< label > &nearestFace) const
Correct all cells connected to boundary (via face). Sets values in.
const labelListList & pointFaces() const
Return point-face addressing.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Definition: boundaryMesh.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
label getPointNeighbours(const primitivePatch &, const label patchFaceI, labelList &) const
Get faces sharing point with face on patch.
const labelListList & pointCells() const
const polyMesh & mesh() const
Access mesh.
Definition: cellDistFuncs.H:95
const vectorField & cellCentres() const
dynamicFvMesh & mesh
Namespace for OpenFOAM.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
scalar smallestDist(const point &p, const polyPatch &patch, const label nWallFaces, const labelList &wallFaces, label &meshFaceI) const
Calculate smallest true distance (and face index)
Definition: cellDistFuncs.C:83
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
#define forAll(list, i)
Definition: UList.H:421
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A list of faces which address into the list of points.
labelHashSet patchSet(const UList< wordRe > &patchNames, const bool warnNotFound=true, const bool usePatchGroups=true) const
Return the set of patch IDs corresponding to the given names.
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
label sumPatchSize(const labelHashSet &patchIDs) const
Sum of patch sizes (out of supplied subset of patches).
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
const Field< PointType > & points() const
Return reference to global points.
labelHashSet getPatchIDs() const
Get patchIDs of/derived off certain type (e.g. &#39;processorPolyPatch&#39;)
defineTypeNameAndDebug(combustionModel, 0)
label maxPatchSize(const labelHashSet &patchIDs) const
Size of largest patch (out of supplied subset of patches)