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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2018 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "regionSide.H"
30#include "meshTools.H"
31#include "primitiveMesh.H"
32
33
34// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35
36namespace Foam
37{
38
39defineTypeNameAndDebug(regionSide, 0);
40
41}
42
43// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44
45// Step across edge onto other face on cell
47(
48 const primitiveMesh& mesh,
49 const label celli,
50 const label facei,
51 const label edgeI
52)
53{
54 label f0I, f1I;
55
56 meshTools::getEdgeFaces(mesh, celli, edgeI, f0I, f1I);
57
58 if (f0I == facei)
59 {
60 return f1I;
61 }
62 else
63 {
64 return f0I;
65 }
66}
67
68
69// Step across point to other edge on face
70Foam::label Foam::regionSide::otherEdge
71(
72 const primitiveMesh& mesh,
73 const label facei,
74 const label edgeI,
75 const label pointi
76)
77{
78 const edge& e = mesh.edges()[edgeI];
79
80 // Get other point on edge.
81 label freePointi = e.otherVertex(pointi);
82
83 const labelList& fEdges = mesh.faceEdges()[facei];
84
85 forAll(fEdges, fEdgeI)
86 {
87 const label otherEdgeI = fEdges[fEdgeI];
88 const edge& otherE = mesh.edges()[otherEdgeI];
89
90 if
91 (
92 (
93 otherE.start() == pointi
94 && otherE.end() != freePointi
95 )
96 || (
97 otherE.end() == pointi
98 && otherE.start() != freePointi
99 )
100 )
101 {
102 // otherE shares one (but not two) points with e.
103 return otherEdgeI;
104 }
105 }
106
108 << "Cannot find other edge on face " << facei << " that uses point "
109 << pointi << " but not point " << freePointi << endl
110 << "Edges on face:" << fEdges
111 << " verts:" << UIndirectList<edge>(mesh.edges(), fEdges)
112 << " Vertices on face:"
113 << mesh.faces()[facei]
114 << " Vertices on original edge:" << e << abort(FatalError);
115
116 return -1;
117}
118
119
120// Step from facei (on side celli) to connected face & cell without crossing
121// fenceEdges.
122void Foam::regionSide::visitConnectedFaces
123(
124 const primitiveMesh& mesh,
125 const labelHashSet& region,
126 const labelHashSet& fenceEdges,
127 const label celli,
128 const label facei,
129 labelHashSet& visitedFace
130)
131{
132 if (!visitedFace.found(facei))
133 {
134 if (debug)
135 {
136 Info<< "visitConnectedFaces : celli:" << celli << " facei:"
137 << facei << " isOwner:" << (celli == mesh.faceOwner()[facei])
138 << endl;
139 }
140
141 // Mark as visited
142 visitedFace.insert(facei);
143
144 // Mark which side of face was visited.
145 if (celli == mesh.faceOwner()[facei])
146 {
147 sideOwner_.insert(facei);
148 }
149
150
151 // Visit all neighbouring faces on faceSet. Stay on this 'side' of
152 // face by doing edge-face-cell walk.
153 const labelList& fEdges = mesh.faceEdges()[facei];
154
155 forAll(fEdges, fEdgeI)
156 {
157 label edgeI = fEdges[fEdgeI];
158
159 if (!fenceEdges.found(edgeI))
160 {
161 // Step along faces on edge from cell to cell until
162 // we hit face on faceSet.
163
164 // Find face reachable from edge
165 label otherFacei = otherFace(mesh, celli, facei, edgeI);
166
167 if (mesh.isInternalFace(otherFacei))
168 {
169 label otherCelli = celli;
170
171 // Keep on crossing faces/cells until back on face on
172 // surface
173 while (!region.found(otherFacei))
174 {
175 visitedFace.insert(otherFacei);
176
177 if (debug)
178 {
179 Info<< "visitConnectedFaces : celli:" << celli
180 << " found insideEdgeFace:" << otherFacei
181 << endl;
182 }
183
184
185 // Cross otherFacei into neighbouring cell
186 otherCelli =
188 (
189 mesh,
190 otherCelli,
191 otherFacei
192 );
193
194 otherFacei =
196 (
197 mesh,
198 otherCelli,
199 otherFacei,
200 edgeI
201 );
202 }
203
204 visitConnectedFaces
205 (
206 mesh,
207 region,
208 fenceEdges,
209 otherCelli,
210 otherFacei,
211 visitedFace
212 );
213 }
214 }
215 }
216 }
217}
218
219
220// From edge on face connected to point on region (regionPointi) cross
221// to all other edges using this point by walking across faces
222// Does not cross regionEdges so stays on one side
223// of region
224void Foam::regionSide::walkPointConnectedFaces
225(
226 const primitiveMesh& mesh,
227 const labelHashSet& regionEdges,
228 const label regionPointi,
229 const label startFacei,
230 const label startEdgeI,
231 labelHashSet& visitedEdges
232)
233{
234 // Mark as visited
235 insidePointFaces_.insert(startFacei);
236
237 if (debug)
238 {
239 Info<< "walkPointConnectedFaces : regionPointi:" << regionPointi
240 << " facei:" << startFacei
241 << " edgeI:" << startEdgeI << " verts:"
242 << mesh.edges()[startEdgeI]
243 << endl;
244 }
245
246 // Cross facei i.e. get edge not startEdgeI which uses regionPointi
247 label edgeI = otherEdge(mesh, startFacei, startEdgeI, regionPointi);
248
249 if (!regionEdges.found(edgeI))
250 {
251 if (!visitedEdges.found(edgeI))
252 {
253 visitedEdges.insert(edgeI);
254
255 if (debug)
256 {
257 Info<< "Crossed face from "
258 << " edgeI:" << startEdgeI << " verts:"
259 << mesh.edges()[startEdgeI]
260 << " to edge:" << edgeI << " verts:"
261 << mesh.edges()[edgeI]
262 << endl;
263 }
264
265 // Cross edge to all faces connected to it.
266
267 const labelList& eFaces = mesh.edgeFaces()[edgeI];
268
269 forAll(eFaces, eFacei)
270 {
271 label facei = eFaces[eFacei];
272
273 walkPointConnectedFaces
274 (
275 mesh,
276 regionEdges,
277 regionPointi,
278 facei,
279 edgeI,
280 visitedEdges
281 );
282 }
283 }
284 }
285}
286
287
288// Find all faces reachable from all non-fence points and staying on
289// regionFaces side.
290void Foam::regionSide::walkAllPointConnectedFaces
291(
292 const primitiveMesh& mesh,
293 const labelHashSet& regionFaces,
294 const labelHashSet& fencePoints
295)
296{
297 //
298 // Get all (internal and external) edges on region.
299 //
300 labelHashSet regionEdges(4*regionFaces.size());
301
302 for (const label facei : regionFaces)
303 {
304 const labelList& fEdges = mesh.faceEdges()[facei];
305
306 regionEdges.insert(fEdges);
307 }
308
309
310 //
311 // Visit all internal points on surface.
312 //
313
314 // Storage for visited points
315 labelHashSet visitedPoint(4*regionFaces.size());
316
317 // Insert fence points so we don't visit them
318 for (const label pointi : fencePoints)
319 {
320 visitedPoint.insert(pointi);
321 }
322
323 labelHashSet visitedEdges(2*fencePoints.size());
324
325
326 if (debug)
327 {
328 Info<< "Excluding visit of points:" << visitedPoint << endl;
329 }
330
331 for (const label facei : regionFaces)
332 {
333 // Get side of face.
334 label celli;
335
336 if (sideOwner_.found(facei))
337 {
338 celli = mesh.faceOwner()[facei];
339 }
340 else
341 {
342 celli = mesh.faceNeighbour()[facei];
343 }
344
345 // Find starting point and edge on face.
346 const labelList& fEdges = mesh.faceEdges()[facei];
347
348 forAll(fEdges, fEdgeI)
349 {
350 label edgeI = fEdges[fEdgeI];
351
352 // Get the face 'perpendicular' to facei on region.
353 label otherFacei = otherFace(mesh, celli, facei, edgeI);
354
355 // Edge
356 const edge& e = mesh.edges()[edgeI];
357
358 if (!visitedPoint.found(e.start()))
359 {
360 Info<< "Determining visibility from point " << e.start()
361 << endl;
362
363 visitedPoint.insert(e.start());
364
365 //edgeI = otherEdge(mesh, otherFacei, edgeI, e.start());
366
367 walkPointConnectedFaces
368 (
369 mesh,
370 regionEdges,
371 e.start(),
372 otherFacei,
373 edgeI,
374 visitedEdges
375 );
376 }
377 if (!visitedPoint.found(e.end()))
378 {
379 Info<< "Determining visibility from point " << e.end()
380 << endl;
381
382 visitedPoint.insert(e.end());
383
384 //edgeI = otherEdge(mesh, otherFacei, edgeI, e.end());
385
386 walkPointConnectedFaces
387 (
388 mesh,
389 regionEdges,
390 e.end(),
391 otherFacei,
392 edgeI,
393 visitedEdges
394 );
395 }
396 }
397 }
398}
399
400
401// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
402
403// Construct from components
405(
406 const primitiveMesh& mesh,
407 const labelHashSet& region, // faces of region
408 const labelHashSet& fenceEdges, // outside edges
409 const label startCelli,
410 const label startFacei
411)
412:
413 sideOwner_(region.size()),
414 insidePointFaces_(region.size())
415{
416 // Storage for visited faces
417 labelHashSet visitedFace(region.size());
418
419 // Visit all faces on this side of region.
420 // Mark for each face whether owner (or neighbour) has been visited.
421 // Sets sideOwner_
422 visitConnectedFaces
423 (
424 mesh,
425 region,
426 fenceEdges,
427 startCelli,
428 startFacei,
429 visitedFace
430 );
431
432 //
433 // Visit all internal points on region and mark faces visitable from these.
434 // Sets insidePointFaces_.
435 //
436
437 labelHashSet fencePoints(fenceEdges.size());
438
439 for (const label edgei : fenceEdges)
440 {
441 const edge& e = mesh.edges()[edgei];
442
443 fencePoints.insert(e.start());
444 fencePoints.insert(e.end());
445 }
446
447 walkAllPointConnectedFaces(mesh, region, fencePoints);
448}
449
450
451// ************************************************************************* //
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
Determines the 'side' for every face and connected to a singly-connected (through edges) region of fa...
Definition: regionSide.H:64
static label otherFace(const primitiveMesh &mesh, const label celli, const label excludeFacei, const label edgeI)
Step across edge onto other face on cell.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
label otherCell(const primitiveMesh &mesh, const label celli, const label facei)
Return cell on other side of face. Throws error.
Definition: meshTools.C:579
label otherEdge(const primitiveMesh &mesh, const labelList &edgeLabels, const label thisEdgeI, const label thisVertI)
Return label of other edge (out of candidates edgeLabels)
Definition: meshTools.C:521
label otherFace(const primitiveMesh &mesh, const label celli, const label facei, const label edgeI)
Return face on cell using edgeI but not facei. Throws error.
Definition: meshTools.C:555
void getEdgeFaces(const primitiveMesh &mesh, const label celli, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
Definition: meshTools.C:479
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
error FatalError
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333