FreeFOAM The Cross-Platform CFD Toolkit
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) 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 "cellDistFuncs.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/wallPolyPatch.H>
30 
31 namespace Foam
32 {
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 defineTypeNameAndDebug(cellDistFuncs, 0);
37 
38 }
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 // Find val in first nElems elements of elems.
44 Foam::label Foam::cellDistFuncs::findIndex
45 (
46  const label nElems,
47  const labelList& elems,
48  const label val
49 )
50 {
51  for (label i = 0; i < nElems; i++)
52  {
53  if (elems[i] == val)
54  {
55  return i;
56  }
57  }
58  return -1;
59 }
60 
61 
62 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
63 
64 // Construct from mesh
65 Foam::cellDistFuncs::cellDistFuncs(const polyMesh& mesh)
66 :
67  mesh_(mesh)
68 {}
69 
70 
71 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
72 
73 // Get patch ids of named patches
75 (
76  const wordList& patchNames
77 ) const
78 {
80 
81  // Construct Set of patchNames for easy checking if included
82  HashSet<word> patchNamesHash(patchNames.size());
83 
84  forAll(patchNames, patchI)
85  {
86  patchNamesHash.insert(patchNames[patchI]);
87  }
88 
89  // Loop over all patches and check if patch name in hashset
90 
91  labelHashSet patchIDs(bMesh.size());
92 
93  forAll(bMesh, patchI)
94  {
95  const polyPatch& patch = bMesh[patchI];
96 
97  if (patchNamesHash.found(patch.name()))
98  {
99  patchIDs.insert(patchI);
100  }
101  }
102  return patchIDs;
103 }
104 
105 
106 // Return smallest true distance from p to any of wallFaces.
107 // Note that even if normal hits face we still check other faces.
108 // Note that wallFaces is untruncated and we explicitly pass in size.
110 (
111  const point& p,
112  const polyPatch& patch,
113  const label nWallFaces,
114  const labelList& wallFaces,
115  label& minFaceI
116 ) const
117 {
118  const pointField& points = patch.points();
119 
120  scalar minDist = GREAT;
121  minFaceI = -1;
122 
123  for(label wallFaceI = 0; wallFaceI < nWallFaces; wallFaceI++)
124  {
125  label patchFaceI = wallFaces[wallFaceI];
126 
127  pointHit curHit = patch[patchFaceI].nearestPoint(p, points);
128 
129  if (curHit.distance() < minDist)
130  {
131  minDist = curHit.distance();
132  minFaceI = patch.start() + patchFaceI;
133  }
134  }
135 
136  return minDist;
137 }
138 
139 
140 // Get point neighbours of faceI (including faceI). Returns number of faces.
141 // Note: does not allocate storage but does use linear search to determine
142 // uniqueness. For polygonal faces this might be quite inefficient.
144 (
145  const primitivePatch& patch,
146  const label patchFaceI,
147  labelList& neighbours
148 ) const
149 {
150  label nNeighbours = 0;
151 
152  // Add myself
153  neighbours[nNeighbours++] = patchFaceI;
154 
155  // Add all face neighbours
156  const labelList& faceNeighbours = patch.faceFaces()[patchFaceI];
157 
158  forAll(faceNeighbours, faceNeighbourI)
159  {
160  neighbours[nNeighbours++] = faceNeighbours[faceNeighbourI];
161  }
162 
163  // Remember part of neighbours that contains edge-connected faces.
164  label nEdgeNbs = nNeighbours;
165 
166 
167  // Add all point-only neighbours by linear searching in edge neighbours.
168  // Assumes that point-only neighbours are not using multiple points on
169  // face.
170 
171  const face& f = patch.localFaces()[patchFaceI];
172 
173  forAll(f, fp)
174  {
175  label pointI = f[fp];
176 
177  const labelList& pointNbs = patch.pointFaces()[pointI];
178 
179  forAll(pointNbs, nbI)
180  {
181  label faceI = pointNbs[nbI];
182 
183  // Check for faceI in edge-neighbours part of neighbours
184  if (findIndex(nEdgeNbs, neighbours, faceI) == -1)
185  {
186  neighbours[nNeighbours++] = faceI;
187  }
188  }
189  }
190 
191 
192  if (debug)
193  {
194  // Check for duplicates
195 
196  // Use hashSet to determine nbs.
197  labelHashSet nbs(4*f.size());
198 
199  forAll(f, fp)
200  {
201  const labelList& pointNbs = patch.pointFaces()[f[fp]];
202 
203  forAll(pointNbs, i)
204  {
205  nbs.insert(pointNbs[i]);
206  }
207  }
208 
209  // Subtract ours.
210  for (label i = 0; i < nNeighbours; i++)
211  {
212  label nb = neighbours[i];
213 
214  if (!nbs.found(nb))
215  {
216  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
217  << "getPointNeighbours : patchFaceI:" << patchFaceI
218  << " verts:" << f << endl;
219 
220  forAll(f, fp)
221  {
222  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
223  << "point:" << f[fp] << " pointFaces:"
224  << patch.pointFaces()[f[fp]] << endl;
225  }
226 
227  for (label i = 0; i < nNeighbours; i++)
228  {
229  SeriousErrorIn("Foam::cellDistFuncs::getPointNeighbours")
230  << "fast nbr:" << neighbours[i]
231  << endl;
232  }
233 
234  FatalErrorIn("getPointNeighbours")
235  << "Problem: fast pointNeighbours routine included " << nb
236  << " which is not in proper neigbour list " << nbs.toc()
237  << abort(FatalError);
238  }
239  nbs.erase(nb);
240  }
241 
242  if (nbs.size())
243  {
244  FatalErrorIn("getPointNeighbours")
245  << "Problem: fast pointNeighbours routine did not find "
246  << nbs.toc() << abort(FatalError);
247  }
248  }
249 
250  return nNeighbours;
251 }
252 
253 
254 // size of largest patch (out of supplied subset of patches)
256  const
257 {
258  label maxSize = 0;
259 
260  forAll(mesh().boundaryMesh(), patchI)
261  {
262  if (patchIDs.found(patchI))
263  {
264  const polyPatch& patch = mesh().boundaryMesh()[patchI];
265 
266  maxSize = Foam::max(maxSize, patch.size());
267  }
268  }
269  return maxSize;
270 }
271 
272 
273 // sum of patch sizes (out of supplied subset of patches)
275  const
276 {
277  label sum = 0;
278 
279  forAll(mesh().boundaryMesh(), patchI)
280  {
281  if (patchIDs.found(patchI))
282  {
283  const polyPatch& patch = mesh().boundaryMesh()[patchI];
284 
285  sum += patch.size();
286  }
287  }
288  return sum;
289 }
290 
291 
292 // Gets nearest wall for cells next to wall
294 (
295  const labelHashSet& patchIDs,
296  scalarField& wallDistCorrected,
297  Map<label>& nearestFace
298 ) const
299 {
300  // Size neighbours array for maximum possible (= size of largest patch)
301  label maxPointNeighbours = maxPatchSize(patchIDs);
302 
303  labelList neighbours(maxPointNeighbours);
304 
305 
306  // Correct all cells with face on wall
307  const vectorField& cellCentres = mesh().cellCentres();
308  const labelList& faceOwner = mesh().faceOwner();
309 
310  forAll(mesh().boundaryMesh(), patchI)
311  {
312  if (patchIDs.found(patchI))
313  {
314  const polyPatch& patch = mesh().boundaryMesh()[patchI];
315 
316  // Check cells with face on wall
317  forAll(patch, patchFaceI)
318  {
319  label nNeighbours = getPointNeighbours
320  (
321  patch,
322  patchFaceI,
323  neighbours
324  );
325 
326  label cellI = faceOwner[patch.start() + patchFaceI];
327 
328  label minFaceI = -1;
329 
330  wallDistCorrected[cellI] = smallestDist
331  (
332  cellCentres[cellI],
333  patch,
334  nNeighbours,
335  neighbours,
336  minFaceI
337  );
338 
339  // Store wallCell and its nearest neighbour
340  nearestFace.insert(cellI, minFaceI);
341  }
342  }
343  }
344 
345 }
346 
347 
348 // Correct all cells connected to wall (via point) and not in nearestFace
350 (
351  const labelHashSet& patchIDs,
352  scalarField& wallDistCorrected,
353  Map<label>& nearestFace
354 ) const
355 {
356  // Correct all (non-visited) cells with point on wall
357 
358  const labelListList& pointCells = mesh().pointCells();
359  const vectorField& cellCentres = mesh().cellCentres();
360 
361  forAll(mesh().boundaryMesh(), patchI)
362  {
363  if (patchIDs.found(patchI))
364  {
365  const polyPatch& patch = mesh().boundaryMesh()[patchI];
366 
367  const labelList& meshPoints = patch.meshPoints();
368  const labelListList& pointFaces = patch.pointFaces();
369 
370  forAll(meshPoints, meshPointI)
371  {
372  label vertI = meshPoints[meshPointI];
373 
374  const labelList& neighbours = pointCells[vertI];
375 
376  forAll(neighbours, neighbourI)
377  {
378  label cellI = neighbours[neighbourI];
379 
380  if (!nearestFace.found(cellI))
381  {
382  const labelList& wallFaces = pointFaces[meshPointI];
383 
384  label minFaceI = -1;
385 
386  wallDistCorrected[cellI] = smallestDist
387  (
388  cellCentres[cellI],
389  patch,
390  wallFaces.size(),
391  wallFaces,
392  minFaceI
393  );
394 
395  // Store wallCell and its nearest neighbour
396  nearestFace.insert(cellI, minFaceI);
397  }
398  }
399  }
400  }
401  }
402 }
403 
404 
405 // ************************ vim: set sw=4 sts=4 et: ************************ //