FreeFOAM The Cross-Platform CFD Toolkit
regionSide.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 "regionSide.H"
27 #include <meshTools/meshTools.H>
28 #include <OpenFOAM/primitiveMesh.H>
29 
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 defineTypeNameAndDebug(regionSide, 0);
37 
38 }
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 // Step across edge onto other face on cell
44 (
45  const primitiveMesh& mesh,
46  const label cellI,
47  const label faceI,
48  const label edgeI
49 )
50 {
51  label f0I, f1I;
52 
53  meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
54 
55  if (f0I == faceI)
56  {
57  return f1I;
58  }
59  else
60  {
61  return f0I;
62  }
63 }
64 
65 
66 // Step across point to other edge on face
68 (
69  const primitiveMesh& mesh,
70  const label faceI,
71  const label edgeI,
72  const label pointI
73 )
74 {
75  const edge& e = mesh.edges()[edgeI];
76 
77  // Get other point on edge.
78  label freePointI = e.otherVertex(pointI);
79 
80  const labelList& fEdges = mesh.faceEdges()[faceI];
81 
82  forAll(fEdges, fEdgeI)
83  {
84  label otherEdgeI = fEdges[fEdgeI];
85 
86  const edge& otherE = mesh.edges()[otherEdgeI];
87 
88  if
89  (
90  (
91  otherE.start() == pointI
92  && otherE.end() != freePointI
93  )
94  || (
95  otherE.end() == pointI
96  && otherE.start() != freePointI
97  )
98  )
99  {
100  // otherE shares one (but not two) points with e.
101  return otherEdgeI;
102  }
103  }
104 
106  (
107  "regionSide::otherEdge(const primitiveMesh&, const label, const label"
108  ", const label)"
109  ) << "Cannot find other edge on face " << faceI << " that uses point "
110  << pointI << " but not point " << freePointI << endl
111  << "Edges on face:" << fEdges
112  << " verts:" << UIndirectList<edge>(mesh.edges(), fEdges)()
113  << " Vertices on face:"
114  << mesh.faces()[faceI]
115  << " Vertices on original edge:" << e << abort(FatalError);
116 
117  return -1;
118 }
119 
120 
121 // Step from faceI (on side cellI) to connected face & cell without crossing
122 // fenceEdges.
123 void Foam::regionSide::visitConnectedFaces
124 (
125  const primitiveMesh& mesh,
126  const labelHashSet& region,
127  const labelHashSet& fenceEdges,
128  const label cellI,
129  const label faceI,
130  labelHashSet& visitedFace
131 )
132 {
133  if (!visitedFace.found(faceI))
134  {
135  if (debug)
136  {
137  Info<< "visitConnectedFaces : cellI:" << cellI << " faceI:"
138  << faceI << " isOwner:" << (cellI == mesh.faceOwner()[faceI])
139  << endl;
140  }
141 
142  // Mark as visited
143  visitedFace.insert(faceI);
144 
145  // Mark which side of face was visited.
146  if (cellI == mesh.faceOwner()[faceI])
147  {
148  sideOwner_.insert(faceI);
149  }
150 
151 
152  // Visit all neighbouring faces on faceSet. Stay on this 'side' of
153  // face by doing edge-face-cell walk.
154  const labelList& fEdges = mesh.faceEdges()[faceI];
155 
156  forAll(fEdges, fEdgeI)
157  {
158  label edgeI = fEdges[fEdgeI];
159 
160  if (!fenceEdges.found(edgeI))
161  {
162  // Step along faces on edge from cell to cell until
163  // we hit face on faceSet.
164 
165  // Find face reachable from edge
166  label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
167 
168  if (mesh.isInternalFace(otherFaceI))
169  {
170  label otherCellI = cellI;
171 
172  // Keep on crossing faces/cells until back on face on
173  // surface
174  while (!region.found(otherFaceI))
175  {
176  visitedFace.insert(otherFaceI);
177 
178  if (debug)
179  {
180  Info<< "visitConnectedFaces : cellI:" << cellI
181  << " found insideEdgeFace:" << otherFaceI
182  << endl;
183  }
184 
185 
186  // Cross otherFaceI into neighbouring cell
187  otherCellI =
189  (
190  mesh,
191  otherCellI,
192  otherFaceI
193  );
194 
195  otherFaceI =
196  otherFace
197  (
198  mesh,
199  otherCellI,
200  otherFaceI,
201  edgeI
202  );
203  }
204 
205  visitConnectedFaces
206  (
207  mesh,
208  region,
209  fenceEdges,
210  otherCellI,
211  otherFaceI,
212  visitedFace
213  );
214  }
215  }
216  }
217  }
218 }
219 
220 
221 // From edge on face connected to point on region (regionPointI) cross
222 // to all other edges using this point by walking across faces
223 // Does not cross regionEdges so stays on one side
224 // of region
225 void Foam::regionSide::walkPointConnectedFaces
226 (
227  const primitiveMesh& mesh,
228  const labelHashSet& regionEdges,
229  const label regionPointI,
230  const label startFaceI,
231  const label startEdgeI,
232  labelHashSet& visitedEdges
233 )
234 {
235  // Mark as visited
236  insidePointFaces_.insert(startFaceI);
237 
238  if (debug)
239  {
240  Info<< "walkPointConnectedFaces : regionPointI:" << regionPointI
241  << " faceI:" << startFaceI
242  << " edgeI:" << startEdgeI << " verts:"
243  << mesh.edges()[startEdgeI]
244  << endl;
245  }
246 
247  // Cross faceI i.e. get edge not startEdgeI which uses regionPointI
248  label edgeI = otherEdge(mesh, startFaceI, startEdgeI, regionPointI);
249 
250  if (!regionEdges.found(edgeI))
251  {
252  if (!visitedEdges.found(edgeI))
253  {
254  visitedEdges.insert(edgeI);
255 
256  if (debug)
257  {
258  Info<< "Crossed face from "
259  << " edgeI:" << startEdgeI << " verts:"
260  << mesh.edges()[startEdgeI]
261  << " to edge:" << edgeI << " verts:"
262  << mesh.edges()[edgeI]
263  << endl;
264  }
265 
266  // Cross edge to all faces connected to it.
267 
268  const labelList& eFaces = mesh.edgeFaces()[edgeI];
269 
270  forAll(eFaces, eFaceI)
271  {
272  label faceI = eFaces[eFaceI];
273 
274  walkPointConnectedFaces
275  (
276  mesh,
277  regionEdges,
278  regionPointI,
279  faceI,
280  edgeI,
281  visitedEdges
282  );
283  }
284  }
285  }
286 }
287 
288 
289 // Find all faces reachable from all non-fence points and staying on
290 // regionFaces side.
291 void Foam::regionSide::walkAllPointConnectedFaces
292 (
293  const primitiveMesh& mesh,
294  const labelHashSet& regionFaces,
295  const labelHashSet& fencePoints
296 )
297 {
298  //
299  // Get all (internal and external) edges on region.
300  //
301  labelHashSet regionEdges(4*regionFaces.size());
302 
303  for
304  (
305  labelHashSet::const_iterator iter = regionFaces.begin();
306  iter != regionFaces.end();
307  ++iter
308  )
309  {
310  label faceI = iter.key();
311 
312  const labelList& fEdges = mesh.faceEdges()[faceI];
313 
314  forAll(fEdges, fEdgeI)
315  {
316  regionEdges.insert(fEdges[fEdgeI]);
317  }
318  }
319 
320 
321  //
322  // Visit all internal points on surface.
323  //
324 
325  // Storage for visited points
326  labelHashSet visitedPoint(4*regionFaces.size());
327 
328  // Insert fence points so we don't visit them
329  for
330  (
331  labelHashSet::const_iterator iter = fencePoints.begin();
332  iter != fencePoints.end();
333  ++iter
334  )
335  {
336  visitedPoint.insert(iter.key());
337  }
338 
339  labelHashSet visitedEdges(2*fencePoints.size());
340 
341 
342  if (debug)
343  {
344  Info<< "Excluding visit of points:" << visitedPoint << endl;
345  }
346 
347  for
348  (
349  labelHashSet::const_iterator iter = regionFaces.begin();
350  iter != regionFaces.end();
351  ++iter
352  )
353  {
354  label faceI = iter.key();
355 
356  // Get side of face.
357  label cellI;
358 
359  if (sideOwner_.found(faceI))
360  {
361  cellI = mesh.faceOwner()[faceI];
362  }
363  else
364  {
365  cellI = mesh.faceNeighbour()[faceI];
366  }
367 
368  // Find starting point and edge on face.
369  const labelList& fEdges = mesh.faceEdges()[faceI];
370 
371  forAll(fEdges, fEdgeI)
372  {
373  label edgeI = fEdges[fEdgeI];
374 
375  // Get the face 'perpendicular' to faceI on region.
376  label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
377 
378  // Edge
379  const edge& e = mesh.edges()[edgeI];
380 
381  if (!visitedPoint.found(e.start()))
382  {
383  Info<< "Determining visibility from point " << e.start()
384  << endl;
385 
386  visitedPoint.insert(e.start());
387 
388  //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.start());
389 
390  walkPointConnectedFaces
391  (
392  mesh,
393  regionEdges,
394  e.start(),
395  otherFaceI,
396  edgeI,
397  visitedEdges
398  );
399  }
400  if (!visitedPoint.found(e.end()))
401  {
402  Info<< "Determining visibility from point " << e.end()
403  << endl;
404 
405  visitedPoint.insert(e.end());
406 
407  //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.end());
408 
409  walkPointConnectedFaces
410  (
411  mesh,
412  regionEdges,
413  e.end(),
414  otherFaceI,
415  edgeI,
416  visitedEdges
417  );
418  }
419  }
420  }
421 }
422 
423 
424 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
425 
426 // Construct from components
428 (
429  const primitiveMesh& mesh,
430  const labelHashSet& region, // faces of region
431  const labelHashSet& fenceEdges, // outside edges
432  const label startCellI,
433  const label startFaceI
434 )
435 :
436  sideOwner_(region.size()),
437  insidePointFaces_(region.size())
438 {
439  // Storage for visited faces
440  labelHashSet visitedFace(region.size());
441 
442  // Visit all faces on this side of region.
443  // Mark for each face whether owner (or neighbour) has been visited.
444  // Sets sideOwner_
445  visitConnectedFaces
446  (
447  mesh,
448  region,
449  fenceEdges,
450  startCellI,
451  startFaceI,
452  visitedFace
453  );
454 
455  //
456  // Visit all internal points on region and mark faces visitable from these.
457  // Sets insidePointFaces_.
458  //
459 
460  labelHashSet fencePoints(fenceEdges.size());
461 
462  for
463  (
464  labelHashSet::const_iterator iter = fenceEdges.begin();
465  iter != fenceEdges.end();
466  ++iter
467  )
468  {
469  const edge& e = mesh.edges()[iter.key()];
470 
471  fencePoints.insert(e.start());
472  fencePoints.insert(e.end());
473  }
474 
475  walkAllPointConnectedFaces(mesh, region, fencePoints);
476 }
477 
478 
479 // ************************ vim: set sw=4 sts=4 et: ************************ //