orientFaceZone.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) 2013-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26Application
27 orientFaceZone
28
29Group
30 grpMeshManipulationUtilities
31
32Description
33 Corrects the orientation of faceZone.
34
35 - correct in parallel - excludes coupled faceZones from walk
36 - correct for non-manifold faceZones - restarts walk
37
38\*---------------------------------------------------------------------------*/
39
40#include "argList.H"
41#include "Time.H"
42#include "syncTools.H"
44#include "PatchEdgeFaceWave.H"
45#include "orientedSurface.H"
46#include "globalIndex.H"
47
48using namespace Foam;
49
50// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51
52int main(int argc, char *argv[])
53{
54 argList::addNote
55 (
56 "Corrects the orientation of faceZone"
57 );
58 #include "addRegionOption.H"
59 argList::addArgument("faceZone");
60 argList::addArgument("point", "A point outside of the mesh");
61
62 #include "setRootCase.H"
63 #include "createTime.H"
64 #include "createNamedPolyMesh.H"
65
66 const word zoneName = args[1];
67 const point outsidePoint = args.get<point>(2);
68
69 Info<< "Orienting faceZone " << zoneName
70 << " such that " << outsidePoint << " is outside"
71 << nl << endl;
72
73
74 const faceZone& fZone = mesh.faceZones()[zoneName];
75
76 if (fZone.checkParallelSync())
77 {
79 << "Face zone " << fZone.name()
80 << " is not parallel synchronised."
81 << " Any coupled face also needs its coupled version to be included"
82 << " and with opposite flipMap."
83 << exit(FatalError);
84 }
85
86 const labelList& faceLabels = fZone;
87
89 (
90 IndirectList<face>(mesh.faces(), faceLabels),
91 mesh.points()
92 );
93
94
95
96 const bitSet isMasterFace(syncTools::getMasterFaces(mesh));
97
98
99 // Data on all edges and faces
100 List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
101 List<patchFaceOrientation> allFaceInfo(patch.size());
102
103 // Make sure we don't walk through
104 // - slaves of coupled faces
105 // - non-manifold edges
106 {
107 const polyBoundaryMesh& bm = mesh.boundaryMesh();
108
109 label nProtected = 0;
110
111 forAll(faceLabels, facei)
112 {
113 const label meshFacei = faceLabels[facei];
114 const label patchi = bm.whichPatch(meshFacei);
115
116 if
117 (
118 patchi != -1
119 && bm[patchi].coupled()
120 && !isMasterFace[meshFacei]
121 )
122 {
123 // Slave side. Mark so doesn't get visited.
124 allFaceInfo[facei] = orientedSurface::NOFLIP;
125 nProtected++;
126 }
127 }
128
129 Info<< "Protected from visiting "
130 << returnReduce(nProtected, sumOp<label>())
131 << " slaves of coupled faces" << nl << endl;
132 }
133 {
134 // Number of (master)faces per edge
135 labelList nMasterFaces(patch.nEdges(), Zero);
136
137 forAll(faceLabels, facei)
138 {
139 const label meshFacei = faceLabels[facei];
140
141 if (isMasterFace[meshFacei])
142 {
143 const labelList& fEdges = patch.faceEdges()[facei];
144 forAll(fEdges, fEdgeI)
145 {
146 nMasterFaces[fEdges[fEdgeI]]++;
147 }
148 }
149 }
150
151 syncTools::syncEdgeList
152 (
153 mesh,
154 patch.meshEdges(mesh.edges(), mesh.pointEdges()),
155 nMasterFaces,
157 label(0)
158 );
159
160
161 label nProtected = 0;
162
163 forAll(nMasterFaces, edgeI)
164 {
165 if (nMasterFaces[edgeI] > 2)
166 {
167 allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
168 nProtected++;
169 }
170 }
171
172 Info<< "Protected from visiting "
173 << returnReduce(nProtected, sumOp<label>())
174 << " non-manifold edges" << nl << endl;
175 }
176
177
178
179 DynamicList<label> changedEdges;
181
182 const scalar tol = PatchEdgeFaceWave
183 <
186 >::propagationTol();
187
188 int dummyTrackData;
189
190 globalIndex globalFaces(patch.size());
191
192 while (true)
193 {
194 // Pick an unset face
195 label unsetFacei = labelMax;
196 forAll(allFaceInfo, facei)
197 {
198 if (allFaceInfo[facei] == orientedSurface::UNVISITED)
199 {
200 unsetFacei = globalFaces.toGlobal(facei);
201 break;
202 }
203 }
204
205 reduce(unsetFacei, minOp<label>());
206
207 if (unsetFacei == labelMax)
208 {
209 break;
210 }
211
212 label proci = globalFaces.whichProcID(unsetFacei);
213 label seedFacei = globalFaces.toLocal(proci, unsetFacei);
214 Info<< "Seeding from processor " << proci
215 << " face " << seedFacei << endl;
216
217 if (proci == Pstream::myProcNo())
218 {
219 // Determine orientation of seedFace
220
221 vector d = outsidePoint-patch.faceCentres()[seedFacei];
222 const vector& fn = patch.faceNormals()[seedFacei];
223
224 // Set information to correct orientation
225 patchFaceOrientation& faceInfo = allFaceInfo[seedFacei];
226 faceInfo = orientedSurface::NOFLIP;
227
228 if ((fn&d) < 0)
229 {
230 faceInfo.flip();
231
232 Pout<< "Face " << seedFacei << " at "
233 << patch.faceCentres()[seedFacei]
234 << " with normal " << fn
235 << " needs to be flipped." << endl;
236 }
237 else
238 {
239 Pout<< "Face " << seedFacei << " at "
240 << patch.faceCentres()[seedFacei]
241 << " with normal " << fn
242 << " points in positive direction (cos = " << (fn&d)/mag(d)
243 << ")" << endl;
244 }
245
246
247 const labelList& fEdges = patch.faceEdges()[seedFacei];
248 forAll(fEdges, fEdgeI)
249 {
250 label edgeI = fEdges[fEdgeI];
251
252 patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
253
254 if
255 (
256 edgeInfo.updateEdge<int>
257 (
258 mesh,
259 patch,
260 edgeI,
261 seedFacei,
262 faceInfo,
263 tol,
264 dummyTrackData
265 )
266 )
267 {
268 changedEdges.append(edgeI);
269 changedInfo.append(edgeInfo);
270 }
271 }
272 }
273
274
275 if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
276 {
277 break;
278 }
279
280
281
282 // Walk
284 <
287 > calc
288 (
289 mesh,
290 patch,
291 changedEdges,
292 changedInfo,
293 allEdgeInfo,
294 allFaceInfo,
295 returnReduce(patch.nEdges(), sumOp<label>())
296 );
297 }
298
299
300 // Push master zone info over to slave (since slave faces never visited)
301 {
302 const polyBoundaryMesh& bm = mesh.boundaryMesh();
303
304 labelList neiStatus(mesh.nBoundaryFaces(), orientedSurface::UNVISITED);
305
306 forAll(faceLabels, i)
307 {
308 const label meshFacei = faceLabels[i];
309 if (!mesh.isInternalFace(meshFacei))
310 {
311 neiStatus[meshFacei-mesh.nInternalFaces()] =
312 allFaceInfo[i].flipStatus();
313 }
314 }
315 syncTools::swapBoundaryFaceList(mesh, neiStatus);
316
317 forAll(faceLabels, i)
318 {
319 const label meshFacei = faceLabels[i];
320 const label patchi = bm.whichPatch(meshFacei);
321
322 if
323 (
324 patchi != -1
325 && bm[patchi].coupled()
326 && !isMasterFace[meshFacei]
327 )
328 {
329 // Slave side. Take flipped from neighbour
330 label bFacei = meshFacei-mesh.nInternalFaces();
331
332 if (neiStatus[bFacei] == orientedSurface::NOFLIP)
333 {
334 allFaceInfo[i] = orientedSurface::FLIP;
335 }
336 else if (neiStatus[bFacei] == orientedSurface::FLIP)
337 {
338 allFaceInfo[i] = orientedSurface::NOFLIP;
339 }
340 else
341 {
343 << "Incorrect status for face " << meshFacei
344 << abort(FatalError);
345 }
346 }
347 }
348 }
349
350
351 // Convert to flipmap and adapt faceZones
352
353 boolList newFlipMap(allFaceInfo.size(), false);
354 label nChanged = 0;
355 forAll(allFaceInfo, facei)
356 {
357 if (allFaceInfo[facei] == orientedSurface::NOFLIP)
358 {
359 newFlipMap[facei] = false;
360 }
361 else if (allFaceInfo[facei] == orientedSurface::FLIP)
362 {
363 newFlipMap[facei] = true;
364 }
365 else
366 {
368 << "Problem : unvisited face " << facei
369 << " centre:" << mesh.faceCentres()[faceLabels[facei]]
370 << abort(FatalError);
371 }
372
373 if (fZone.flipMap()[facei] != newFlipMap[facei])
374 {
375 nChanged++;
376 }
377 }
378
379 reduce(nChanged, sumOp<label>());
380 if (nChanged > 0)
381 {
382 Info<< "Flipping " << nChanged << " out of "
383 << globalFaces.size() << " faces." << nl << endl;
384
385 mesh.faceZones()[zoneName].resetAddressing(faceLabels, newFlipMap);
386 if (!mesh.faceZones().write())
387 {
389 << "Failed writing faceZones" << exit(FatalError);
390 }
391 }
392
393 Info<< "End\n" << endl;
394
395 return 0;
396}
397
398
399// ************************************************************************* //
reduce(hasMovingMesh, orOp< bool >())
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A List with indirect addressing.
Definition: IndirectList.H:119
Wave propagation of information along patch. Every iteration information goes through one layer of fa...
A list of faces which address into the list of points.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
virtual bool checkParallelSync(const bool report=false) const
Check whether all procs have faces synchronised.
Definition: faceZone.C:516
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Transport of orientation for use in PatchEdgeFaceWave.
void flip()
Reverse the orientation.
bool updateEdge(const polyMesh &mesh, const indirectPrimitivePatch &patch, const label edgeI, const label facei, const patchFaceOrientation &faceInfo, const scalar tol, TrackingData &td)
Influence of face on edge.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
const word & name() const noexcept
The zone name.
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
dynamicFvMesh & mesh
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const std::string patch
OpenFOAM patch number as a std::string.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
constexpr label labelMax
Definition: label.H:61
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333